15.4 Arbitrary mesh mappings

Next: 15.5 Irrotational flow around Up: 15 SPECIAL FEATURES Previous: 15.3 Numerical errors in

15.4 Arbitrary mesh mappings

The mapmesh function is used particularly for free surface flows. A useful special form of this application is the ability to map a mesh smoothly onto a reasonably similar shape, specified by boundary constraint functions.

As an example of the use of mapmesh, we set up an idealised problem. Suppose we want to map the initial mesh onto a region specified by

$\begin{displaymath} F_i(x) = 0\mbox{ with tags } i = 0,1,\ldots\end{displaymath}$

We suppose the region is occupied by a fluid of kinematic viscosity $\nu$ which, for simplicity, has no resistance to compression. The Reynolds number is assumed to be very low, so momentum effects can be ignored. Consequently, the pressure can also be assumed to be zero, and the equation for fluid motion is a balance of acceleration and shear stress. The equation, discretised in time by a fully implicit method, is

$\begin{displaymath} \frac{v_i-v_i^{\ast}}{\delta t}-\frac{\partial}{\partial x_j... ...artial x_j}+ \frac{\partial v_j}{\partial x_i}\right)\right]=0\end{displaymath}$

where $\mbox{\boldmath{$v$}}$ is the velocity at the end of the timestep, and $\mbox{\boldmath{$v$}}^{\ast}$ is the starting velocity. The full stress expression is needed here because there is no incompressibility. [

For discussion of the second term inside the square brackets above, see Section 14.1.]

The boundary condition on the boundary tagged will be a normal stress proportional to the constraint $F_i(\mbox{\boldmath{$x$}})$ . It is applied in the sense that it is positive if the mesh needs to be squeezed inwards to conform to the boundary.

The nodes of the mesh will be moved in accordance with the resulting velocity, except for the surface nodes which will be moved in the direction normal to the surface. So, ignoring changes in the normal direction over the timestep, the boundary stress at the end of the step on the boundary segment with tag will be

$\begin{displaymath} \left[F_i+(\nabla F_i\cdot\mbox{\boldmath{$n$}})(\mbox{\bold... ...elta v$}}\cdot\mbox{\boldmath{$n$}}\right]\mbox{\boldmath{$n$}}\end{displaymath}$

We now assume that the velocities are small and the timestep is very long, so the time derivative can be dropped. The stresses will be applied as a natural boundary condition. The unknown in the equations is the displacement $\mbox{\boldmath{$v$}}\delta t$ . The viscosity is made small to ensure the boundary displacement is large in the timestep. With a larger value, the convergence will be slower, but more stable.

To illustrate the use of mapmesh, we convert a rectangle to conform with the curved shape with boundaries

1.
2.
3.
4.

The initial rectangle is $1.6\times 1$ with square elements side 0.1 as shown above. The mesh specification for this initial rectangle is contained in file distort.msh which was created using the tool Unit as described in Section 16.6. distort.msh is quite a large file since it contains a full listing of nodes and tags, as well as information on node connectivity.

The Fasttalk code that distorts the rectangle into the curved shape is distort.prb . This code starts by defining the parameter nu and assigning the values 1-4 to the parameter all. This use of the D statement saves time and reduces typing errors.

The main problem is called stress. This expresses the calculations at the heart of the mesh distortion algorithm. The various terms in stress are explained below.

Macro run controls the distortion. The first mapmesh command transfers the coordinates of the original mesh onto V101 and V102. In the original domain as generated by Unit, the rectangle had dimensions $1.6\times 1.0$ .The commands

   extent = max V101
   V101 = V101/extent

scale the

-dimension of the mesh to 1.0, whilst the second mapmesh command then creates a unit rectangle by transferring V101, V102 onto the mesh storage locations. Finally, macro run makes three calls to macro step whose effect is explained below.

We now look at the various statements inside problem stress.

   A   stress
   e   - D_j{nu}D_j U1_i - D_j{nu}D_iU1_j
   b   1   V301 = -X2 + X1*(1-X1)
   b   2   V301 = -X1 + 1 + X2*(1-X2)
   b   3   V301 =  X2 - 1  + X1*(1-X1)
   b   4   V301 =  X1 + X2*(1-X2)

The A statement declares that the problem is called stress. The e statement calculates the stress terms in the equation that defines the timestepping procedure. The four b statements define the four constraints

on the four boundary segments $i=1,\ldots,4$ .

   b   all  V200 = normal
   b   all  V400 = [grad] {V301}
   b   all  V302 = V400@V200

These statements make calculations required for the boundary stress term

$\begin{displaymath} \left[F_i+(\nabla F_i\cdot\mbox{\boldmath{$n$}})(\mbox{\bold... ...lta v$}}\cdot\mbox{\boldmath{$n$}})\right]\mbox{\boldmath{$n$}}\end{displaymath}$

on all boundaries. Firstly the two components of the normal to the boundary are stored in vectors V201, V202. The next statement calculates $\nabla F_i$ and stores the result on V401, V402. Then the scalar product $\nabla F_i\cdot\mbox{\boldmath{$n$}}$ is calculated and stored in V302. The next two statements

   b   all  V500 = V302, V302
   b   all  V500 = V500 * V200

calculate and store $(\nabla F_i\cdot\mbox{\boldmath{$n$}})\mbox{\boldmath{$n$}}$ on V501, V502.

   b   all  V302=V301
   b   all  {V200*V300}_i+{V500 t V200}_ijU1_j

These statements copy the constraints $F_i(\mbox{\boldmath{$x$}})$ from V301 to V302, and then complete the calculation of the boundary stress term

$\begin{displaymath} \left[F_i+(\nabla F_i\cdot\mbox{\boldmath{$n$}})(\mbox{\bold... ...lta v$}}\cdot\mbox{\boldmath{$n$}})\right]\mbox{\boldmath{$n$}}\end{displaymath}$

   < step
       stress
       show L1 V100
       solve
       mapm 200
       V200=V200+V100
       mapm 200
       clear
       prim
       popp
   >

Macro step first assembles the finite element representation of the distortion problem all and then prints out the L1 norm of V100 which gives a measure of the boundary stress (which is in turn proportional to the deviations of the boundary constraints). The sparse matrix system is then solved by Fastflo's direct solver, and the current mesh coordinates are transferred to V201, V202 by the command mapm 200. The new mesh location is calculated and transferred back to the required special storage locations for meshes. Finally, the graphics screen is cleared, the new mesh displayed, and the vector stack popped ready for the next application of macro step.

The output of the computations is

   L1 V100 =  0.662031
   L1 V100 =  0.000357158
   L1 V100 =  3.80661e-06

and the distorted mesh is shown below.

Different physical models can be adopted. A degree of resistance to compression (a gas law, for example) might be desirable. The boundary nodes could be allowed to move sideways to minimise stress, although this requires differentiating the boundary constraint expressions.