Spatial discretization

The unknown functions of velocity and pressure are approximated using equal order interpolation for all variables in the final configuration

$$u_i = \mbox{\boldmath$N$} ^T(X,t) \mbox{\boldmath$U$} _i$$ (2.53)

$$p = \mbox{\boldmath$N$} ^T(X,t) \mbox{\boldmath$P$} \;.$$ (2.54)

By applying the Galerkin weighted residual method on the splitted governing equations, following system of linear algebraic equations is obtained

$$\mbox{\boldmath$M$} \mbox{\boldmath$U$} ^* = \mbox{\boldmath$M$} \mbox{\boldmath$U$} ^n + \Delta t \mbox{\boldmath$F$} - \frac{\Delta t\mu}{\rho} \mbox{\boldmath$K$} \mbox{\boldmath$U$} ^{n+\theta}\;,$$ (2.55)

$$\mbox{\boldmath$L$} \mbox{\boldmath$P$} ^{n+1} = \frac{\rho}{\Delta t}\left(\mbox{\boldmath$G$}^T\mbox{\boldmath$U$}^*-\hat{\mbox{\boldmath$U$}}\right) \;,$$ (2.56)

$$\mbox{\boldmath$M$} \mbox{\boldmath$U$} ^{n+1} = \mbox{\boldmath$M$} \mbox{\boldmath$U$} ^* - \frac{\rho}{\Delta t} \mbox{\boldmath$G$} \mbox{\boldmath$P$} ^{n+1} \;.$$ (2.57)

The matrix $M$ denotes the mass matrix in a lumped form, whereas the vector $F$ stands for the load vector. The matrix $G$ represents the gradient operator, which is the transposition of the divergence operator denoted simply as $G$$^T$. Matrices $K$ and $L$ are build in a similar way however noted differently. Both mean the Laplacian operator. Due to its common use in computational mechanics, the classical notation of the stiffness matrix $K$ is used. Prescribed velocity components are enclosed in vector $\hat{\mbox{\boldmath $U$}}$.

In each computational time step, an iteration is performed until the equilibrium is reached. Depending on the value of $\theta$ used, the equation system for the components of the auxiliary velocity $U^*_i$ (2.55) can be solved either explicitly $\theta=0$ or implicitly $\theta \neq 0$. Then, the calculated values of the auxiliary velocity are used as input for the pressure computation (2.56). The last system of equations (2.57) determines the velocity values at the end of the time step, taking auxiliary velocities and pressure or pressure increments into account.

Let us summarize the iterative step. The position of the particles at the end of the previous time step is known, as well as the the value of the velocity $u^n$ and pressure $p^n$. The set of governing equations is build up for the unknowns at the end of the solution step $\theta^{n+1}$, however based on the geometry of the previous step. The changes in the position are neglected. Once the convergence is reached, the final position is computed from the old one modified by the displacement due to obtained velocity. After that, solution can proceed to the next time step.