Mesh Mess

FEM, CFD, CAE, and all that.

Tag: potential flow

Energy dissipation of water waves

This is a interesting, and subtle argument, pointed out to me by Harry Yeh. For a linear wave, energy dissipation rate, averaged over length and time time, is

2 \mu \frac{1}{LT} \int_0^T\int_0^L\int_{-h}^0{\epsilon_{ij}}\,dydxdt=2 \mu \frac{1}{LT} \int_0^T\int_0^L\int_{-h}^0{\frac{\partial^2 \phi}{\partial x\partial y}^2}\,dydxdt=2 \mu k \frac{a^2 \omega^2}{\tanh k h}

and it’s apparently non-zero (\epsilon_{ij} is the strain rate tensor). However, for mechanical energy E(t), we have

E(t)=\int_{D(t)}{\frac{1}{2}\rho |\nabla \phi|^2+\rho g y}\,dV

and its material derivative is

\frac{D E}{D t}=\int_D{\rho \nabla \phi\cdot\nabla\phi_t}\,dV+\int_{\partial D}{(\frac{1}{2}\rho|\nabla\phi|^2+\rho g y)U\cdot n}\,ds

Using \nabla\cdot(\phi_t\nabla\phi)=\nabla\phi\cdot\nabla\phi_t+\phi_t\nabla^2\phi=\nabla\phi\cdot\nabla\phi_t and divergence theorem, we have

\frac{D E}{D t}=\int_{\partial D}{\rho \phi_t(\nabla\phi\cdot n-U\cdot n)-p U\cdot n}\,ds

At impermeable bottom, above integral is zero due to \phi \cdot n=0 and U\cdot n=0; at free surface, kinematic BC and dynamic BC (p=0) indicate above integral is zero as well. Therefore

DE/Dt=0

Question: how to explain the seemingly inconsistency of above two results ? i.e. why we have constant mechanical energy, in the meantime a nonzero dissipation rate?

In order to do this, we use the general form of time derivative of mechanical energy E:

\frac{D}{D t}(\rho \frac{|u|^2}{2}+\rho g y)=\nabla\cdot(u\cdot \tau_{ij})-\tau_{ij}:\nabla u

the 1st term RHS is rate of work done by surroundings to the fluid (material domain), while the 2nd represents the dissipation of mechanical energy. We can see RHS is zero after using divergence theorem and BC’s. But, the key point here is that those two terms are nonzero by themselves, unless flow is inviscid.The last thing to be emphasized is that, the linear momentum for potential flow will have no dissipation term even for viscous flows, due to incompressibility condition.

Now we can see, for a viscous, potential flow, the dissipation rate is non-zero, but its mechanical energy remains unchanged because of the work done at the free surface. This work done, as energy input, is in the form of viscous stress (from air), not covered by DBC p=0. If free surface BC is replaced with non-stress condition (different from non-pressure condition), the flow will constantly lose energy due to net dissipation.

Mixed-Eularian-Lagrangian method

In last post I talked wrote about a time integration for NS equation. As a special case, time update for potential flow is actually easier, due to fact that

1. the governing equation, i.e. the Laplace equation is linear.

2. the solution dependes only on the boundary condition, in which only the free surface boundary condition is the tricky part.

3. Even though the field equation is in Eularian representation, the free surface can be described using Lagrangian representation.

For free surface potential flow, the time updating of the problem is twofold: update the (free surface) geometry and update the boundary conditioin at the free surface. Without mentioning the boundary condition at other part of the domain boundary, we are solving the problem characterized by

Field equation: \nabla^2 \phi=0,\quad \forall x\in \sigma

Free surface kinetic boundary condition(FSKBC): D\vec{r}/Dt=\nabla \phi, \quad \forall x\in \gamma_0

Free surface dynamic boundary condition(FSDBC): \phi_t+1/2 (\nabla\phi)^2+gz+p/\rho=0,\quad \forall x\in \gamma_0

The first boundary condition says that the particle velocity for material particle on the free surface, i.e. the velocity of free surface, equals to the flow velocity defined by velocity potential. The second boundary condition is theĀ  Bernoulli equation, with constant on RHS assigned to 0, which could be looked as consumed by \phi_t. Remembering that

\phi_t=D\phi/Dt-\vec{u}\cdot \nabla \phi=D\phi/Dt-(\nabla \phi)^2

FSDBC can be converted to

D\phi/Dt=1/2(\nabla \phi)^2-g\eta-p/\rho,\quad \forall x\in \gamma_0

The plan of mixed-Eularian-Lagrangian(MEL) method, is to solve \phi^n at time t_n, then use the BC equations to update the boundary condition. One should remember, the LHS of BC equations use material particle as variable, i.e. in Lagrangian representation. Denote particle on the free surface by \zeta, MEL is characterized by two steps:

1. Solve field equation with BCs with variable at t_n, i.e. \phi^n, \vec{r}^n;

2. Update the BCs at the free surface using

(\vec{r}^{n+1}(\zeta)-\vec{r}^n(\zeta))/\Delta t=\nabla \phi^n(\zeta)

(\phi^{n+1}(\zeta)-\phi^n(\zeta))/\Delta t=1/2(\nabla \phi^n(\zeta))^2-g\eta-p/\rho

As mentioned at the end of last post, usually the computing cost is characterized by the first step, while the second step concerns the stability and convergence.