Mesh Mess

FEM, CFD, CAE, and all that.

Tag: PDE

How to categorize different numerical methods for PDE?

All the numerical methods for PDE can be boiled down two things: the projection from infinite dimension to finite dimension where the numerical solution is seeked, and the discretization of the equaiton itself. Quote from page 1 of Hesthaven and Warburton’s book “Nodal Discontinuous Galerkin Methods“:

The construction of any numerical method for solving a partial differential equation requires one to consider the two choices:

  • How does one represent the solution by an approximate solution?
  • In which sense will the approximate solution satisfy the partial differential equation?

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”.

What’s the difference?

Today I was asked for a one-sentence comment of the difference between finite difference method and finite element method. What I came up was:

In finite difference method we approximate the differential expression, or the partial derivative, while in finite element method we approximate the space in which the solution lies.

Hopefully this shed some light on the issue.

Saddle point problem in Navier-Stokes equation

It is well known that the the optimization problem

\min J(x)= \frac{1}{2}(Ax, x)-(f, x)

is equal to solving a linear operator equation

Ax=f

provided that A is symmetric and positive definite. Here x\in X which is a Hilbert space, and f\in X' while A is a bilinear operator.

This is a very general problem considered widely in many areas. However, in the real world, some further condition are to be applied on the space where x lies on. If this is described as constraint equation

Bx=g

then we are dealing with constraint optimization problem ,to which Lagrange multiplier is usually the first one can come up. We then try to solve the stationary point of Lagrangian

L(x,\lambda)=\frac{1}{2}(Ax, x)-(f, x)+\lambda (Bx-g)

Now supposing we set X=R^n, using stationary point condition on above equation, we arrive at

\left[\begin{array}{cc} A & B^T\\ B & O \end{array}\right]\left[\begin{array}{c} x\\ y\end{array}\right]=\left[\begin{array}{c} f\\ g\end{array}\right]

where \lambda is replaced by y indicating the parallel unknown as x in the problem. The stationary point nature of this solution gives the name of saddle point problem, whose general case is

L(x,v)\le L(x,y) \le L(u,y), \forall u,v \Leftrightarrow \min_u\max_v L(u,v)=L(x,y)=\max_v\min_u L(u,v)

So, putting a constraint to a minimization problem, one try to solve a saddle problem instead of a simple linear system with SPD matrix.

For a linear system, one can always decompose it in to the form of

\left[\begin{array}{cc} A & B_1^T\\ B_2 & -C \end{array}\right]\left[\begin{array}{c} x\\ y\end{array}\right]=\left[\begin{array}{c} f\\ g\end{array}\right]

and this is the most general form of saddle problem.

For incompressible Navier-Stokes equation, we have (without mentioning the boundary conditions)

\frac{\partial\mathbf{u}}{\partial t}+\mathbf{u}\cdot\nabla\mathbf{u}-\nu\Delta \mathbf{u}+\nabla p/\rho=\mathbf{f}

\nabla\cdot\mathbf{u}=0

Following \theta-scheme by Glowinski, one has the time discretized plan as:

(\mathbf{u^{n+\theta}}-\mathbf{u}^n)/(\theta\Delta t)-\alpha\nu\Delta\mathbf{u}^{n+\theta}+\nabla p^{n+\theta}=\mathbf{f}^{n+\theta}+\beta\nu\Delta \mathbf{u}^n-(\mathbf{u}^n\cdot\Delta)\mathbf{u}^n

\nabla\cdot\mathbf{u}^{n+\theta}=0

then

(\mathbf{u^{n+1-\theta}}-\mathbf{u}^{n+\theta})/((1-2\theta)\Delta t)-\beta\nu\Delta\mathbf{u}^{n+1-\theta}+(\mathbf{u}^{n+1-\theta}\cdot\Delta)\mathbf{u}^{n+1-\theta}=\mathbf{f}^{n+\theta}+\alpha\nu\Delta \mathbf{u}^{n+\theta}-\nabla p^{n+\theta}

then

(\mathbf{u}^{n+\theta}-\mathbf{u}^{n+1-\theta})/(\theta\Delta t)-\alpha\nu\Delta\mathbf{u}^{n+1}+\nabla p^{n+1}=\mathbf{f}^{n+1}+\beta\nu\Delta \mathbf{u}^{n+1-\theta}-(\mathbf{u}^{n+1-\theta}\cdot\Delta)\mathbf{u}^{n+1-\theta}

\nabla\cdot\mathbf{u}^{n+1}=0

The first and last step in above scheme are to solve Stokes problem like

-\Delta\mathbf{u}+\nabla p=\mathbf{f}

\nabla\cdot \mathbf{u}=0

Multiply \mathbf{u} on both sides of the first equation, by using Green’s formula, we can see it is just the saddle point problem, with A as -\Delta and B as \nabla:

\left[\begin{array}{cc} -\Delta & \nabla\\ \nabla & O \end{array}\right]\left[\begin{array}{c} \mathbf{u}\\ p\end{array}\right]=\left[\begin{array}{c} \mathbf{f}\\ 0\end{array}\right]

Here pressure acts as Lagrange multiplier, and the Lagrangian is

L(\mathbf{u}, p)=\frac{1}{2}|\Delta \mathbf{u}|_{\Omega}^2-(\mathbf{f}, \mathbf{u})+p(\nabla\cdot\mathbf{u}) \Longrightarrow

L(\mathbf{u})=\frac{1}{2}|\Delta \mathbf{u}|_{\Omega}^2-(\mathbf{f}, \mathbf{u})

where |\cdot|_{\Omega} is the seminorm on \Omega and \mathbf{u} is in H_0^1(\Omega).

Lagrangian representation for incompressible flow

Even Eulerian representation gives a field description of flow problems, it’s characterized and troubled by the convection terms. On the other hand, Lagrangian representation does not have this drawback, although the nonlinearity problem posed by convection term is replaced by nonlinearity of mapping between initial spacial cooridnates depending on particle index and spacial coordinates at some time later, as I am about to show.

Denote spacial coordinates, velocity, pressure in Eulerian representation as x_i, u_i(x,t), p(x,t) (symbol without subscript means the vector of proper dimension), one has the governing equations for inviscid incompressible flow:

\partial u_i/\partial x_i=0

\partial u_i/\partial t+u_j\partial u_i/\partial x_j=-\partial p/\partial x_i+\nu \partial^2 u_i/\partial x_jx_j

For some reference time t^0, and material coordinate (particle index) \xi, spacial position is described by a particle trajectory mapping from Lagrangian coordinate to Eulerian coordinate:

x_i=\phi_i(\xi,t)

Mapping \phi between material coordinate \xi and spacial coordinate x is called deformation.

When material coordinate coincides with spacial coordinate at reference time t^0, we have

x(\xi,t^0)=\phi(\xi,t^0)=\xi

in which x, \phi and \xi are all generally in \mathbb{R}^3. Then the Lagrangian representation of unknowns are

x=\phi(\xi,t)

u^t_i(\xi)=u_i(x,t)=u_i(\phi(\xi,t),t)

p^t(\xi)=p(\phi(\xi,t),t)

Here symbols with superscript denotes the variables in Lagrangian representation.  The material derivative and spacial/reference derivative, for some variable q are related by

\partial q/\partial x_i=(\partial q/\partial \xi_j)(\partial \xi_j/\partial x_i)= (\partial q/\partial \xi_j)(\partial \phi_i(\xi,t)/\partial \xi_j)^{-1}=F_{ij}^{-1}\partial q/\partial \xi_j

in which F_{ij} is the Jacobi matrix of \phi, i.e. the deformation gradient:

F=\nabla_{\xi}\phi=\frac{\partial \phi_i}{\partial \xi_j}

A second order tensor’s reference coordinate form and material coordinate form is connected by Piola transform:

T(\xi)=J\tau(x)F^{-t}\Longrightarrow \partial u_i/\partial \xi_j=J\partial u_i/\partial x_k F_{jk}^{-1}

where T is a second order tensor in reference configuration, and \tau is second order tensor corresponding to T  in deformed configuration. The superscript -t means “inverse” & “transpose”. Considering the divergence operation, we have:

\nabla_{\xi}\cdot T=\frac{\partial T_{ij}}{\partial \xi_j}

\nabla_x\cdot \tau=\frac{\partial \tau_{ij}}{\partial x_j}

The divergence theorem leads to

\int_{\Omega_0}\nabla_{\xi}\cdot T=\int_{\Gamma_0}Tn_0

\int_{\Omega}\nabla_x\cdot \tau=\int_{\Gamma}\tau n

Following theorem can be proved for the Piola transform:

\nabla_{\xi}\cdot T(\xi)=J\nabla_x\cdot \tau(x),\forall x=\phi(\xi,t)

T(\xi)n_0dS_0=\tau(x)ndS,\forall x=\phi(\xi,t)

and area derivatives dS_0 and dS are related by

dS=J|F^{-t}n^0|dS_0

Apply results above to the viscosity term in momentum equation, we have

\nu\frac{\partial^2 u_i}{\partial x_jx_j}=\nu J^{-1}\frac{\partial}{\partial \xi_j}(JF_{jk}^{-1}F_{kl}^{-1}\frac{\partial u_i}{\partial \xi_l})

Finally, the Lagrangian formulation for momentum equation is:

\frac{Du_i(\phi(\xi,t),t)}{Dt}=-F_{ij}^{-1}\frac{\partial p(\phi(\xi,t),t)}{\partial \xi_j}+\nu J^{-1}\frac{\partial}{\partial \xi_j}(JF_{jk}^{-1}F_{kl}^{-1}\frac{\partial u_i(\phi(\xi,t),t)}{\partial \xi_l})

J here bears the physical interpretation of volume change ratio, which is constantly one for incompressible flow. Depending on whether the reference time t^0 is taken as some time during the movement or always the initial time with initial configuration, the Lagrangian description can be divided into updated Lagrangian formulation or Total Lagrangian formulation.