water wave | Mesh Mess

05/29/2012

You get what you mesh for

In finite element modeling, in fact, in all the numerical modeling, the resolution of the results depends on the space where the problem is projected. As to PDE’s numerical solution, the projection is from space with infinite dimension to space with finite dimension. So no matter it’s finite element or finite volume or finite difference, the size of the mesh is what you can see at the finest scale. Beyond that, it’s either physical modeling (such as turbulence subscale modeling) or some other educated guess.

Here is an example.

Water surface as wave coming in

This plot shows a finite element solution of Navier-Stokes equation, in a water waves problem. The free surface is tracked using Level Set method. Taking closer look reveals the following not-so-smooth “water surface”:

Free surface, enlarged

The resolution of the free surface is restricted by the size of the element. For any mesh, examing close enough reveals similar wiggles. To put it another way, in the numerical simulation, free surface is not really a surface, it’s more like a layer with finite thickness. If I plot the centroid of the elements that reside on the free surface, it looks like this:

Elements on the free surface

This is what a non-polished numerical result should look like. Smoothing such as with cubic splines sweepts everything under the table, including the message: “you get what you mesh for”.

Leave a comment

07/26/2010

A test station in the Caspian Sea

English Russia has this post about how a naval weapon test station was constructed in the Caspian Sea. Though impaired by the wavy sea (the wave height during a storm in Caspian sea could reach 7.5-8.0 m) and abandoned, the structure remains untoppled after over half a century. The construction of the foundation, is quite a engineering masterpiece even for today’s standard:

The entire underwater part of the gigantic construction, called “Massive”, was built on the shore in the foundation pit with a capacity of 530,000 cubic meters. Sometimes they even dug by hand, but more often dredgers were used.

On the bottom of the foundation pit a huge ferroconcrete “box” was built, measuring 14 meters high. Next, the plate separating the pit from the sea was destroyed and the surfaced “box” was dragged to the building site of the station. There, in 1935, the “box” was filled with water and placed on the specially site, which was made of quarry stone. Thereafter construction of the station was done in a usual way, with people and materials being delivered by ships.

Leave a comment

06/02/2010

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.

Leave a comment

04/07/2010

Oscillatory boundary layer

This was addressed by Lamb, and G.K. Batchelor, and many others. For shallow water wave, linear theory predicts a uniform oscillatory movement of water particles along depth. Assuming a impermeable plane ocean bottom, the boundary layer behavior is similar to that in the fluid with oscillating bottom ( $y=0$ ) and rest at infinity ( $y=\infty$ ). Solution of a 2D case shows a wave-like behavior for horizontal velocity varying in vertical direction within boundary layer, with corresponding wave number $\kappa=\sqrt{\omega/2\nu}$ , where $\omega$ is frequency of outer flow, and $\nu$ is kinematic viscosity. Solution also shows a exponential decay of this wave’s amplitude as $y$ increases. Dimension argument gives the effective range of viscosity at bottom for this kind of flow, for common ocean wave period, as of order of 1cm.

Horizontal velocity agains time

Horizontal velocity profile within quarter of outer flow period

Leave a comment

02/19/2010

Nonlinear shallow water equation

The name of the equation pretty much gives the quantities that matter in this scenario: the nonlinearity and shallowness. For the standard input of the water wave equation with velocity potential as the unknown, we have:

$\nabla^2 \phi=0$ as governing equation;

$\eta_t+\phi_x\eta_x+\phi_y\eta_y-\phi_z=0$ as kinematic boundary condition (KBC) at free surface;

$\phi_t+g\eta+\frac{1}{2}(\nabla \phi)^2+\frac{p}{\rho}=0$ as dynamics boundary condition (DBC) at free surface.

In order to quantify the effect of nonlinearity and shallowness, non-dimensional form of the equation is introduced by using normalized variables:

$(x, y)\leftarrow (x, y)/l,\quad z\leftarrow z/h, \quad t\leftarrow ct/l, \quad \eta\leftarrow \eta/a,\quad \phi\leftarrow h\phi/(alc)$

where $l$ is characteristic length scale, $a$ is wave amplitude, $h$ is water depth, $c$ is shallow water celerity by linear theory $c=\sqrt{gh}$ . The quantity that evaluates nonlinearity is defined as

$\epsilon=a/h$

and that for shallowness is

$\delta=(h/l)^2$

Now the equations become

$\Delta(\phi_{xx}+\phi_{yy})+\phi_{zz}=0$

$\phi_t+\frac{\epsilon}{2}(\phi_x^2+\phi_y^2)+\frac{\epsilon}{2\delta}\phi_z^2+\eta$ , for $z=1+\epsilon\eta$ (DBC)

$\delta[\eta_t+\epsilon(\phi_x\eta_x+\phi_y\eta_y)]-\phi_z=0$ , for $z=1+\epsilon\eta$ (KBC)

If the unknown $\phi$ is expanded with $\delta$ :

$\phi=\phi_0+\delta\phi_1+\delta^2\phi_2+\dots$

One can find for zero order term $\phi_0$ , governing equations gives:

$\phi_{0zz}=0$

combining homogeneous Neumann boundary condition at the bottom, this means $\phi_{0z}=0$ all over the region, and $\phi_0$ is independent of vertical location. If denote

$u(x,y,t)=\phi_{0x}, \quad v(x,y,t)=\phi_{0y}$

we have for $\phi_1$ :

$\phi_1=-z^2/2(u_x+v_y),\quad \phi_2=z^4/24[(\nabla^2u)_x+(\nabla^2v)_y]$

All above results come from the governing equation and bottom boundary condition, and have nothing to do with free surface boundary conditions, i.e. they are concluded by the fact of flow being incompressible, irrotational, and having immiscible bottom. Now, if we put the expansion of $\phi$ in free surface boundary condition equations, and retain only up to 2nd order of $\delta$ , $\epsilon$ in KBC, and up to 1st order of $\delta$ , $\epsilon$ in DBC, combing the definition notation of $(u,v)$ , we have:

$u_t+\epsilon(uu_x+vv_x)+\eta_x-1/2\delta(u_{txx}+v_{txy})=0$

$v_t+\epsilon(uu_y+vv_y)+\eta_y-1/2\delta(u_{txy}+v_{tyy})=0$

$\eta_t+[u(1+\epsilon\eta)]_x+[v(1+\epsilon\eta)]_y=\delta/6[(\nabla^2u)_x+(\nabla^2v)y]$

In above equations if we take $\delta\rightarrow 0$ and $\epsilon\rightarrow 0$ , i.e., assuming linearity and shallowness, we have the linear wave equation for $\eta$ :