As an appendix promised for the article How fast helium balloon floats…, I’m going to explain the hint on how to derive the insights behind the quadratic drag force presents in the damped harmonic oscillations found in the previous articles.

Solving the usual damped oscillations is simple enough because it’s linear. Meanwhile, the quadratic force is not so straight-forward. Originally, I didn’t find any open literature that discuss this. There are several articles found online with the keyword “Quadratic Drag Force Damping Oscillator”, but most of it is behind a paygate, so I can’t access it. Some of the article I found, only discuss a falling object under drag force. That one is also simple, but way different from the case of quadratic damping found in the oscillator.

So, at that time, I have no choice but to try to solve it myself. Perhaps there is a better way? In that case, let me know.

Looking at the equation of motions

The equation of motions can be found from a simple Newton’s Law. But there is a tricky catch.

We first introduce an oscillator with position $x(t)$ measured from it’s equilibrium position.

This is a damping oscillation with all linear order of position $x(t)$:

$\ddot x=-\omega^2 x - D_1\dot x$

Meanwhile, this is a quadratic damping of the same oscillation:

$\ddot x=-\omega^2 x - D_2\dot{x}^2$

All the second terms corresponds to the drag force, with each drag constant ($D_1$ and $D_2$) so that we remember that the dimension is different between each equations. However, we should be careful because in the case of linear damping, the information about the direction of the drag force is conveyed entirely in the expression. Meanwhile for the quadratic damping above, it is only valid if the direction of the movement is $\dot{x}\ge0$ or into the positive axis.

It is not complete. We should also observe the equation when the velocity is pointing to the opposite direction.

$\ddot x= -\omega^2 x + D_2 \dot{x}^2$

As you can see, the sign alternates. More compactly it can be written:

$\ddot x = -\omega^2 x - D_2 \dot{x} \left| \dot{x} \right|$

But we still need to solve both separately.

Separating the variables

Both quadratic damping equation is difficult to solve directly.

Usually, from the equation of motion we convert it into an energy equation like this:

$\ddot x=-\omega^2 x - D_2\dot{x}^2 \\ \dot x \cdot d\dot x=-\omega^2 x \cdot dx - D_2 \dot{x}^2 dx$

But the second term in the right hand side is difficult to integrate. Even if we can, we still need to separate the variables to further integrate it to retrieve the time function of x $t(x)$.

Instead of that, I offer an alternative approach. We know that the derivative of $\dot{x}^2$ will be somewhat linear in power (not quadratic). So instead of integrating the equation of motions, we took a derivative instead.

$\frac{d\ddot{x}}{dt} = -\omega^2 \dot x - 2D_2 \dot{x} \ddot{x}$

I couldn’t find a way to express triple dot in Katex, so the above notation will do. Notice that the right hand side has common $\dot{x}$ and it could also be “chain-rule”’d so that we eliminate $dt$. Here’s the differential form.

$d\ddot{x}=(-\omega^2-2D_2\ddot{x})dx$

The above form is integrable and separable if we move the right hand terms with $\ddot{x}$ to left hand side. The above expression is for $\dot{x}>0$. For $\dot{x}<0$, it’s similar with following the same reasoning.

$d\ddot{x}=(-\omega^2+2D_2\ddot{x})dx$

Inferring the phase space diagram

From the previous differential form, by integration we have the following relations.

$\dot{x} > 0; k_+=e^{2D_2x}(\omega^2+2D_2\ddot{x}) \\ \dot{x} < 0; k_-=e^{-2D_2x}(\omega^2-2D_2\ddot{x})$

I put the constant in the left hand side intentionally to emphasize more that the equation looks like a contour diagram. For example in the equation of circle. But, in this case, the amplitude of the contour (distance from the origin), is controlled by the ratio $e^{2D_2x}$.

A constant $k_+$ is used when the object oscillates with positive velocity, up to a half cycle. After that the velocity is zero and it will return back with negative velocity, in which we use $k_-$ for the remaining half cycle. To get these constants, we can use initial condition. But since each initial condition is calculated from the start of the half cycles, we can infer that each constants will have different values and probably have a pattern for even-half-cycles and odd-half cycles.

We can retrieve relations between (only) $x$ and $\dot{x}$ by substituting it back into the previous equation of motions.

$\dot{x} > 0; k_+=e^{2D_2x}(\omega^2 - 2D_2 \omega^2 x-2 D_2^2 \dot{x}^2)\\ \dot{x} < 0; k_-=e^{-2D_2x}(\omega^2 +2D_2 \omega^2 x-2 D_2^2 \dot{x}^2)$

Because we have relations between $x$ and $\dot{x}$, we might understand more if we create a phase space diagram. By drawing a phase diagram, you can understand more about the cycles of the oscillator.

If we treat $x$ as the horizontal axis, and then $\dot{x}$ as the vertical axis, then we can set that the oscillations happens when maximum displacement starts. Let’s start with negative x.

Before we proceed, I just wanted to note that I should’ve made a good diagram, but I haven’t found a good math graph library yet for React + MDX, so maybe I’ll do that later. Anyway, it should be straightforward to draw the phase diagram if you set some initial conditions.

We can avoid setting initial conditions by using unitary quantity. Let’s just introduce a new dimensionless variables for both $x$ as $X$ and $\dot{x}$ as $V$. The new variables have maximum absolute value of one. Because of this, the constant will also change unit/dimension to match, so we renamed it as $K_+$ and $K_-$.

We substitute this:

$X=2D_2x \\ V^2=\frac{2D_2^2}{\omega^2}\dot{x}^2 \\ \alpha =\frac{2D_2}{\omega^2}\ddot{x} \\ K_+=\frac{k_+}{\omega^2} \\ K_-=\frac{k_-}{\omega^2}$

The phase space relations conveniently becomes:

$V > 0; K_+=e^X(1+\alpha) \\ V < 0; K_-=e^{-X}(1-\alpha) \\ V > 0; K_+=e^X(1-X-V^2) \\ V < 0; K_-=e^{-X}(1+X-V^2)$

The relations is now more simple, although we need to remember the corresponding signs for each half cycles.

So, to proceed with understanding the cycles. We first start with initial position $X=-1$. The velocity here must be 0, so that means $K_+=\frac{2}{e}$.

By knowing $K_+$, we can proceed to know what the maximum acceleration at the initial position. The maximum acceleration is $\frac{\omega^2}{2D_2}$. We infer that for each initial half cycles, the initial acceleration magnitude will be less than that. I have to say that I’m surprised the value entirely determined by $\omega^2$ and $D_2$, since we set that the maximum initial displacement in the phase space to be $\left|X\right|=1$.

Next, in the quarter cycle, the acceleration must be $\alpha=0$, so that velocity is at maximum. We got $X=\ln\left|\frac{2}{e}\right|$ . Then $\dot{x}=\frac{\omega}{D_2}\sqrt{\frac{1}{2}\ln{\frac{e}{2}}}$.

Continuing to the end at the half cycle. This one is a little bit tricky. We have to know when the $K_+$ relation intersects the X axis again. From the relation formula, we intuitively know that it must cross X axis twice at maximum. So you had to solve $K_+e^{-X}+X-1=0$, probably numerically. From this, you will have the next $X$ for using the next half cycle phase. Let’s just call our initial $X$ for the previous cycle to be $X_0$ and $X_1$ for the next cycle.

After this, you just have to repeat the process. This time, using $X_1$ to get the value for $K_-$. Once you got both $K_+$ and $K_-$, you can reuse the constant for every corresponding half cycle. So what I imagine the phase diagram will be is a similar shape for the corresponding half cycle with the same parity (odd with odd, even with even), but very different between odd and even half cycles.

Using a different set of coordinate origin

In the case of my original problem, the floating balloon. I am using a different reference coordinate. Instead of using the final equilibrium position as coordinate $x=0$, I’m using the initial displacement as $x=0$. This can tremendously help to guess the trajectory.

Let’s say that the distance between the initial position $x=0$ and the equilibrium is it’s maximum displacement $A$. Then by using Newton’s Law, when the damping is finished, whatever form of $\ddot{y}(t)$ it has to corresponds with the Newton’s Law.

Equation of motion without drag force:

$\ddot y =-\omega^2 (y - A)$

Equation of motion with drag force:

$\ddot y = -\omega^2 (y-A) -D_2 \dot{y}^2$

Note that, due to the origin shift froom coordinate $x$ to $y$, the original equation $\ddot x=-\omega^2 x$ were translated according to the $x=y-A$. It’s a simple translation so the derivative stays the same.

However, because of these coordinate shift it’s easier to imagine that the end result is the same, regardless if there is a drag force or not. From the Newton’s First Law, we can immediately calculate $A$. Instead of relying from the equation to derive initial condition.

You can also see that if the equation we got from Newton’s Law is non-homogenous (it has constant term), we can convert it into a homogenous equation with just coordinate shift. For example, one might start with equation:

$\ddot y = C-\omega^2 y -D_2 \dot{y}^2$

We can translate the origin point such that $C=\omega^2 A$, and convert the equation into:

$x=y-A \\ \ddot x = -\omega^2 x - D_2 \dot{x}^2$

As you can see, just by doing coordinate transform, we are back to the initial equation we discuss above.

Solving the trajectory of the motion

From the phase space relations between $X$ and $V$ we can solve for $t$ because the form is integrable and separable.

$\dot{X}=\frac{V}{\sqrt{2}\omega} \\ V>0; \omega \Delta t=\int{\frac{1}{\sqrt{2(1-X-K_+e^{-X})}}}dX \\ V<0; \omega \Delta t=\int{\frac{1}{\sqrt{2(1+X-K_-e^{X})}}}dX$

The left hand side is analog to an angular phase from a given angular speed. Just like in the case of harmonic oscillator. For each upward/downward equation, $\Delta t$ corresponds to the time it took to reach half of that cycle.

I don’t have a way to integrate it analytically, but it should be possible to integrate it numerically. The key insight that we have is the half-cycle period is some sort of function of $K$.

To get the position as a function of time (the trajectory of $x(t)$). You can invert the $t(x)$ function of each half cycle.

I personally think that the integral terms here: $\int{\frac{1}{\sqrt{2(1-X-K_+e^{-X})}}}dX$ have a very generic form that perhaps someone should invent a new function that have that terms as it derivative. Kind of like the Jacobi elliptic function that is not an elementary function but people invent it because it is related with orbit of celestial bodies.

Conclusions

In this article we discussed some hints to solve for quadratic dampening on harmonic oscillator. Not every question is answered, but it should provide enough insights to the dynamic of the system. The phase space that the system generates is very unique to this type of problem.