15.6 Stefan and moving boundary problems

Next: 16 MESH GENERATION (ADVANCED) Up: 15 SPECIAL FEATURES Previous: 15.5 Irrotational flow around

15.6 Stefan and moving boundary problems

There are various ways in Fastflo of handling moving boundary problems. One is to simply move the node points of the mesh to follow the fluid, producing a Lagrangian frame of reference. Another is the Arbitrary Lagrangian-Eulerian (ALE) approach, in which the boundary fits the new free surface, but a mesh smoothness condition regulates the movement of interior nodes and the velocity difference between fixed and transformed meshes appears as a convective velocity in the equations.

In this section, we develop the ALE approach in Fastflo. Analysis is done on a fixed reference mesh; a mapping function computed on this mesh describes the relation to the real space with moving free surface. Large deformations can be accommodated without undue mesh distortion. The method has widespread application, and is illustrated for a 2D Stefan problem in which we compute the freezing of a fluid region.

Let the true space coordinate vector be x, and let

$\begin{displaymath} {\bf X} = {\bf X}({\bf x},t)\end{displaymath}$

be a smooth mapping into a domain with fixed boundaries, which is meshed.

The X variables can be regarded as a set of time-varying curvilinear coordinates, and various differential operators encountered in partial differential equations are transformed as in the following examples

in which $W_{ij} = (\partial X_j / \partial x_i )$ and $J = ({\rm det}~W)^{-1}$

A time derivative can be written

$\begin{displaymath} \left( \frac{\partial A}{\partial t}\right)_{x} = \left( \... ...rtial t}\right)_{X} \frac{\partial }{\partial X_k}(JW_{jk}A)\end{displaymath}$

and the generic convective term can be written as

$\begin{displaymath} \left( \frac{\partial A}{\partial t}\right) _{x} + \frac{\... ...x_j}{\partial t} \right) _{X} \right) JW_{jk}A \right] \right]\end{displaymath}$

With these transformations, the entire set of equations can be solved in the fixed-mesh X-domain. The mapping from x into the fixed geometry X can be chosen to preserve internal smoothness, together with some necessary boundary conditions.

By way of illustration, we wish to solve a Stefan problem in which the dominant equation is the time-varying heat equation

$\begin{displaymath} \rho c_p \frac {\partial T}{\partial t} = \nabla_x.(k\nabla_x T)\end{displaymath}$

on a varying domain $\Omega$ in the x-plane. We make the mapping to the fixed X domain, in which case the equation to be solved is

$\begin{displaymath} \rho c_p \left[ J \left( \frac{\partial T}{\partial t}\right... ...j}\left(JW_{mj}W_{mk} \frac{\partial A}{\partial X_k}\right) \end{displaymath}$

A key point of the ALE method is that the motion of the internal mesh nodes can be assigned arbitrarily. For convenience, we choose to calculate the ``velocity''

$\begin{displaymath} {\bf v} = \left( \frac{\partial {\bf x}}{\partial t}\right) _{X}\end{displaymath}$

by solving the elasticity-like equation

$\begin{displaymath} \frac{\partial}{\partial X_j} \left(\frac{\partial v_i}{\partial X_j}+ \frac{\partial v_j}{\partial X_i}\right)=0\end{displaymath}$

This equation is simply added to the overall set. The required boundary condition is that the resolute of v normal to the surface in the x-plane be equal to the normal velocity required on physical grounds.

Finally, the update of the mapping

$\begin{displaymath} {\bf x} = {\bf x}({\bf X},t)\end{displaymath}$

is given by the approximate time integral

$\begin{displaymath} {\bf x}({\bf X}, t+\delta t) = {\bf x}({\bf X}, t) + {\bf v} ({\bf X}, t)\delta t\end{displaymath}$

[

Some caution is required here! This update is not fully implicit, whereas the mainstream finite element methods used in Fastflo are implicit. Stability restrictions on timesteps can be expected to occur in the ALE method.]

We now consider a two-dimensional problem.

The domain chosen in the X-plane is the unit square, with the free surface at . The condition at the free boundary is given by the Stefan condition:

$\begin{displaymath} k ~{\bf n}.\nabla T = L\rho ~{\bf n}.\frac{d{\bf x}}{dt}\end{displaymath}$

Normally one would think of n as a unit normal, but because it appears on both sides, it can be of any magnitude.

In general, the free surface may be defined by a relation

$\begin{displaymath} F({\bf X}) = 0\end{displaymath}$

so a normal is given by

$\begin{displaymath} \frac{\partial F}{\partial x_i} = W_{ij}\frac{\partial F}{\partial X_j}\end{displaymath}$

and the Stefan condition becomes

$\begin{displaymath} k \frac{\partial F}{\partial X_j}W_{mj}W_{mk} \frac{\partia... ...L\rho W_{ij}\frac{\partial F}{\partial X_j}\frac{d{\bf x}}{dt}\end{displaymath}$

In our case,

, and

$\begin{displaymath} W_{ij}\frac{\partial F}{\partial X_j}= W_{i1}\end{displaymath}$

Suppose the original -domain was the thin film of ice and , with an initial temperature distribution given by

$\begin{displaymath} T({\bf x}, t) = (T_0+x_2^2)(50x_1 -1)\end{displaymath}$

The temperature at

is held at the original quadratic variation

, and the surfaces at

and

are insulating. To complete the problem, imagine that there is an infinite pool of water at freezing temperature occupying the region adjacent to the free boundary. In this case, no heat conduction takes place in the fluid, and the region does not figure in the computations.

The mesh file stefan2D.msh is produced using Unit and consists of a unit square in the X-plane. It has 64 equal 8-noded squares as shown in the figures at the end of this section. The boundary segments are tagged as follows:

$\begin{displaymath} 0 < X_1 < 1, ~~X_2 = 0, ~~~ {\rm tag~3 ~~(insulating)}~~~~~~~~~~~~~~ \end{displaymath}$

$\begin{displaymath} X_1 = 1, ~~0 < X_2 < 1, ~~~ {\rm tag~2 ~~(free~boundary)}~~~~~~~~~~~ \end{displaymath}$

$\begin{displaymath} 0 < X_1 < 1, ~~X_2 = 1, ~~~ {\rm tag~3 ~~(insulating)}~~~~~~~~~~~~~~ \end{displaymath}$

$\begin{displaymath} X_1 = 0, ~~0 < X_2 < 1, ~~~ {\rm tag~1 ~~(prescribed~temperature)} \end{displaymath}$

The problem file is stefan2D.prb . At the start of this file, certain parameters are defined. These parameters prescribe respectively the stopping time t1 for the computation (in seconds), the timestep dt (in seconds), the latent heat of fusion

, the thermal diffusivity

, the density $\rho$ , the specific heat

, and the temperature

. The values are for ice and SI units are used.

The next part of the file specifies two problems, heat and spread. To understand this part of the code, it is instructive to review the equations to be solved and the transformations (1-3) at the start of this section.

Problem heat is the transformed heat equation. V200 is the current value of x, and the grad operation puts into the locations V400, V500 the values of the matrix $Q_{ij}=( \partial x_j / \partial X_i)$ . Note that the size of the object returned, namely a 2 $\times$ 2 matrix, is determined by the size of V200 and the operation performed on it, so the results are placed in V400 and the next consecutive locations. Then the inverse of the matrix $T_{ij} = (Q^{-1})_{ij} = (\partial X_j /\partial x_i)$ is placed in V600, V700. V800 is just formed as the transpose of , and so the product of with its transpose, $T_{mi}T_{mj}$ is placed in V1000, V1100. V1201 contains the determinant , and V1300 is the product $J (\partial {\bf x} / \partial t)$ .

With all that decided, the transformed heat equation is written down with a fully implicit time discretisation. The timestepping procedure is the backward Euler method (see Section 7.3.1). The boundary conditions at and simply say that the new temperature is the same as the old one. At the other boundary segments, the absence of an explicit condition ensures that there is no normal heat flux.

The problem spread describes the equation for updating the transformation, by timestepping the function ${\bf x}({\bf X},t)$ . In fact it computes $( d{\bf x} /dt)$ . The transformation matrices are computed as before, although they are now required only for the boundary condition at the free surface. The components of $( d{\bf x} /dt)$ normal to the fixed surfaces are set to zero, and the Stefan condition is applied, as discussed above, to the free surface.

The remainder of the specification consists of a macro run. This successively updates the temperature (using the old mesh velocity) and then the mesh (using the new temperature). The temperature and mesh are not strongly coupled. Note that, when heat is assembled, various registers in the vector stack need to be passed from the global to the local position. This explains the additional addresses mentioned after the name heat. After solving both heat and spread, the stack is popped to bring the registers back to their original designation. The command mapm 200 in each timestep swaps the contents of V200 with the mesh locations and enable a contour plot of the temperature at the physical locations to be drawn The second mapm 200 then restores V200 to the transformed mesh.

The sequence is run over 4,000,000 seconds (a little more than ``40 days and 40 nights'') with 200 equal timesteps, and the extent of the frozen region and the temperature contours at the end are as shown below.

The ALE method has the virtue of preserving the front without diffusion, and of allowing more physics to be included at the interface if necessary. For example, surface tension would not be difficult to include. A limitation of the ALE method is that it cannot easily handle changes of topology, such is the meeting of two fronts. The method described here can be straightforwardly generalised to three dimensions.

Finally, we return to the comment made earlier about timestep restrictions. If the timestep is made too large, the explicit scheme for updating the mesh breaks down.