Summary of Scheme, Boundary Conditions and Implementation

The splitting scheme involves the following four substeps:

This has been implemented in the code [30] for triangular and tetrahedral elements and in PRISM [40] for quadrilateral (nodal) elements. has now been recoded in an object oriented way so that it now handles all the element types we have described in chapter 2. The first step we call the non-linear step. We assume that the velocity fields are stored in terms of coefficients of the shape functions we outlined in 3. We then evaluate the velocity fields at quadrature points for each elements using a transform based on tensor product summation of the modes; this has been detailed in [30]. We now have a ``physical'' representation of the velocity at a set of points. For the non-linear step we take the spatial derivatives of these fields at the quadrature points and form the the non-linear convection terms by multiplying the velocity and its derivatives together at these points. We then take the current values of the non-linear terms and either one, two or three non-linear terms from previous times steps and add them together with a set of coefficients detailed in [39]. We now have the intermediate velocity as a field at the quadrature points.

For the second step we go through the list of boundary conditions for the elements and calculate the Neumann boundary condition for pressure using the compatability condition from equation 5.3.

For the third step we take the divergence of the intermediate velocity, using the same spatial element-wise derivative matrices as before. We then take the inner product over each element with that element's basis functions. We now have a set of coefficients for the right-hand-side of the pressure equation. Next we calculate the pressure field using the implicit solver we outlined in chapter 4, using the Schur complement of the global discrete Laplace operator.

For the last step we need to evaluate the pressure field at the quadrature points and then take its gradient also. From this we calculate the right-hand-side for the velocity solves. We take the inner product of the righthand side with the set of elemental basis functions. Again we use the implicit solvers to invert the Helmholtz equations for the velocity.