PDE « Mesh Mess

What’s the difference?

Posted in FEM by Yi Zhang on 01/04/2010

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.

Tagged with: FEM, Finite difference, numerical, PDE, Sobolev space

leave a comment

Saddle point problem in Navier-Stokes equation

Posted in Flows by Yi Zhang on 12/21/2009

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)$ .

Tagged with: CFD, numerical, PDE

leave a comment

Lagrangian representation for incompressible flow

Posted in Flows by Yi Zhang on 07/24/2009

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.

Tagged with: CFD, PDE

leave a comment

Sobolev space: I

Posted in FEM by Yi Zhang on 07/04/2009

Sobolev space $W_p^l$ is based on generalize derivatives and the corresponding norm induced. As a special case, space $H^l$ has the meaning of finite “intensity” in general physical sense. For example, if $u(x)$ is the flow velocity in spatial domain $\Omega$ , then $u\in H^1(\Omega)$ means the function and its first order derivatives are all in $L^2(\Omega)$ , i.e.

$\int_\Omega{\nabla u \cdot \nabla u+u^2}d\Omega$

Physically it means the total intensity of the sources within $\Omega$ is finite, so is the total kinetic energy within the domain. Obviously this requirement is fundamental for any engineering analysis.

Meyers and Serrin proved that $W_p^l\bigcap C^\infty$ is dense in $W_p^l$ under the Sobolev space norm. This means that the function in Sobolev space can be approximated using smooth functions, and in many occasions the properties of Sobolev space is obtained by first demonstrate it in $C^\infty$ .

The fundamental question of approximating the variational form of a PDE is twofold:

How much the smoothness/differentiability one can have from the weak solution?

How close the boundary condition can be approximated?

The first question is referred to as regularity problem. The properties of embedding, and traces of the Sobolev spaces are fundamental to answer questions above. It must be emphasized that the regularity of a weak solution $Lu=f$ depends on the data $f$ , the geometry of domain, and the boundary condition. All three factors must be examined in practical finite element convergence analysis.

Tagged with: CFD, FEM, PDE, Sobolev space

leave a comment

Time integration: incompressible flow

Posted in Flows by Yi Zhang on 06/28/2009

I try to consolidate last post using a scheme for incompressible flow under Lagrangian description. Those two condition make somehow make the plan much easier against general case. Since for the moment I am reading this paper, I will take some of its ideas, even though this whole modern time splitting scheme by projection methods owns its origin to Chorin and Temam in 1960′s.

Consider the governing equations:

$\frac{D\rho}{Dt}=-\rho\frac{\partial u_i}{\partial x_i}$ (MEL-in mixed-Eulerian-Lagrangian form) or $\frac{\partial u_i}{\partial x_i}=0$ (Eulerian form)

$\frac{Du_i}{Dt}=-\frac{1}{\rho}\frac{\partial p}{\partial x_i}+\mu\frac{1}{\rho}\frac{\partial^2 u_i}{\partial x_j\partial x_j}+f_i$

For a specific material particle $\xi$ , both the displacement and velocity are functions of time:

$u=u(\xi,t)\quad x=x(\xi,t)$

Now, let’s suppose the variables at time $t_i$ are known as $\phi^n=\phi(\xi,t_n)$ . For velocity we have the general time integration:

$\frac{u_i^{n+1}-u_i^n}{\Delta t}=(1-\theta)(-\frac{1}{\rho^n}\frac{\partial p^n}{\partial x_i^n}+\mu\frac{1}{\rho^n}\frac{\partial^2 u_i^n}{\partial x_j^n\partial x_j^n}+f_i^n)$

$+\theta(-\frac{1}{\rho^{n+1}}\frac{\partial p^{n+1}}{\partial x_i^{n+1}}+\mu\frac{1}{\rho^{n+1}}\frac{\partial^2 u_i^{n+1}}{\partial x_j^{n+1}\partial x_j^{n+1}}+f_i^{n+1})$

In equation above, partial derivative by $x_i^{n+1}$ means take the derivative within the spatial configuration of time $t_{n+1}$ . This spatial domain is obtained using velocity updating:

$x(t_{n+1})=x(t_n)+u(t_n)\Delta t$

Parameter $\theta$ varies between 0 and 1. When $\theta=0$ , RHS is based on time $t_{n+1}$ , i.e. implicit scheme, while when $\theta=0$ we have explicit scheme. Generally, implicit schemes are prefered for its stability merit, which I adopt in this post. So we have

$\frac{u_i^{n+1}-u_i^n}{\Delta t}=-\frac{1}{\rho^{n+1}}\frac{\partial p^{n+1}}{\partial x_i^{n+1}}+\mu\frac{1}{\rho^{n+1}}\frac{\partial^2 u_i^{n+1}}{\partial x_j^{n+1}\partial x_j^{n+1}}+f_i^{n+1}$

There are two choices to implement the discrete governing equations. One is two resolve $u^{n+1}_1,u^{n+1}_2,u^{n+1}_3,p^{n+1}$ all together, by coupling equation above with mass conservation equation of Eulerian form $\nabla\cdot \vec{u}=0$ . The other one is try to update velocity separately, for this purpose we have to transform the mass equation, because it is this equation that couples velocities in different directions. On the other hand, we also want to decouple the pressure and velocity. Before introduce the scheme, let’s first take a look at the physical interpretation of the coordnate system.

When solving the dynamical problem numerically, we hope at one time step, we take care of time derivatives and space derivatives separately. Suppose I am examining a fixed time step $t_n\rightarrow t_{n+1}$ , the material particle indexed by $\xi$ can actually be indexed by the spatial position (in Eulerian coordinates) $\vec{x}$ . In other words, when we talk about some specific material particle, it can be labeled with location $\vec{x}$ at time $t_n$ . So $x^{n+1}(\xi)=x^{n+1}(x)$ simply means the $t_{n+1}$ location of some particle which at $t_n$ occupies $x$ . By this notation, the momentum equation, assuming constant density and body force, is

$\frac{u_i^{n+1}(x^{n+1})-u_i^n(x)}{\Delta t}=-\frac{1}{\rho}\frac{\partial p^{n+1}}{\partial x_i^{n+1}}+\mu\frac{1}{\rho}\frac{\partial^2 u_i^{n+1}}{\partial x_j^{n+1}\partial x_j^{n+1}}+f_i$

and the mass equation is

$\frac{\rho^{n+1}(x^{n+1})-\rho^{n}(x)}{\Delta t}=-\rho\frac{\partial u_i}{x_i}$

In order to decouple velocity and pressure in momentum equation, we split it:

$\frac{u_i^{n+1}(x^{n+1})-u^{\ast}(x)}{\Delta t}+\frac{u^{\ast}(x)-u_i^n(x)}{\Delta t}=[-\frac{1}{\rho}\frac{\partial (p^{n+1}-p^n)}{\partial x_i}]+[-\frac{1}{\rho}\frac{\partial p^n}{\partial x_i}+\mu\frac{1}{\rho}\frac{\partial^2 u_i^{\ast}}{\partial x_j\partial x_j}+f_i]$

The two terms on both side form pairs accordingly. The second terms gives an implicit scheme for intermediat variable $u^{\ast}$ . Similar split in mass equation gives

$\frac{\rho}{\Delta t}\frac{\partial u^{\ast}_i}{\partial x_i}=\frac{\partial^2 p^{n+1}}{\partial x_i^2}$

which gives $p^{n+1}$ , then the first terms in split momentum equation is used to obtain $u^{n+1}_i$ :

$u^{n+1}_i=u^{\ast}_i-\frac{\Delta t}{\rho}\frac{\partial (p^{n+1}-p^n)}{\partial x_i}$

The three equations above, are to be used to find $u^{\ast},p^{n+1},u^{n+1}$ sequentially. This means to solve two spatial domain PDEs for the first two equations, and take one derivative for the last.

Tagged with: CFD, FEM, Math, numerical, PDE, time integration

2 comments

Older Posts »

M	T	W	T	F	S	S
« Mar
	1	2	3	4	5	6
7	8	9	10	11	12	13
14	15	16	17	18	19	20
21	22	23	24	25	26	27
28	29	30	31

Mesh Mess

What’s the difference?

Saddle point problem in Navier-Stokes equation

Lagrangian representation for incompressible flow

Sobolev space: I

Time integration: incompressible flow

Recent Posts

Tags

Meta

Blogroll

Follow “Mesh Mess”