Temporal Discretization

In order to describe in simple terms the solution algorithm, it is convenient to write equation (5.8) in matrix form. Thus, we introduce the mass matrix $\mm \equiv (\phi({\bf x}), \phi({\bf x}))$, the derivative matrix $\dm \equiv (\phi, \phi_{\bf x})$, and the stiffness (Laplacian) matrix $\lm \equiv (\phi_{\bf x}, \phi_{\bf x})$, where the basis functions $\phi({\bf x})$ are defined in section 3. We also denote the nonlinear contributions by $\nl(U,W)$ to denote dependence on the two velocities, i.e. the flow velocity field U and the mesh velocity W. Using this notation, equation (5.8) becomes  

$$\frac{d}{d t} (\mm U_i) + \nl_i(U,W) = \dm_i^T P - \lm_{ij}U_j + F_i.$$ (3)

and the incompressibility condition (5.6) is  

$$\dm_i U_i = 0.$$ (4)

The mesh velocity is defined based on the mesh coordinates from  

$$\frac{d X_i}{d t} = W_i.$$ (5)

In the above equations we have denoted the fields with capital letters to denote the discrete version of the continuous field, e.g. $u_i \mapsto U_i$, etc.

Before we proceed with the boundary conditions, we discretize equations (5.9, 5.10) and (5.11) in time using a splitting scheme and a third-order stiffly stable integration scheme (see reference [55]).

We solve for time step (n+1); first treating explicitly the nonlinear and the mesh velocity terms:

Next we treat the elliptic terms implicitly

where we compute the intermediate field from

$$\hat{\hat{U}}_i = \mm \hat{U}_i + \dm_i^T P^{n+1}.$$

The constants $\alpha_q, \beta_q$ are integration weights and are defined in [55]. The mesh velocity is in general arbitrary, and it can be specified explicitly or be obtained from a Laplace equation as done here following [54]. More recent work in [53] suggests a modified approach where a variable coefficient is used in the Laplacian in order to provide enhanced smoothing and thus prevent sudden distortions in the mesh.