18.3 A basic operator-splitting algorithm

Next: 18.4 Flow in a Up: 18 SOLUTION OF 3D Previous: 18.2 Solution of a

18.3 A basic operator-splitting algorithm

The governing equations for incompressible laminar fluid flows are the Navier-Stokes momentum equation and the continuity equation. We will use the Fasttalk syntax to set up a basic operator-splitting algorithm to solve the 3D governing equations.

The Navier-Stokes momentum equation is:

$\begin{displaymath} \frac{\partial\mbox{\boldmath{$v$}}}{\partial t}+(\mbox{\bol... ... -\frac{1}{\mbox{Re}}\nabla\cdot[\nabla\mbox{\boldmath{$v$}}]=0\end{displaymath}$ (1)

where Re is the Reynolds number and $\mbox{\boldmath{$v$}}$ is the velocity vector.

For the unsteady Navier-Stokes equation, we solve (1) as an initial value problem, that is, at time step , we know the velocity value $\mbox{\boldmath{$v$}}^n$ , and we obtain the velocity at time step by solving for the velocity increment $\delta\mbox{\boldmath{$v$}}$ so that $\mbox{\boldmath{$v$}}^{n+1}=\mbox{\boldmath{$v$}}^n+\delta\mbox{\boldmath{$v$}}$ . From the Navier-Stokes equation, we see that $\delta\mbox{\boldmath{$v$}}$ satisfies:

$\begin{displaymath} \frac{\partial\delta \mbox{\boldmath{$v$}}}{\partial t}+(\mb... ... p +\frac{1}{\mbox{Re}}\nabla\cdot[\nabla\mbox{\boldmath{$v$}}]\end{displaymath}$ (2)

where due to the nonlinear convection term, the $\mbox{\boldmath{$v$}}$ on the RHS is the latest updated velocity value from $\mbox{\boldmath{$v$}}^n$ ; that is, $\mbox{\boldmath{$v$}}=\mbox{\boldmath{$v$}}^n+\delta \mbox{\boldmath{$v$}}$ .The pressure will be solved through the continuity equation in a way that is explained later.

The left-hand side of (1) can be split into two principal operators as follows:

$\begin{displaymath} \frac{\partial\delta \mbox{\boldmath{$v$}}}{\partial t}+(\mb... ... p +\frac{1}{\mbox{Re}}\nabla\cdot[\nabla\mbox{\boldmath{$v$}}]\end{displaymath}$ (3)

where $L_c(\delta \mbox{\boldmath{$v$}})$ and $L_d(\delta \mbox{\boldmath{$v$}})$ are the convective and diffusive operators, respectively.

This equation for velocity increment can be split as follows:

$\begin{displaymath} \frac{\partial\delta \mbox{\boldmath{$v$}}^{\ast}}{\partial ... ... p +\frac{1}{\mbox{Re}}\nabla\cdot[\nabla\mbox{\boldmath{$v$}}]\end{displaymath}$ (4)

$\begin{displaymath} \mbox{\boldmath{$v$}}^{\ast}=\mbox{\boldmath{$v$}}+\delta \mbox{\boldmath{$v$}}^{\ast}\end{displaymath}$

$\begin{displaymath} \frac{\partial\delta \mbox{\boldmath{$v$}}}{\partial t}-L_d(... ...ac{1}{\mbox{Re}}\nabla\cdot[\nabla\mbox{\boldmath{$v$}}^{\ast}]\end{displaymath}$

(5)

We can set up numerical schemes to integrate equations (4) and (5) respectively. We define the numerical scheme for equation (4) as

and the scheme for equation (5) as

. Equation (3) is advanced from timestep

by the two-stage convection-diffusion split:

$\begin{displaymath} \mbox{\boldmath{$v$}}^{n+1}=\mbox{\boldmath{$v$}}^n+N_dN_c\delta \mbox{\boldmath{$v$}}\end{displaymath}$

Alternatively we can set up a symmetric sequence of the numerical schemes to have

$\begin{displaymath} \mbox{\boldmath{$v$}}^{n+2}=\mbox{\boldmath{$v$}}^n+N_cN_dN_dN_c\delta \mbox{\boldmath{$v$}}\end{displaymath}$

(6)

In the Navier-Stokes equation, the convective operator

is nonlinear. Therefore, a local iteration loop is needed for equation (6).

The numerical scheme will be a fully implicit backward Crank-Nicolson method. The three components of velocity vector $\mbox{\boldmath{$v$}}$ can be solved separately by . In this way, less computer memory is required for inverting the linear matrices. For the -component , can be shown to be

$\begin{displaymath} \frac{\delta u_y^k}{DT}+(\mbox{\boldmath{$v$}}\cdot\nabla)\d... ...\partial y} +\frac{1}{\mbox{Re}}\nabla\cdot[\nabla u_y^{n+1,k}]\end{displaymath}$ (7)

where denotes the th iterative step for the nonlinear convective operator, and

$\begin{displaymath} u_y^{n+1,k+1}=u_y^{n+1,k}+\delta u_y^k\end{displaymath}$

Equation (7) is coded in Fasttalk as follows:

A   cony                  % name of the equation
% Nc for Uy Component
e  {1/DT} U1            \ % first term on LHS
  +{V100}_jD_j U1       \ % convective operator on LHS
  -D_j{1.0/Re}D_j U1    \ % diffusive term on LHS
 = {-1.0/DT} {V302}     \ % first term on RHS
  -{V100}_jD_j {V102}   \ % pressure gradient, Y dir'n
  +D_j{1.0/Re}D_j{V102} \ % diffusive term on RHS
%
b Uy=0      U1 = {0.0}        % b.c.: Uy = 0.0
b Uy=1      U1 = {1.0}        % b.c.: Uy = 1.0
b Uy=-1     U1 = {-1.0)       % b.c.: Uy =-1.0
b Uy=outlet n.grad U1 = {0.0} % b.c.:  Neumann Condition
%

where vector V100 contains the velocity vector $\mbox{\boldmath{$v$}}$ , the

velocity component resides in vector V102 and V201 contains the pressure value. Vector V300 contains the velocity vector increment $\delta \mbox{\boldmath{$v$}}^k$ , with the $\delta u_y^k$ value residing in vector V302.

Here the boundary conditions use symbolic tags such as Uy = 0. These have to be related to the numerical tags by previous D statements such as

   D     Uy=0    2
   D     Uy=0    3

A symbolic tag can stand for several numerical tags, and each numerical tag can have several symbolic tags associated with it.

The numerical scheme for solving the convection equation is constructed in Fasttalk as follows:

< solcony
!# convection operator for Uy     % print out message
cony 1 100 V100 200 V200 300 V300 % assemble Y conv eqn
M_rhs_L1 = L1 V101                % L1 norm of RHS
show M_rhs_L1                     % show RHS L1 norm
itsolve 1                         % call GMRES solver
show V101                         % show  solution value
V202 = V202 + V101                % U^{n+1}=U^{n+1}+dU^k
V402 = V101 + V402                % dU^{k+1}=dU^{k-1}+dU^k
popp                              % push vector stack back
>

The call to assemble the problem cony uses explicit argument passing since more than the default two levels are required. After the statement cony 1, which says the output is to appear starting at location 1, the pair 100 V100 says that global vector V100 is to be passed to the local level, and called V100 there. If we had wanted to call it V500 locally, we would have written 500 V100. The remaining pairs have similar effects.

Similarly, the numerical scheme can be set up for the component in a fully implicit backward Crank-Nicolson approach:

$\begin{displaymath} \frac{\delta u_y^k}{DT} -\frac{1}{\mbox{Re}}\nabla\cdot[\nab... ...\partial y} +\frac{1}{\mbox{Re}}\nabla\cdot[\nabla u_y^{n+1,k}]\end{displaymath}$ (8)

where $u_y^{n+1}$ is always updated at each step through $u_y^{n+1,k+1}=u_y^{n+1,k}+\delta u_y^k$ . This equation (8) is coded in Fasttalk as follows:

A  dify                        % name of the equation
% Nd for Uy Component
e   {1/DT} U1               \ % first term on LHS
   -D_j{1.0/Re}D_j U1       \ % diffusive term on LHS.
 = {-1.0/DT} {V302}         \ % first term on RHS
   -{V100}_j D_j {V102}     \ % convective operator on RHS
   -{0.0,1,0,0.0}_jD_j{V201}\ % pressure gradient, Y dir'n
   +D_j{1.0/Re} D_j {V102}  \ % diffusive term on RHS
%
b  Uy=0      U1 = {0.0}        % b.c.: Uy = 0.0
b  Uy=1      U1 = {1.0}        % b.c.: Uy = 1.0
b  Uy=-1     U1 = {-1.0}       % b.c.: Uy =-1.0
b  Uy=outlet n.grad U1 = {0.0} % b.c.: Neumann condition
%

The numerical scheme

for solving equation (5) is contained in macro soldify (see below). Note that this macro also makes use of explicit message passing in the same way as macro solcony.

< soldify
!# diffusion operator for Uy      % print out message
dify 1 100 V100 200 V200 300 V300 % Y diffusive eqn
M_rhs_L1 = L1 V101                % L1 norm of RHS
show M_rhs_L1                     % show RHS L1 norm
itsolve 1                         % call GMRES solver
show V101                         % show solution value
V202 = V202 + V101                % U^{n+1}=U^{n+1}+dU^k
V402 = V101 + V402                % dU^{k+1}=dU^{k-1}+dU^k
popp                              % push vector stack back
>

After we have built up

and

for the velocity component

, we can similarly code up the equations and solution procedures for the

and

components of the velocity vector $\mbox{\boldmath{$v$}}$ . Actually, the schemes for

and

are almost identical to that of

, the only difference being the direction of the pressure gradient.

After all the macros for the velocity components are constructed, they can be put together to produce a single macro for all velocity components, namely

< solcon
    solconx
    solcony
    solconz
>

where solconx and solconz are the macros for solving the convective equations of

and

respectively. For the diffusive equations, we construct macro command soldif given by

< soldif
    soldifx
    soldify
    soldifz
>

where soldifx and soldifz are the macros for solving the diffusive equations of

and

respectively.

As a summary, in accordance to equation (6), the operator-splitting algorithm for the momentum Navier-Stokes equation can be written in Fasttalk in the following macro opsplit.

< opsplit
    solcon
    soldif
    soldif
    solcon
>

In this algorithm, we have assumed that pressure

is known, and mass conservation is satisfied. The pressure

needs to be obtained through an equation that is derived from the mass conservation law. For incompressible fluid flows, this mass conservation law produces the continuity equation

$\begin{displaymath} \nabla\cdot\mbox{\boldmath{$v$}}=0\end{displaymath}$

An artificial compressibility method is used here to derive the equation for pressure. From the continuity equation, we can assume that the pressure

satisfies a pseudo-transient state

$\begin{displaymath} \frac{\partial p}{\partial\tau}=-\beta\nabla\cdot\mbox{\boldmath{$v$}}\end{displaymath}$

where parameter $\tau$ is the pseudo-time. Through a second-order Taylor expansion of the term $\partial p/\partial\tau$ and substitution of the Navier-Stokes momentum equation, we can arrive at a Poisson equation for pressure:

$\begin{displaymath} \delta p-\beta\frac{(D\tau)^2}{2}\nabla^2\delta p=-D\tau\beta\nabla\cdot\mbox{\boldmath{$v$}}\end{displaymath}$

where $D\tau$ is the pseudo-time step. This equation can be coded in Fasttalk as the following problem.

A  masscon                       % name of the cty eqn
e  {1.0}U1                      \% first term of eqn
   -D_j{0.5*Dtau*Dtau*beta}D_jU1\% Laplace operator
 = -{beta*Dtau}D_j{V100}_j       % div of velocity V100
b   inlet  U1 = {0.0}            % bc for flow inlet
b   outlet U1 = {0.0}            % bc for flow outlet
b   wall   n.grad U1 = {0.0}     % bc for wall

where velocity vector $\mbox{\boldmath{$v$}}$ is stored in the V100 vector stack location, as it was for the momentum equation. The pressure value

will be stored in the vector location V201.

After the pressure equation has been defined as shown above in Fasttalk, we can set up a Fasttalk macro to solve this equation. The macro command for solving the continuity equation is given the name solmass and it is defined in Fasttalk as follows:

< solmass
masscon 1 100 V100           % continuity eqn assembled
CON_rhs_L1 = L1 V101         % L1 norm of RHS of eqn
CON_rhs_L1 = CON_rhsL1/Nodes % average residual per node
show CON_rhs_L1              % show L1 norm of RHS 
itsolve 3                    % call iterative CG solver
show V101                    % show  max/min values of
                             % pressure increment
V301 = V101 + V301           % add increment to pressure
popp                         % push vector stack down
>

This macro solmass can be inserted into the macro opsplit for the Navier-Stokes momentum equation, so that the mass conservation law is satisfied when solving for the velocity. Macro opsplit for the operator-splitting algorithm for incompressible laminar fluid flows takes its final form as follows.

< opsplit
    solcon   % solving  the  convective operator
    zero     % after first iteration, bc's are set zero
    solmass  % solving for pressure
    soldif   % solving for diffusive operator
    soldif   % solving for diffusive operator
    solcon   % solving the convective operator
    solmass  % solving for pressure
>

You can invoke the operator-splitting code to predict a 3D incompressible flow by simply typing opsplit. This will produce a solution at the first timestep. Since the convective operator is nonlinear, an iterative process is needed to reach a converged solution for the timestep. A very simple loop can be set up to achieve this in Fasttalk:

< one_step                    % name of this macro
    while(CON_rhs_L1 > 0.0001)% cty cnvgnce criterion
      opsplit                 % call o-s algorithm
    endwhile                  % end of while loop
    show Time                 % print on screen the time
    store Time                % store Time in rhs.out
    store V100                % store vel vector V100 in
                              % file rhs.out
>

Timestepping from time

to final time

is accomplished by macro run.

< run                % name of this macro
  Time = T0          % starting time
  while (Time < Tn)
    one_step         % calls operator-splitting algorithm
                     % at this time step
    Time = Time + DT % a new time step is started
    V300=0.0         % velocity increment is reset to
                     % zero at new time step
  endwhile
>

We have built up an operator-splitting algorithm for incompressible fluid flows in 3D. This algorithm is in finite element formulation, and it can be used on 3D mesh types that are supported by Fastflo. For example, they can be 8-node or 20-node brick elements, or 4-node or 10-node tetrahedron elements.