For a test case, I want to determine the velocity profile of a viscously damped standing wave.

By linearizing the density ($\rho=\rho_0+\rho'$) and velocity ($ux=ux'$), the continuity and Navier-Stokes equations result in, respectively:

\begin{align} \partial_t\rho' + \rho_0\partial_xu_x' &= 0 \tag{1} \\ \partial_t^2\rho' &= \partial_x^2\rho'c_s^2 + \nu\partial_t\partial_x\rho' \tag{2} \end{align}

The $c_s$ is just a constant indicating we are dealing with an ideal pressure term ($p=\rho c_s^2$)

A solution for the density to $(2)$ is given by:


where $$k_x=2\pi/n_x, \quad \omega_r=\frac{1}{2}k_x^2\nu, \quad \omega_i=k_xc_s\sqrt{1-\left(\frac{1}{2}\frac{k_x\nu}{c_s} \right)^2} \, .$$

Now I want to determine the velocity; it would seem straightforward to use $(1)$ to get


and integrate to get


where $K$ is an integration constant. My approach was to determine $K$ by setting the velocity zero at a antinode (at $x=n_x/4$), to get

$$u_{x}'=-\frac{1}{k_{x}}\frac{\triangle\rho}{\rho_{0}}\cos\left(k_{x}x\right)\left[\omega_{r}\cos\left(\omega_{i}t\right)-\omega_{i}\sin\left(\omega_{i}t\right)\right]\exp\left(-\omega_{r}t\right) \, .$$

However, comparing the simulation with the analytical solution it seems that the amplitude of the velocity is much larger in the simulation.

Is my approach described above at all correct?

  • $\begingroup$ No one has any ideas or tips? $\endgroup$
    – nluigi
    Aug 3, 2015 at 10:55
  • 1
    $\begingroup$ Just some thoughts. If $u'_x$ is a perturbation, why would $K$ be non-zero? Do I understand correctly that you have periodic conditions? I assume you did not linearize your simulation, so how convinced are you that the linearization is justified? Is amplitude really the only difference between your analytical and numerical analysis (i.e, normalized you have agreement?) $\endgroup$
    – Bernhard
    Aug 4, 2015 at 15:52
  • $\begingroup$ @Bernhard, you are right that $K$ is zero... i hadn't realized $\cos(k_xn_x/4)=0$. I indeed have periodic boundary conditions. In my simulation i enforce $\Delta\rho/\rho0\ll1$, which i assume is justification enough to linearize the model. $\endgroup$
    – nluigi
    Aug 6, 2015 at 14:35
  • $\begingroup$ To get the solution to equation 2 are you just assuming $\partial_{t} \rightarrow -i \omega$ and $\partial_{t} \rightarrow i k_{x}$? Also, how do you define $\Delta \rho$? Is it the difference between the perturbed and unperturbed or the total and the unperturbed? $\endgroup$ Aug 8, 2015 at 20:09
  • $\begingroup$ I assume $\rho'$ is separable into $\rho'=f(x)g(t)$ to get $\frac{1}{f+\partial_tf}\partial_t^2f=\frac{1}{g}\partial_x^2g=-k^2$. For $g(x)$ it follows: $g(x)=A\sin(k_xx)+B\cos(k_xx)$. For $f(t)$ i assume a form $f(t)=\exp(-\omega t)$, which leads to a quadratic equation for $\omega$: $\omega^2-k^2(\omega-1)=0$. Since it is a standing wave, I can decompose $\omega$ into: $\omega=\omega_r+i\omega_i$ and find that $omega_r=\frac{1}{2}k^2$ and $\omega_i=\pm\sqrt{1-(\frac{1}{2}k)^2}$. Since $\omega$ is complex, it follows that: $f(t)=\cos(\omega_it)\exp(-\omega_rt)$ $\endgroup$
    – nluigi
    Aug 11, 2015 at 10:30

2 Answers 2


Make sure you've normalized everything correctly in both the analytical and numerical solutions so you're comparing apples to apples. Is $n_x$ the wavelength? If so, then the factor of $\cos\left(k_x\frac{n_x}{4} \right)$ is just 0. That seems right, since $u'_x$ is then $\pi/2$ out of phase with $\Delta \rho$, and the velocity perturbation is symmetrical. Try setting $\omega_r = 0$ for a more transparent solution. Otherwise, everything looks fine with the analytical solution.

  • $\begingroup$ When normalized, the analytical solution for the density looks the same except after some time the amplitude decreases quicker for the simulation. I will post the results in abit to make it more clear. Unfortunately due to the nature of the simulation i cannot set $\omega_r=0$ without also setting $\omega_i=0$. $\endgroup$
    – nluigi
    Aug 6, 2015 at 16:19
  • $\begingroup$ It seems unlikely to me that the time dependence would be anything other than the $exp(-\omega_r t)$ you have in the analytical solution, so I would look for a bug in the code. Also check the long-time behavior - nonlinear effects should become minimal later on even if they change the time dynamics early. Why can't you set $\nu = 0$? Alternatively, can you try eliminating the nonlinear term in the code and see if you get better agreement? $\endgroup$
    – pwf
    Aug 6, 2015 at 17:31
  • $\begingroup$ It seems unlikely to me aswell, the only thing it can be is the derived form of $\omega_r$ is incorrect. I am quite sure, non-linear effects (u mean due to $u\partial_x u$ right?) are absent. I am actually not solving the continuity and Navier-Stokes equations using standard CFD tools but using Lattice Boltzmann methods. Inherently, I cannot set $\nu=0$ because in that limit the dynamics are undefined. $\endgroup$
    – nluigi
    Aug 11, 2015 at 10:22

I would go about the theoretical discussion differently.

First consider irrotational inviscid waves, governed by Laplace's equation in the interior, i.e. $$\nabla^2 \phi = 0$$ where $\phi=\phi(x,z,t)$ is the velocity potential, $z$ is the vertical and $x$ the horizontal direction.

The boundary conditions for water waves, in the absence of viscosity, are

$$ \phi_t+\frac{1}{2}(\nabla\phi)^2 +gz =0;\quad \eta_t+\phi_x\eta_x = \phi_z \quad @z=\eta;$$

where $\eta$ is the free surface displacement, and both boundary conditions are evaluated at the free surface. Finally, we have $\phi_z=0$ at $z=-h$, with $h$ the depth of the water, taken to tend towards infinity. Note, the governing equation is linear, but the boundary conditions are nonlinear, and more strikingly, are evaluated at a dependent variable of the system. The later of which is the reason why these equations are very difficult to solve.

Now, a linear standing wave can be thought of as two oppositely traveling waves of equal amplitude, frequency. Take $\omega$ to be the angular frequency, and $k$ the wavenumber, then it is easy to show that the linearized boundary conditions together with Laplace's equation implies

$$\phi = \frac{a\omega}{k} \cos kx \sin\omega t\ e^{kz}; \quad \eta=a\cos kx \cos \omega t$$

where $\omega=\sqrt{gk}$.

Now, to add viscosity.

The motion described above can exist even with viscosity if we apply the following normal and tangential surface stresses:

$$\vec{F}\cdot \hat{z} =-p +2\mu \frac{\partial v}{\partial y} = -p+2\mu k^2\frac{a\omega}{k}\cos kx \sin \omega t\ e^{kz}$$ and $$\vec{F}\cdot \hat{x} = \mu\left(\frac{\partial v}{\partial x}+\frac{\partial u}{\partial y}\right)= -2\mu k^2 a\frac{\omega}{k} \sin kx \sin \omega t\ e^{kz}$$

with $\mu$ the dynamic viscosity and $(u,v)=(\phi_x,\phi_z)$. Therefore, the work done by these forces, averaged over a wavelength and over a wave period

$$\frac{1}{T}\frac{1}{\lambda} \int_0^T\int_0^{\lambda} \vec{u}\cdot \vec{F} dx =\mu k^2 a^2, $$

where we assume that the wave amplitude varies slowly compared with the frequency of the wave.

Next, the total energy in a standing wave is $E=\frac{1}{2} \rho g a^2$, so that in the absence of surface forces, we must have

$$\frac{d}{dt}E = 2\mu k^2 a^2 \omega \implies a=a_o^{-\nu k^2t}.$$

We note that kinematic viscosity $\nu =\mu/\rho$ is much more effective at annihilating shorter waves, as we would naively expect from the form of the dissipation given in the Navier Stokes equations.

The theoretical free surface profile and velocity potential are

$$\eta= a_o^{-\nu k^2t}\cos kx \cos \omega ; \quad \phi =\frac{a_o^{-\nu k^2t} \omega}{k}\cos kx \sin \omega .$$

and finally the velocity fields are

$$(u,v) = a_o^{-\nu k^2t} \omega\sin \omega t \ e^{kz} (-\sin kx, cos kx).$$

Reference: Lamb (1932, $\S$ 348)

  • $\begingroup$ Thank you for the reference and i like your solution to two-dimensional standing surface waves! However, while it is about standing waves, it is only remotely related to my question which deals with a one-dimensional viscously damped standing wave based on the linearized continuity and navier-stokes equations. In my case there is also now gravity term present. $\endgroup$
    – nluigi
    Aug 11, 2015 at 10:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.