Time integration: incompressible flow

by Yi Zhang

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.