Mesh deformation using the biharmonic operator - Semantic Scholar
INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING Int. J. Numer. Meth. Engng 2001; 00:1–30 Prepared using nmeauth.cls [Version: 2000/01/19 v2.0]
Mesh deformation using the biharmonic operator Brian T. Helenbrook
∗
Mechanical and Aeronautical Engineering Department Clarkson University Potsdam, NY 13699-5725
SUMMARY The use of the biharmonic operator for deforming a mesh in an arbitrary-Lagrangian-Eulerian simulation is investigated. The biharmonic operator has the advantage that two conditions can be specified on each boundary of the mesh. This allows both the position and the normal mesh spacing along a boundary to be controlled, which is important for two-fluid interfaces and periodic boundaries. At these boundaries, we can simultaneously fix the position of the boundary and ensure that the normal mesh spacing is continuous across the boundary. In addition, results for deforming surfaces show that greater surface deformation can be tolerated when using biharmonic equations compared to approaches using second-order partial differential equations. A final advantage is that with the biharmonic operator, the integrity of a grid in a moving boundary layer can be preserved as the boundary moves. The main disadvantage of the approach is it’s increased computational expense. c 2001 John Wiley & Sons, Ltd. Copyright key words: adaptation
arbitrary-Lagrangian-Eulerian (ALE), mesh deformation, mesh movement, r-
Introduction In a typical arbitrary-Lagrangian-Eulerian (ALE) flow simulation, a boundary of a mesh tracks some moving/deforming feature of interest. This is usually a Lagrangian feature such as a cylinder surface or a free-surface wave. As the boundary of the mesh moves, the interior region of the mesh is deformed “arbitrarily” to maintain a quality mesh for a Eulerian flow solution. Because the mesh deformation is somewhat arbitrary, there have been several different methods proposed for doing this. These methods typically compromise between the amount of boundary deformation that can be tolerated and the computational expense of the procedure.
∗ Correspondence
to: Mechanical and Aeronautical Engineering Department Clarkson University Potsdam, NY 13699-5725 e-mail: [email protected] phone: 315-268-2204 fax: 315-268-6438
c 2001 John Wiley & Sons, Ltd. Copyright
2
B. T. HELENBROOK
Most methods use some form of a second-order partial differential equation (PDE). Typical examples are the Laplacian method [1, 2], the spring method [3, 4, 5, 6], and the elastic deformation method [7, 8, 9]. A limitation of a second-order method is that at boundaries either the mesh position or the normal mesh spacing can be specified, but it is impossible to specify both. A typical example where we might want to control both is a moving cylinder. At the surface of the cylinder, we would like to control the position of the boundary as well as the boundary-layer resolution i.e. the normal mesh spacing. With a second-order PDE, this is impossible. Another example is the simulation of two-fluid flows. In this case, the edges of two meshes meet to define the position of an interface between two different fluids. This mesh boundary is deformed to track the position of the interface and thus is a Dirichlet boundary condition for the mesh deformation PDE. However, we would also like to guarantee that the two meshes maintain a smooth gradation in mesh spacing across the interface. This requires that we specify both the position and the normal spacing at the interface boundary, which is impossible with a second-order PDE. To overcome this limitation, the use of a fourth-order PDE is proposed for determining the interior mesh deformation. Fourth-order generalizations of any of the three second-order deformation methods can be created, however this work focuses on a generalization of the Laplacian method because it is the most straightforward. The first section of this paper describes the use of the biharmonic operator as a fourth-order generalization of the Laplacian operator for mesh deformation. In this section, the appropriate choice of boundary conditions for modeling different types of boundaries in an ALE flow simulation is also discussed. The following section presents a multigrid solution method and examines the relative cost of the biharmonic versus the Laplacian methods. Results are obtained using both methods for several cases including a deforming free-surface, a deforming two-fluid interface, a moving cylinder in a channel, and a moving cylinder with a refined boundary layer mesh.
Governing Equations A typical method used to deform a mesh is the following system of Laplacian equations ∇2 X = 0 ∇2 Y = 0
(1)
where X and Y are the new mesh coordinates and ∇2 is the Laplacian operator defined on the mesh prior to deformation. Along any boundary we can either fix X(Y ) which is Dirichlet boundary condition or ∇X(Y ) · ~n where ~n is the normal to the unperturbed domain. The latter is a Neumann boundary condition. A typical example of this is shown in Figures 1 and 2. The first figure shows an unstructured mesh discretization of the domain x, y ∈ [0, 1] × [−1, 0]. We want to find a deformed mesh for the same domain except with the top surface defined by y = 0.1 cos(2πx). To find the deformed mesh, for X we apply Dirichlet boundary conditions along all boundaries with X = x. Since x is a solution to equation (1) which satisfies these boundary conditions, there is no change in the x position of any of the mesh points. For Y , on the top surface we set Y = 0.1 cos(2πx) which is a Dirichlet boundary condition. On the left and right boundaries periodic conditions are enforced. Along the bottom boundary, a Dirichlet boundary condition is enforced with Y = −1. This results in the deformed mesh shown in Figure 2. c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
MESH DEFORMATION USING THE BIHARMONIC OPERATOR
3
Examining Figure 2, we see that there is a large variation of the normal mesh spacing along the top surface. For most flow calculations, this variation is undesirable. However, because it is impossible to control both the boundary position and the normal mesh spacing, there is no way to avoid this problem when using a second-order PDE. For this reason, we propose using the biharmonic operator. ∇4 X = 0 (2) ∇4 Y = 0 Because these are fourth-order PDE’s, we can specify two boundary conditions along each boundary of the domain. To determine the appropriate combinations, we examine the weak formulation of Equation 2. Since the equations for X and Y are uncoupled, we examine only the equation for X. Multiplying by a test function φ, integrating over the domain, and integrating by parts twice, we arrive at Z Z Z 2 2 3 ∇ φ∇ XdΩ + φ∇ X · ~ndΓ − ∇φ · ~n∇2 XdΓ = 0 (3) Ω
Γ
Γ
where Ω is the domain and Γ is the boundary of the domain. The first integral over Γ reveals that we can either fix X on the boundary or ∇3 X · ~n. Similarly from the second integral we see we can either fix ∇X · ~n or ∇2 X. Thus, there are four possible combinations. However, the combination of ∇2 X and ∇3 X · ~n is not well posed. This can be seen by writing the equation ∇4 X = 0 as −∇2 X = S (4) −∇2 S = 0 Fixing ∇3 X · ~n and ∇2 X is then equivalent to fixing S and ∇S · ~n. The equation for S and its boundary conditions are then uncoupled from that for X. From this we see that using this combination of boundary conditions is an attempt to set both the value and the normal derivative on the boundary for a second-order PDE. The three possible boundary conditions are then to fix X and ∇2 X; ∇X · ~n and ∇3 X · ~n; or X and ∇X · ~n. With homogeneous boundary conditions for ∇2 X and ∇3 X · ~n, use of the first two of these boundary conditions will give exactly the same results as the Laplacian. This can be seen by examining Equation 4. With homogeneous boundary conditions, S will everywhere be zero, leaving ∇2 X = 0. However, with the biharmonic operator, we now have the option of a third boundary condition which is to specify both the position of the boundary and the normal spacing. At interfaces between different meshes or at periodic boundaries, several variants of the above three conditions arise. For example at an x-periodic boundary, the x positions of the vertices must be fixed, but the y positions can be allowed to change. In addition, we would like the mesh spacing to be smooth across the periodic boundary. To accomplish this, we fix the values of X on both the left and right periodic boundaries and enforce continuity of ∇X · ~n and ∇2 X. This fixes the x position of the boundaries but ensures that the normal mesh spacing is continuous across the boundary. For Y , we enforce continuity of Y , ∇Y · ~n, ∇2 Y and ∇Y 3 · ~n. This further ensures mesh continuity across the boundary while allowing the vertical position of the vertices on the boundaries to change. For both X and Y , we have four boundary conditions for the two boundaries which is consistent. c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
4
B. T. HELENBROOK
For an interface between different meshes, we have two possibilities. We can either fix (X, Y ) and ∇(X, Y ) · ~n independently on either side of the interface or we can fix (X, Y ) and enforce continuity of ∇(X, Y ) · ~n and ∇2 (X, Y ). In the first case, the normal mesh spacing on both meshes is fixed at the interface. In the second, the normal mesh spacing can change but is continuous across the interface. The advantages and disadvantages of these are examined in the Results section.
Numerical Implementation The numerical implementation of these equations follows the form given in Equation 4. We demonstrate the implementation for a two-dimensional unstructured mesh of triangles, but the method could just as easily be applied to structured meshes and in three dimensions. The implementation follows a Galerkin finite-element procedure. The discrete equations for X are given by Z Z Z ∇φ · ∇XdΩ − φ∇X · ~ndΓ = φSdΩ (5) Ω
Γ
Z
Ω
Z ∇φ · ∇SdΩ −
Ω
φ∇S · ~ndΓ = 0
(6)
Γ
where φ, X, and S are all approximated using the standard piecewise-linear and continuous space of basis functions. To simplify the formulation, the integral on the right hand side of the first equation is “mass-lumped” (off-diagonal matrix entries are added to the main diagonal and then set to zero). This allows S to be determined directly from X with no matrix inversion. Similar equations apply for Y . The boundary conditions X fixed and ∇2 X = 0 are enforced by setting the values of X and S at the vertices on the boundary. To apply the boundary condition, ∇3 X · ~n = 0 which is equivalent to ∇S · ~n = 0 we set the boundary flux term in Equation 6 to zero. The boundary condition ∇X · ~n fixed is applied in an indirect manner by enforcing the condition that the initial mesh with undeformed boundaries exactly satisfy the equations: On the undeformed mesh, S is calculated from Equation 5 assuming ∇X · ~n is zero. The error in the solution to Equation 6 is then calculated. This negative of the error is treated as a source term in Equation 6 such that the initial, undeformed mesh exactly satisfies the equations. In this way, we implicitly determine appropriate boundary fluxes for ∇X · ~n.
Solution Method For a periodic one-dimensional domain, the eigenvalues of the biharmonic equations scale as k 4 where k is the wavenumber of the eigenmode. For the Laplacian equations, the eigenvalues scale as k 2 . Thus for a given range of possible wavenumbers as defined by the mesh resolution, the biharmonic equations are a stiffer system than the Laplacian equations. It is well known that a multigrid iterative method can avoid the difficulties of stiffness associated with the eigenvalue distribution of the Laplacian equations. A multigrid method is also an effective approach for obtaining solutions to the biharmonic equations. c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
MESH DEFORMATION USING THE BIHARMONIC OPERATOR
5
We use geometric multigrid defined by a sequence of nested unstructured meshes. Each successively coarser mesh has approximately one half the mesh resolution of the previous mesh in any direction. The method of creating these meshes given an initial mesh is described in [10]. The interpolation between meshes is performed using linear operators. This is described in [11]. For a relaxation scheme, a Jacobi iteration is used. This is implemented as follows: Given an initial guess for X, S is calculated from Equation 5. This is then inserted into Equation 6 with the constant source term added to indirectly apply the ∇X · ~n boundary conditions. The error in the solution is used to drive the change in the guess for X as follows ∆Xi =
C Ei Di
(7)
where ∆Xi is the change to the initial guess at any vertex i and Ei is the error in the solution associated with that vertex. C is a relaxation factor between 0 and 1, and Di is the diagonal component of the matrix form of equations 5 and 6 after the S variables have been eliminated from the system. This is given by " # PNj Ni 2 X kij k jk Di = + kij k=1 (8) Vj Vj j=1 where Ni is the number of sides connected to vertex i and Vi is 1/3 the area of all triangles connected to vertex i. kij is the matrix entry at row i column j associated with the operator ∇2 as discretized by equations 5 or 6. Since the Laplacian operator is symmetric, kij = kji . Also, kij is only non-zero for vertices which are connected by a triangle side. The values of kij for each side can be calculated as follows:
c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
6
B. T. HELENBROOK
for every side: kij = 0 for every triangle: for every side of triangle: kij -= length of side / area of triangle two opposing sides: kij (opposing) += length of side / area of triangle endloop endloop
The error Ei is calculated using the kij by looping over the sides. This is more efficient than doing an element by element integration [1]. Since this has to be done twice, once to calculate S and once to calculate E, one Jacobi iteration for the biharmonic is approximately twice as expensive as a Jacobi iteration for the Laplacian.
Results We begin by re-examining the example problem shown in Figures 1 and 2. The boundary conditions for X are X and ∇X · ~n fixed on the top surface. On the left and right boundaries, X is fixed and ∇X · ~n and ∇2 X are forced to be periodic. On the bottom surface we fix X and set ∇2 X = 0. For Y , on the top surface we fix Y and ∇Y · ~n. On the periodic boundaries, we enforce continuity of Y , ∇Y · ~n, ∇2 Y , and ∇3 Y · ~n. On the bottom boundary, we fix Y and set ∇2 Y = 0. These boundary conditions guarantee that there is no jump in mesh spacing across the periodic boundaries. This would not be possible with a second-order PDE because we would not be able to control both X and ∇X · ~n. Similarly on the top boundary, we can now control both the boundary position and the normal mesh spacing. Figure 3 shows the results of the mesh deformation when the free surface is perturbed by Y = 0.1 cos(2πx). There is very little variation in normal mesh spacing along the top surface which is a significant improvement over the Laplacian result shown in Figure 2. The importance of this increases as the top surface becomes more deformed. Figures 4 and 5 show results for a perturbation amplitude of 0.25. In this case, the Laplacian method allows mesh lines to cross and negative triangle areas to occur. The biharmonic method produces a usable mesh. However, there is some variation in normal mesh spacing along the top surface. This occurs because there is a large change in normal direction of the boundary relative to the initial mesh. At the top boundary, ∇X ·~n is near zero on the initial mesh and ∇Y ·~n determines the normal mesh spacing. As the boundary rotates from horizontal, the normal mesh spacing reduces because an increasing weight is given to ∇X · ~n in determining the normal mesh spacing. In the next example we show that for most flow calculations, this variation in normal spacing can easily be eliminated. In the above test, the deformation from the initial mesh to the final mesh is accomplished in one step. In most simulations however, the deformation occurs over many time steps. To see whether this improves the final mesh, we perform a simulation with the same initial and final states as in the previous calculations, however the change in the upper boundary is spread over 10 steps. The upper boundary position at each intermediate step is given by y = 0.25 i/10 cos(2πx) where i is the step number. At each step, the Laplacian constants kij c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
MESH DEFORMATION USING THE BIHARMONIC OPERATOR
7
are reevaluated and the implicit process of determining ∇X(Y ) · ~n is repeated. Figure 6 shows the result when this is done using a Laplacian method and Figure 7 shows the results for the biharmonic method. On the last deformation interval, the Laplacian method allows mesh lines to cross and negative volumes to occur. The biharmonic method produces a mesh with smoothly varying triangles of good quality. By dividing the boundary movement into intervals, we also eliminate the variation in normal mesh spacing which is seen in Figure 5. To determine the relative cost of the Laplacian method versus the biharmonic method, we examine the multigrid convergence rates for the previous problem. A w-cycle [11] multigrid iteration is used for both methods with 1 Jacobi iteration on each mesh before coarsening to the next mesh. No iterations are done after moving to a finer mesh. Six multigrid levels are used such that the coarsest mesh consists of two triangles. For both methods, we use a relaxation factor of 0.9. For this problem, this gives slightly better results than adding the full correction. Figure 8 shows the maximum change in position for any vertex after a fine grid Jacobi iteration versus the number of w-cycle iterations. Clearly the Laplacian equations converge faster. How much this affects the efficiency depends on how the method is used; there is no need to solve the mesh-movement equations to high precision because the mesh positions are somewhat arbitrary. To determine the minimum number of w-cycle iterations needed, we repeat the prior test in which the upper boundary deformation is spread over ten intervals. For the biharmonic method approximately 20 w-cycles each interval are needed to produce a mesh visually identical to Figure 7. For the Laplacian method, 10 w-cycles per interval are needed to reproduce Figure 6. Including the fact that a Jacobi iteration of the biharmonic method is twice as expensive as the Laplacian method, we see that the biharmonic method is around 4 times the cost of the Laplacian method. However, in many cases this may be offset by the fact that the cost of the mesh-movement is smaller than the flow solver. Such is the case in spectral/hp simulations [2, 12] where the flow resolution is much greater than the number of elements in the mesh. For implicit finite-volume solvers, this cost can also be negligible because it is more difficult to solve the implicit Navier-Stokes equations [4]. Furthermore, because the meshes are of better quality than with the Laplacian method, this method will reduce the need to restructure the mesh which can improve the accuracy and efficiency of a calculation. One of the main motivations for trying the biharmonic method is two-fluid problems. For these problems, the mesh tracks a moving interface across which one would like to ensure that the normal mesh spacing is continuous. As discussed previously, this can be accomplished in two ways: The first is to enforce constraints on the position (X, Y ) and the normal spacing ∇(X, Y ) · ~n of each mesh separately. To make the mesh continuous and smooth across the interface, these constraints must be the same on either side of the interface. The second is to fix the interface position (X, Y ) and enforce continuity of both ∇(X, Y ) · ~n and ∇2 X. The first method directly determines the mesh spacing near the interface for both meshes. The second method only guarantees that the mesh spacing will be continuous. For most two-fluid problems, there is usually a boundary-layer along the interface which needs to be resolved. Thus, it is usually most appropriate to directly control the normal mesh spacing. To demonstrate the advantage of the biharmonic method for interface problems, we examine the mesh movement associated with a simple rising surface. The initial mesh is as shown in Figure 1 with an upper mesh reflected over the x-axis. We then raise the midplane from 0 to 0.5. Figure 9 shows the top half of the mesh for the Laplacian method. In this case there is a discontinuity in mesh spacing across the interface. The result obtained using the biharmonic c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
8
B. T. HELENBROOK
method is shown in Figure 10. The biharmonic method produces a smooth gradation in mesh spacing at the interface. For this problem, either of the two different interface boundary conditions produce identical results. This can be shown with a simple one-dimensional analysis. One drawback of the biharmonic method is that even in one dimension, Y is not bounded by the values on the boundary of the domain. In the previous problem, if we increase the height of the interface from 0.5 to 0.75, the solution for Y in the upper domain increases from 0.75 at the interface to values greater than one then decreases to the boundary value of one. With the second-order method, the values of Y are always between the minimum and maximum values on the boundary. In spite of this, the biharmonic method produces better meshes for most problems involving interface deformation as shown in the next example. Figure 11 shows results for a deforming interface with the two different interface boundary conditions. The initial mesh is the same as described in the previous problem. The final interface position is given by y = 3/8 cos(2πx), and the interface is deformed over 15 steps. The left half of Figure 11 shows the resulting mesh for the boundary conditions (X, Y ) and ∇(X, Y ) · ~n fixed separately on either side of the interface. (The result is reflected over the y-axis.) The right half shows the resulting mesh for (X, Y ) fixed and ∇(X, Y ) ·~n and ∇2 (X, Y ) continuous across the interface. Both meshes are fairly similar; Because of the stretching of the triangles near the interface, there is a significant compression of the triangles at the upper and lower boundaries. For this problem, it is not clear which boundary condition performs better, however we still favor the first boundary condition for the reasons discussed above. We also note that this calculation could not be completed with the Laplacian method because negative triangle areas occur. The problem of a moving cylinder in a channel was used as a test problem in [9] for an elastic deformation method. The initial mesh is shown in Figure 12. For the first calculation, we fix (X, Y ) on the boundary of the channel. We then move the cylinder one diameter to the right. Figure 13 shows the result for the Laplacian method. On the right half of the cylinder, the mesh lines are crossed and there are triangles with negative areas. Figure 14 shows the result for the biharmonic method. On the boundary of the channel ∇2 (X, Y ) = 0 and on the cylinder surface ∇(X, Y ) · ~n is fixed. In this case, the mesh remains well connected. Also note that the cells on the surface of the cylinder remain relatively undistorted. As a final test, we examine the problem of a moving cylinder with a refined boundary layer. Because of it’s relevance in fluid dynamics, similar problems have been examined by many previous authors [1, 4, 5, 13]. Figure 15 shows the initial mesh with an inset showing the refined region around the cylinder. We then move the cylinder 4.5 diameters to the left and upward and solve for the deformed mesh. Figure 16 shows the mesh after being deformed by the biharmonic method. The boundary conditions used are the same as for the previous problem of a cylinder in a channel. From the inset, we see that the mesh near the cylinder remains almost unchanged. The far-field mesh distorts, but none of the triangles develop high aspect ratios. With the Laplacian method, mesh lines cross and negative volumes occur so the problem cannot be completed. To fix this problem, several authors have proposed a modification to the Laplace method in which a variable diffusivity is introduced. Near the cylinder the diffusivity is increased making the mesh stiffer. There are several different implementations of this technique [1, 3, 4, 8]. In several of the studies it was found that increasing the diffusivity as the inverse of the square of the mesh length scale was the most effective for boundary layers [4, 8]. To compare the biharmonic method to this approach, we implement the “spring” method used in c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
MESH DEFORMATION USING THE BIHARMONIC OPERATOR
9
[3, 4]: Each side of the mesh is treated as a spring with a spring constant inversely proportional to the square of it’s length. Figure 17 shows the deformed mesh calculated using the spring method. For this problem, the spring method is as effective as the biharmonic method. The large variable diffusivity in the modified Laplacian/spring methods increases the stiffness of the equations to be solved compared to the standard Laplacian method. For the above problem however, the iteration convergence rate with multigrid and a Jacobi iteration is still faster than for the biharmonic equations. In addition, since the biharmonic equations are roughly twice as expensive as the modified Laplacian/spring approach, there is no reason to use the biharmonic for this problem. However for problems without a significant gradation in the mesh length scale (such as for most of the examples in this paper), the modified Laplacian/spring approaches will perform as poorly as the Laplacian method. Furthermore, for problems with boundary deformation in a localized region of near uniform mesh resolution (typical of an ALE two-phase flow computation), the biharmonic method should perform better.
Conclusions The ability to control both the position of boundaries and the normal mesh spacing is a significant advantage of the biharmonic method over second-order methods for deforming a mesh. This not only allows us to guarantee that mesh spacings will be continuous at interfacial or periodic boundaries, but it also allows the simulation of problems with greater boundary deformation than second-order methods without restructuring the mesh. This is valuable for both free-surface and interfacial flow calculations. Furthermore for problems with boundary layers, the biharmonic method can help preserve the spacing of the mesh in the boundary-layer so that restructuring is only needed in regions exterior to the boundary layer. This can also be accomplished by a variable diffusivity Laplacian method, which is computationally cheaper. However for problems of multiple bodies with interacting boundary layers or with localized boundary deformation, the biharmonic method may perform better. The main disadvantage of the biharmonic method is the cost. Using multigrid to solve the system of equations, the Laplacian method is about 1/4 the cost of the biharmonic method. Thus, the biharmonic method will not be optimal for all problems. However, for simulations in which the mesh movement is a small fraction of the the total cost, this may not be important. Furthermore, the fact that the need for mesh restructuring is reduced may also outweigh the increased cost.
ACKNOWLEDGEMENTS
Mesh generation for the example problems was done using EasyMesh which is free software written by Bojan Niceno. This work was supported in part by the Department of Energy under the ASCI/ASAP contract to Stanford University. Opinions expressed here are those of the authors and not of Stanford University or the Department of Energy.
REFERENCES c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
10
B. T. HELENBROOK
1. R. L¨ ohner and C. Yang. Improved ALE mesh velocities for moving boundaries. Comm. Num. Meth. Eng., 12:599–608, 1996. 2. I. Robertson and S. Sherwin. Free-surface flow simulation using hp/spectral elements. J. Comp. Phys., 155:26–53, 1999. 3. J. Batina. Unsteady euler algoithm with unstructured dynamic mesh for complex-aircraft aerodynamic analysis. AIAA Journal, 29(3):327–333, 1991. 4. V. Venkatakrishnan and D. J. Mavriplis. Implicit method for the computation of unsteady flows on unstructured grids. AIAA, 95-1705, 1995. 5. O. Hassan, E. J. Probert, and K. Morgan. Unstructured mesh procedures for the simulation of three dimensional transient compressible inviscid flows with moving boundary components. Int. J. Num. Meth. Fluids, 27:41–55, 1998. 6. C. Farhat, C. Degand, B. Koobus, and M. Lesoinee. Torsional springs for two-dimensional dynamic unstructured fluid meshes. Comp. Meth. Appl. Mech. Engrg., 163:231–245, 1998. 7. D. R. Lynch. Unified approach to simulation on deforming elements with application to phase change problems. J. Comput. Phys., 47:387–411, 1982. 8. T. E. Tezduyar, M. Behr, S. Mittal, and A. A. Johnson. Computation of unsteady incompressible flows with the stabilized finite element methods–space-time formulations, iterative strategies and massively parallel implementations. In P. Smolinksi, editor, New Methods in Transient Analysis, pages 7–24. ASME, New York, NY, 1992. 9. T. J. Baker. Mesh modification for solution adaptation and time evolving domains. In 7th International Conference on Numerical Grid Generation in Computational Field Simulations, 2000. 10. P. Lin, L. Martinelli, T. J. Baker, and A. Jameson. Two-dimensional implicit time dependent calculations for incompressible flows on adaptive unstructured meshes. AIAA, 01-2655, 2001. 11. D. J. Mavriplis. Multigrid techniques for unstructured meshes. Technical Report 95-27, ICASE, April 1995. 12. B. T. Helenbrook. A two-fluid spectral element method. Comp. Meth. Appl. Mech. Eng., 191:273–294, 2001. 13. T. E. Tezduyar, M. Behr, and J. Liou. A new strategy for finite element computations involving moving boundaries and interfaces—the deforming-spatial-domain / space-time procedure: Ii. computation of freesurface flows, two-liquid flows, and flows with drifting cylinders. Comp. Meth. App. Mech. Eng., 94:353– 371, 1992.
c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
MESH DEFORMATION USING THE BIHARMONIC OPERATOR
11
0.1 0 -0.1 -0.2 -0.3
y
-0.4 -0.5 -0.6 -0.7 -0.8 -0.9 -1
0
0.25
0.5
0.75
1
x Figure 1. Initial free-surface mesh
c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
12
B. T. HELENBROOK
0.1 0 -0.1 -0.2 -0.3
y
-0.4 -0.5 -0.6 -0.7 -0.8 -0.9 -1
0
0.25
0.5
0.75
1
x Figure 2. Free-surface mesh deformed with the Laplacian operator. Amplitude = 0.1
c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
MESH DEFORMATION USING THE BIHARMONIC OPERATOR
13
0.1 0 -0.1 -0.2 -0.3
y
-0.4 -0.5 -0.6 -0.7 -0.8 -0.9 -1
0
0.25
0.5
0.75
1
x Figure 3. Free-surface mesh deformed with the biharmonic operator. Amplitude = 0.1
c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
14
B. T. HELENBROOK
0.2 0.1 0 -0.1 -0.2
y
-0.3 -0.4 -0.5 -0.6 -0.7 -0.8 -0.9 -1
0
0.25
0.5
0.75
1
x Figure 4. Free-surface mesh deformed with the Laplacian operator. Amplitude = 0.25
c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
MESH DEFORMATION USING THE BIHARMONIC OPERATOR
15
0.2 0.1 0 -0.1 -0.2
y
-0.3 -0.4 -0.5 -0.6 -0.7 -0.8 -0.9 -1
0
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
1
x Figure 5. Free-surface mesh deformed with the biharmonic operator. Amplitude = 0.25
c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
16
B. T. HELENBROOK
0.2 0.1 0 -0.1 -0.2
y
-0.3 -0.4 -0.5 -0.6 -0.7 -0.8 -0.9 -1
0
0.25
0.5
0.75
1
x Figure 6. Free-surface mesh deformed with the Laplacian operator over ten intermediate steps.
c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
MESH DEFORMATION USING THE BIHARMONIC OPERATOR
17
0.2 0.1 0 -0.1 -0.2
y
-0.3 -0.4 -0.5 -0.6 -0.7 -0.8 -0.9 -1
0
0.25
0.5
0.75
1
x Figure 7. Free-surface mesh deformed with the biharmonic operator over ten intermediate steps.
c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
18
Maximum Correction to Vertex Position
B. T. HELENBROOK
10
-1
10
-2
10-3
Biharmonic 10-4 10
-5
10
-6
10-7
Laplacian 10-8 0
20
40
60
80
w-cycles Figure 8. Multigrid convergence rates for the Laplacian and Biharmonic methods.
c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
MESH DEFORMATION USING THE BIHARMONIC OPERATOR
19
1
y
0.75
0.5
0.25
0 0
0.25
0.5
0.75
1
x Figure 9. Laplacian method for deforming a mesh with a rising interface
c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
20
B. T. HELENBROOK
1
y
0.75
0.5
0.25
0 0
0.25
0.5
0.75
1
x Figure 10. Biharmonic method for deforming a mesh with a rising interface
c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
MESH DEFORMATION USING THE BIHARMONIC OPERATOR
21
1
0.75
y
0.5
0.25
0
-0.25
-0.5
0
0.5
x Figure 11. Interface deformation using the biharmonic method. Left: upper and lower meshes are treated separately. Right: mesh continuity enforced.
c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
22
B. T. HELENBROOK
y
1 0
-1
0
2
x 4
6
8
Figure 12. Undeformed mesh for cylinder movement in a channel
c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
23
MESH DEFORMATION USING THE BIHARMONIC OPERATOR
y
1 0
-1
0
2
x 4
6
8
Figure 13. Deformed mesh calculated with the Laplacian method for cylinder movement in a channel
c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
24
B. T. HELENBROOK
y
1 0
-1
0
2
x 4
6
8
Figure 14. Deformed mesh calculated with the biharmonic method for cylinder movement in a channel
c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
25
MESH DEFORMATION USING THE BIHARMONIC OPERATOR
y
10
0
-10
-10
0
x
10
20
Figure 15. Undeformed mesh for cylinder movement with a boundary layer. The inset shows the mesh around the cylinder
c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
26
B. T. HELENBROOK
y
10
0
-10
-10
0
x
10
20
Figure 16. Deformed mesh calculated using the biharmonic method for cylinder movement with a boundary layer. The inset shows the mesh around the cylinder
c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
27
MESH DEFORMATION USING THE BIHARMONIC OPERATOR
y
10
0
-10
-10
0
x
10
20
Figure 17. Deformed mesh calculated using the spring method for cylinder movement with a boundary layer. The inset shows the mesh around the cylinder
c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
Mesh deformation using the biharmonic operator Brian T. Helenbrook
∗
Mechanical and Aeronautical Engineering Department Clarkson University Potsdam, NY 13699-5725
SUMMARY The use of the biharmonic operator for deforming a mesh in an arbitrary-Lagrangian-Eulerian simulation is investigated. The biharmonic operator has the advantage that two conditions can be specified on each boundary of the mesh. This allows both the position and the normal mesh spacing along a boundary to be controlled, which is important for two-fluid interfaces and periodic boundaries. At these boundaries, we can simultaneously fix the position of the boundary and ensure that the normal mesh spacing is continuous across the boundary. In addition, results for deforming surfaces show that greater surface deformation can be tolerated when using biharmonic equations compared to approaches using second-order partial differential equations. A final advantage is that with the biharmonic operator, the integrity of a grid in a moving boundary layer can be preserved as the boundary moves. The main disadvantage of the approach is it’s increased computational expense. c 2001 John Wiley & Sons, Ltd. Copyright key words: adaptation
arbitrary-Lagrangian-Eulerian (ALE), mesh deformation, mesh movement, r-
Introduction In a typical arbitrary-Lagrangian-Eulerian (ALE) flow simulation, a boundary of a mesh tracks some moving/deforming feature of interest. This is usually a Lagrangian feature such as a cylinder surface or a free-surface wave. As the boundary of the mesh moves, the interior region of the mesh is deformed “arbitrarily” to maintain a quality mesh for a Eulerian flow solution. Because the mesh deformation is somewhat arbitrary, there have been several different methods proposed for doing this. These methods typically compromise between the amount of boundary deformation that can be tolerated and the computational expense of the procedure.
∗ Correspondence
to: Mechanical and Aeronautical Engineering Department Clarkson University Potsdam, NY 13699-5725 e-mail: [email protected] phone: 315-268-2204 fax: 315-268-6438
c 2001 John Wiley & Sons, Ltd. Copyright
2
B. T. HELENBROOK
Most methods use some form of a second-order partial differential equation (PDE). Typical examples are the Laplacian method [1, 2], the spring method [3, 4, 5, 6], and the elastic deformation method [7, 8, 9]. A limitation of a second-order method is that at boundaries either the mesh position or the normal mesh spacing can be specified, but it is impossible to specify both. A typical example where we might want to control both is a moving cylinder. At the surface of the cylinder, we would like to control the position of the boundary as well as the boundary-layer resolution i.e. the normal mesh spacing. With a second-order PDE, this is impossible. Another example is the simulation of two-fluid flows. In this case, the edges of two meshes meet to define the position of an interface between two different fluids. This mesh boundary is deformed to track the position of the interface and thus is a Dirichlet boundary condition for the mesh deformation PDE. However, we would also like to guarantee that the two meshes maintain a smooth gradation in mesh spacing across the interface. This requires that we specify both the position and the normal spacing at the interface boundary, which is impossible with a second-order PDE. To overcome this limitation, the use of a fourth-order PDE is proposed for determining the interior mesh deformation. Fourth-order generalizations of any of the three second-order deformation methods can be created, however this work focuses on a generalization of the Laplacian method because it is the most straightforward. The first section of this paper describes the use of the biharmonic operator as a fourth-order generalization of the Laplacian operator for mesh deformation. In this section, the appropriate choice of boundary conditions for modeling different types of boundaries in an ALE flow simulation is also discussed. The following section presents a multigrid solution method and examines the relative cost of the biharmonic versus the Laplacian methods. Results are obtained using both methods for several cases including a deforming free-surface, a deforming two-fluid interface, a moving cylinder in a channel, and a moving cylinder with a refined boundary layer mesh.
Governing Equations A typical method used to deform a mesh is the following system of Laplacian equations ∇2 X = 0 ∇2 Y = 0
(1)
where X and Y are the new mesh coordinates and ∇2 is the Laplacian operator defined on the mesh prior to deformation. Along any boundary we can either fix X(Y ) which is Dirichlet boundary condition or ∇X(Y ) · ~n where ~n is the normal to the unperturbed domain. The latter is a Neumann boundary condition. A typical example of this is shown in Figures 1 and 2. The first figure shows an unstructured mesh discretization of the domain x, y ∈ [0, 1] × [−1, 0]. We want to find a deformed mesh for the same domain except with the top surface defined by y = 0.1 cos(2πx). To find the deformed mesh, for X we apply Dirichlet boundary conditions along all boundaries with X = x. Since x is a solution to equation (1) which satisfies these boundary conditions, there is no change in the x position of any of the mesh points. For Y , on the top surface we set Y = 0.1 cos(2πx) which is a Dirichlet boundary condition. On the left and right boundaries periodic conditions are enforced. Along the bottom boundary, a Dirichlet boundary condition is enforced with Y = −1. This results in the deformed mesh shown in Figure 2. c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
MESH DEFORMATION USING THE BIHARMONIC OPERATOR
3
Examining Figure 2, we see that there is a large variation of the normal mesh spacing along the top surface. For most flow calculations, this variation is undesirable. However, because it is impossible to control both the boundary position and the normal mesh spacing, there is no way to avoid this problem when using a second-order PDE. For this reason, we propose using the biharmonic operator. ∇4 X = 0 (2) ∇4 Y = 0 Because these are fourth-order PDE’s, we can specify two boundary conditions along each boundary of the domain. To determine the appropriate combinations, we examine the weak formulation of Equation 2. Since the equations for X and Y are uncoupled, we examine only the equation for X. Multiplying by a test function φ, integrating over the domain, and integrating by parts twice, we arrive at Z Z Z 2 2 3 ∇ φ∇ XdΩ + φ∇ X · ~ndΓ − ∇φ · ~n∇2 XdΓ = 0 (3) Ω
Γ
Γ
where Ω is the domain and Γ is the boundary of the domain. The first integral over Γ reveals that we can either fix X on the boundary or ∇3 X · ~n. Similarly from the second integral we see we can either fix ∇X · ~n or ∇2 X. Thus, there are four possible combinations. However, the combination of ∇2 X and ∇3 X · ~n is not well posed. This can be seen by writing the equation ∇4 X = 0 as −∇2 X = S (4) −∇2 S = 0 Fixing ∇3 X · ~n and ∇2 X is then equivalent to fixing S and ∇S · ~n. The equation for S and its boundary conditions are then uncoupled from that for X. From this we see that using this combination of boundary conditions is an attempt to set both the value and the normal derivative on the boundary for a second-order PDE. The three possible boundary conditions are then to fix X and ∇2 X; ∇X · ~n and ∇3 X · ~n; or X and ∇X · ~n. With homogeneous boundary conditions for ∇2 X and ∇3 X · ~n, use of the first two of these boundary conditions will give exactly the same results as the Laplacian. This can be seen by examining Equation 4. With homogeneous boundary conditions, S will everywhere be zero, leaving ∇2 X = 0. However, with the biharmonic operator, we now have the option of a third boundary condition which is to specify both the position of the boundary and the normal spacing. At interfaces between different meshes or at periodic boundaries, several variants of the above three conditions arise. For example at an x-periodic boundary, the x positions of the vertices must be fixed, but the y positions can be allowed to change. In addition, we would like the mesh spacing to be smooth across the periodic boundary. To accomplish this, we fix the values of X on both the left and right periodic boundaries and enforce continuity of ∇X · ~n and ∇2 X. This fixes the x position of the boundaries but ensures that the normal mesh spacing is continuous across the boundary. For Y , we enforce continuity of Y , ∇Y · ~n, ∇2 Y and ∇Y 3 · ~n. This further ensures mesh continuity across the boundary while allowing the vertical position of the vertices on the boundaries to change. For both X and Y , we have four boundary conditions for the two boundaries which is consistent. c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
4
B. T. HELENBROOK
For an interface between different meshes, we have two possibilities. We can either fix (X, Y ) and ∇(X, Y ) · ~n independently on either side of the interface or we can fix (X, Y ) and enforce continuity of ∇(X, Y ) · ~n and ∇2 (X, Y ). In the first case, the normal mesh spacing on both meshes is fixed at the interface. In the second, the normal mesh spacing can change but is continuous across the interface. The advantages and disadvantages of these are examined in the Results section.
Numerical Implementation The numerical implementation of these equations follows the form given in Equation 4. We demonstrate the implementation for a two-dimensional unstructured mesh of triangles, but the method could just as easily be applied to structured meshes and in three dimensions. The implementation follows a Galerkin finite-element procedure. The discrete equations for X are given by Z Z Z ∇φ · ∇XdΩ − φ∇X · ~ndΓ = φSdΩ (5) Ω
Γ
Z
Ω
Z ∇φ · ∇SdΩ −
Ω
φ∇S · ~ndΓ = 0
(6)
Γ
where φ, X, and S are all approximated using the standard piecewise-linear and continuous space of basis functions. To simplify the formulation, the integral on the right hand side of the first equation is “mass-lumped” (off-diagonal matrix entries are added to the main diagonal and then set to zero). This allows S to be determined directly from X with no matrix inversion. Similar equations apply for Y . The boundary conditions X fixed and ∇2 X = 0 are enforced by setting the values of X and S at the vertices on the boundary. To apply the boundary condition, ∇3 X · ~n = 0 which is equivalent to ∇S · ~n = 0 we set the boundary flux term in Equation 6 to zero. The boundary condition ∇X · ~n fixed is applied in an indirect manner by enforcing the condition that the initial mesh with undeformed boundaries exactly satisfy the equations: On the undeformed mesh, S is calculated from Equation 5 assuming ∇X · ~n is zero. The error in the solution to Equation 6 is then calculated. This negative of the error is treated as a source term in Equation 6 such that the initial, undeformed mesh exactly satisfies the equations. In this way, we implicitly determine appropriate boundary fluxes for ∇X · ~n.
Solution Method For a periodic one-dimensional domain, the eigenvalues of the biharmonic equations scale as k 4 where k is the wavenumber of the eigenmode. For the Laplacian equations, the eigenvalues scale as k 2 . Thus for a given range of possible wavenumbers as defined by the mesh resolution, the biharmonic equations are a stiffer system than the Laplacian equations. It is well known that a multigrid iterative method can avoid the difficulties of stiffness associated with the eigenvalue distribution of the Laplacian equations. A multigrid method is also an effective approach for obtaining solutions to the biharmonic equations. c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
MESH DEFORMATION USING THE BIHARMONIC OPERATOR
5
We use geometric multigrid defined by a sequence of nested unstructured meshes. Each successively coarser mesh has approximately one half the mesh resolution of the previous mesh in any direction. The method of creating these meshes given an initial mesh is described in [10]. The interpolation between meshes is performed using linear operators. This is described in [11]. For a relaxation scheme, a Jacobi iteration is used. This is implemented as follows: Given an initial guess for X, S is calculated from Equation 5. This is then inserted into Equation 6 with the constant source term added to indirectly apply the ∇X · ~n boundary conditions. The error in the solution is used to drive the change in the guess for X as follows ∆Xi =
C Ei Di
(7)
where ∆Xi is the change to the initial guess at any vertex i and Ei is the error in the solution associated with that vertex. C is a relaxation factor between 0 and 1, and Di is the diagonal component of the matrix form of equations 5 and 6 after the S variables have been eliminated from the system. This is given by " # PNj Ni 2 X kij k jk Di = + kij k=1 (8) Vj Vj j=1 where Ni is the number of sides connected to vertex i and Vi is 1/3 the area of all triangles connected to vertex i. kij is the matrix entry at row i column j associated with the operator ∇2 as discretized by equations 5 or 6. Since the Laplacian operator is symmetric, kij = kji . Also, kij is only non-zero for vertices which are connected by a triangle side. The values of kij for each side can be calculated as follows:
c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
6
B. T. HELENBROOK
for every side: kij = 0 for every triangle: for every side of triangle: kij -= length of side / area of triangle two opposing sides: kij (opposing) += length of side / area of triangle endloop endloop
The error Ei is calculated using the kij by looping over the sides. This is more efficient than doing an element by element integration [1]. Since this has to be done twice, once to calculate S and once to calculate E, one Jacobi iteration for the biharmonic is approximately twice as expensive as a Jacobi iteration for the Laplacian.
Results We begin by re-examining the example problem shown in Figures 1 and 2. The boundary conditions for X are X and ∇X · ~n fixed on the top surface. On the left and right boundaries, X is fixed and ∇X · ~n and ∇2 X are forced to be periodic. On the bottom surface we fix X and set ∇2 X = 0. For Y , on the top surface we fix Y and ∇Y · ~n. On the periodic boundaries, we enforce continuity of Y , ∇Y · ~n, ∇2 Y , and ∇3 Y · ~n. On the bottom boundary, we fix Y and set ∇2 Y = 0. These boundary conditions guarantee that there is no jump in mesh spacing across the periodic boundaries. This would not be possible with a second-order PDE because we would not be able to control both X and ∇X · ~n. Similarly on the top boundary, we can now control both the boundary position and the normal mesh spacing. Figure 3 shows the results of the mesh deformation when the free surface is perturbed by Y = 0.1 cos(2πx). There is very little variation in normal mesh spacing along the top surface which is a significant improvement over the Laplacian result shown in Figure 2. The importance of this increases as the top surface becomes more deformed. Figures 4 and 5 show results for a perturbation amplitude of 0.25. In this case, the Laplacian method allows mesh lines to cross and negative triangle areas to occur. The biharmonic method produces a usable mesh. However, there is some variation in normal mesh spacing along the top surface. This occurs because there is a large change in normal direction of the boundary relative to the initial mesh. At the top boundary, ∇X ·~n is near zero on the initial mesh and ∇Y ·~n determines the normal mesh spacing. As the boundary rotates from horizontal, the normal mesh spacing reduces because an increasing weight is given to ∇X · ~n in determining the normal mesh spacing. In the next example we show that for most flow calculations, this variation in normal spacing can easily be eliminated. In the above test, the deformation from the initial mesh to the final mesh is accomplished in one step. In most simulations however, the deformation occurs over many time steps. To see whether this improves the final mesh, we perform a simulation with the same initial and final states as in the previous calculations, however the change in the upper boundary is spread over 10 steps. The upper boundary position at each intermediate step is given by y = 0.25 i/10 cos(2πx) where i is the step number. At each step, the Laplacian constants kij c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
MESH DEFORMATION USING THE BIHARMONIC OPERATOR
7
are reevaluated and the implicit process of determining ∇X(Y ) · ~n is repeated. Figure 6 shows the result when this is done using a Laplacian method and Figure 7 shows the results for the biharmonic method. On the last deformation interval, the Laplacian method allows mesh lines to cross and negative volumes to occur. The biharmonic method produces a mesh with smoothly varying triangles of good quality. By dividing the boundary movement into intervals, we also eliminate the variation in normal mesh spacing which is seen in Figure 5. To determine the relative cost of the Laplacian method versus the biharmonic method, we examine the multigrid convergence rates for the previous problem. A w-cycle [11] multigrid iteration is used for both methods with 1 Jacobi iteration on each mesh before coarsening to the next mesh. No iterations are done after moving to a finer mesh. Six multigrid levels are used such that the coarsest mesh consists of two triangles. For both methods, we use a relaxation factor of 0.9. For this problem, this gives slightly better results than adding the full correction. Figure 8 shows the maximum change in position for any vertex after a fine grid Jacobi iteration versus the number of w-cycle iterations. Clearly the Laplacian equations converge faster. How much this affects the efficiency depends on how the method is used; there is no need to solve the mesh-movement equations to high precision because the mesh positions are somewhat arbitrary. To determine the minimum number of w-cycle iterations needed, we repeat the prior test in which the upper boundary deformation is spread over ten intervals. For the biharmonic method approximately 20 w-cycles each interval are needed to produce a mesh visually identical to Figure 7. For the Laplacian method, 10 w-cycles per interval are needed to reproduce Figure 6. Including the fact that a Jacobi iteration of the biharmonic method is twice as expensive as the Laplacian method, we see that the biharmonic method is around 4 times the cost of the Laplacian method. However, in many cases this may be offset by the fact that the cost of the mesh-movement is smaller than the flow solver. Such is the case in spectral/hp simulations [2, 12] where the flow resolution is much greater than the number of elements in the mesh. For implicit finite-volume solvers, this cost can also be negligible because it is more difficult to solve the implicit Navier-Stokes equations [4]. Furthermore, because the meshes are of better quality than with the Laplacian method, this method will reduce the need to restructure the mesh which can improve the accuracy and efficiency of a calculation. One of the main motivations for trying the biharmonic method is two-fluid problems. For these problems, the mesh tracks a moving interface across which one would like to ensure that the normal mesh spacing is continuous. As discussed previously, this can be accomplished in two ways: The first is to enforce constraints on the position (X, Y ) and the normal spacing ∇(X, Y ) · ~n of each mesh separately. To make the mesh continuous and smooth across the interface, these constraints must be the same on either side of the interface. The second is to fix the interface position (X, Y ) and enforce continuity of both ∇(X, Y ) · ~n and ∇2 X. The first method directly determines the mesh spacing near the interface for both meshes. The second method only guarantees that the mesh spacing will be continuous. For most two-fluid problems, there is usually a boundary-layer along the interface which needs to be resolved. Thus, it is usually most appropriate to directly control the normal mesh spacing. To demonstrate the advantage of the biharmonic method for interface problems, we examine the mesh movement associated with a simple rising surface. The initial mesh is as shown in Figure 1 with an upper mesh reflected over the x-axis. We then raise the midplane from 0 to 0.5. Figure 9 shows the top half of the mesh for the Laplacian method. In this case there is a discontinuity in mesh spacing across the interface. The result obtained using the biharmonic c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
8
B. T. HELENBROOK
method is shown in Figure 10. The biharmonic method produces a smooth gradation in mesh spacing at the interface. For this problem, either of the two different interface boundary conditions produce identical results. This can be shown with a simple one-dimensional analysis. One drawback of the biharmonic method is that even in one dimension, Y is not bounded by the values on the boundary of the domain. In the previous problem, if we increase the height of the interface from 0.5 to 0.75, the solution for Y in the upper domain increases from 0.75 at the interface to values greater than one then decreases to the boundary value of one. With the second-order method, the values of Y are always between the minimum and maximum values on the boundary. In spite of this, the biharmonic method produces better meshes for most problems involving interface deformation as shown in the next example. Figure 11 shows results for a deforming interface with the two different interface boundary conditions. The initial mesh is the same as described in the previous problem. The final interface position is given by y = 3/8 cos(2πx), and the interface is deformed over 15 steps. The left half of Figure 11 shows the resulting mesh for the boundary conditions (X, Y ) and ∇(X, Y ) · ~n fixed separately on either side of the interface. (The result is reflected over the y-axis.) The right half shows the resulting mesh for (X, Y ) fixed and ∇(X, Y ) ·~n and ∇2 (X, Y ) continuous across the interface. Both meshes are fairly similar; Because of the stretching of the triangles near the interface, there is a significant compression of the triangles at the upper and lower boundaries. For this problem, it is not clear which boundary condition performs better, however we still favor the first boundary condition for the reasons discussed above. We also note that this calculation could not be completed with the Laplacian method because negative triangle areas occur. The problem of a moving cylinder in a channel was used as a test problem in [9] for an elastic deformation method. The initial mesh is shown in Figure 12. For the first calculation, we fix (X, Y ) on the boundary of the channel. We then move the cylinder one diameter to the right. Figure 13 shows the result for the Laplacian method. On the right half of the cylinder, the mesh lines are crossed and there are triangles with negative areas. Figure 14 shows the result for the biharmonic method. On the boundary of the channel ∇2 (X, Y ) = 0 and on the cylinder surface ∇(X, Y ) · ~n is fixed. In this case, the mesh remains well connected. Also note that the cells on the surface of the cylinder remain relatively undistorted. As a final test, we examine the problem of a moving cylinder with a refined boundary layer. Because of it’s relevance in fluid dynamics, similar problems have been examined by many previous authors [1, 4, 5, 13]. Figure 15 shows the initial mesh with an inset showing the refined region around the cylinder. We then move the cylinder 4.5 diameters to the left and upward and solve for the deformed mesh. Figure 16 shows the mesh after being deformed by the biharmonic method. The boundary conditions used are the same as for the previous problem of a cylinder in a channel. From the inset, we see that the mesh near the cylinder remains almost unchanged. The far-field mesh distorts, but none of the triangles develop high aspect ratios. With the Laplacian method, mesh lines cross and negative volumes occur so the problem cannot be completed. To fix this problem, several authors have proposed a modification to the Laplace method in which a variable diffusivity is introduced. Near the cylinder the diffusivity is increased making the mesh stiffer. There are several different implementations of this technique [1, 3, 4, 8]. In several of the studies it was found that increasing the diffusivity as the inverse of the square of the mesh length scale was the most effective for boundary layers [4, 8]. To compare the biharmonic method to this approach, we implement the “spring” method used in c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
MESH DEFORMATION USING THE BIHARMONIC OPERATOR
9
[3, 4]: Each side of the mesh is treated as a spring with a spring constant inversely proportional to the square of it’s length. Figure 17 shows the deformed mesh calculated using the spring method. For this problem, the spring method is as effective as the biharmonic method. The large variable diffusivity in the modified Laplacian/spring methods increases the stiffness of the equations to be solved compared to the standard Laplacian method. For the above problem however, the iteration convergence rate with multigrid and a Jacobi iteration is still faster than for the biharmonic equations. In addition, since the biharmonic equations are roughly twice as expensive as the modified Laplacian/spring approach, there is no reason to use the biharmonic for this problem. However for problems without a significant gradation in the mesh length scale (such as for most of the examples in this paper), the modified Laplacian/spring approaches will perform as poorly as the Laplacian method. Furthermore, for problems with boundary deformation in a localized region of near uniform mesh resolution (typical of an ALE two-phase flow computation), the biharmonic method should perform better.
Conclusions The ability to control both the position of boundaries and the normal mesh spacing is a significant advantage of the biharmonic method over second-order methods for deforming a mesh. This not only allows us to guarantee that mesh spacings will be continuous at interfacial or periodic boundaries, but it also allows the simulation of problems with greater boundary deformation than second-order methods without restructuring the mesh. This is valuable for both free-surface and interfacial flow calculations. Furthermore for problems with boundary layers, the biharmonic method can help preserve the spacing of the mesh in the boundary-layer so that restructuring is only needed in regions exterior to the boundary layer. This can also be accomplished by a variable diffusivity Laplacian method, which is computationally cheaper. However for problems of multiple bodies with interacting boundary layers or with localized boundary deformation, the biharmonic method may perform better. The main disadvantage of the biharmonic method is the cost. Using multigrid to solve the system of equations, the Laplacian method is about 1/4 the cost of the biharmonic method. Thus, the biharmonic method will not be optimal for all problems. However, for simulations in which the mesh movement is a small fraction of the the total cost, this may not be important. Furthermore, the fact that the need for mesh restructuring is reduced may also outweigh the increased cost.
ACKNOWLEDGEMENTS
Mesh generation for the example problems was done using EasyMesh which is free software written by Bojan Niceno. This work was supported in part by the Department of Energy under the ASCI/ASAP contract to Stanford University. Opinions expressed here are those of the authors and not of Stanford University or the Department of Energy.
REFERENCES c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
10
B. T. HELENBROOK
1. R. L¨ ohner and C. Yang. Improved ALE mesh velocities for moving boundaries. Comm. Num. Meth. Eng., 12:599–608, 1996. 2. I. Robertson and S. Sherwin. Free-surface flow simulation using hp/spectral elements. J. Comp. Phys., 155:26–53, 1999. 3. J. Batina. Unsteady euler algoithm with unstructured dynamic mesh for complex-aircraft aerodynamic analysis. AIAA Journal, 29(3):327–333, 1991. 4. V. Venkatakrishnan and D. J. Mavriplis. Implicit method for the computation of unsteady flows on unstructured grids. AIAA, 95-1705, 1995. 5. O. Hassan, E. J. Probert, and K. Morgan. Unstructured mesh procedures for the simulation of three dimensional transient compressible inviscid flows with moving boundary components. Int. J. Num. Meth. Fluids, 27:41–55, 1998. 6. C. Farhat, C. Degand, B. Koobus, and M. Lesoinee. Torsional springs for two-dimensional dynamic unstructured fluid meshes. Comp. Meth. Appl. Mech. Engrg., 163:231–245, 1998. 7. D. R. Lynch. Unified approach to simulation on deforming elements with application to phase change problems. J. Comput. Phys., 47:387–411, 1982. 8. T. E. Tezduyar, M. Behr, S. Mittal, and A. A. Johnson. Computation of unsteady incompressible flows with the stabilized finite element methods–space-time formulations, iterative strategies and massively parallel implementations. In P. Smolinksi, editor, New Methods in Transient Analysis, pages 7–24. ASME, New York, NY, 1992. 9. T. J. Baker. Mesh modification for solution adaptation and time evolving domains. In 7th International Conference on Numerical Grid Generation in Computational Field Simulations, 2000. 10. P. Lin, L. Martinelli, T. J. Baker, and A. Jameson. Two-dimensional implicit time dependent calculations for incompressible flows on adaptive unstructured meshes. AIAA, 01-2655, 2001. 11. D. J. Mavriplis. Multigrid techniques for unstructured meshes. Technical Report 95-27, ICASE, April 1995. 12. B. T. Helenbrook. A two-fluid spectral element method. Comp. Meth. Appl. Mech. Eng., 191:273–294, 2001. 13. T. E. Tezduyar, M. Behr, and J. Liou. A new strategy for finite element computations involving moving boundaries and interfaces—the deforming-spatial-domain / space-time procedure: Ii. computation of freesurface flows, two-liquid flows, and flows with drifting cylinders. Comp. Meth. App. Mech. Eng., 94:353– 371, 1992.
c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
MESH DEFORMATION USING THE BIHARMONIC OPERATOR
11
0.1 0 -0.1 -0.2 -0.3
y
-0.4 -0.5 -0.6 -0.7 -0.8 -0.9 -1
0
0.25
0.5
0.75
1
x Figure 1. Initial free-surface mesh
c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
12
B. T. HELENBROOK
0.1 0 -0.1 -0.2 -0.3
y
-0.4 -0.5 -0.6 -0.7 -0.8 -0.9 -1
0
0.25
0.5
0.75
1
x Figure 2. Free-surface mesh deformed with the Laplacian operator. Amplitude = 0.1
c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
MESH DEFORMATION USING THE BIHARMONIC OPERATOR
13
0.1 0 -0.1 -0.2 -0.3
y
-0.4 -0.5 -0.6 -0.7 -0.8 -0.9 -1
0
0.25
0.5
0.75
1
x Figure 3. Free-surface mesh deformed with the biharmonic operator. Amplitude = 0.1
c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
14
B. T. HELENBROOK
0.2 0.1 0 -0.1 -0.2
y
-0.3 -0.4 -0.5 -0.6 -0.7 -0.8 -0.9 -1
0
0.25
0.5
0.75
1
x Figure 4. Free-surface mesh deformed with the Laplacian operator. Amplitude = 0.25
c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
MESH DEFORMATION USING THE BIHARMONIC OPERATOR
15
0.2 0.1 0 -0.1 -0.2
y
-0.3 -0.4 -0.5 -0.6 -0.7 -0.8 -0.9 -1
0
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
1
x Figure 5. Free-surface mesh deformed with the biharmonic operator. Amplitude = 0.25
c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
16
B. T. HELENBROOK
0.2 0.1 0 -0.1 -0.2
y
-0.3 -0.4 -0.5 -0.6 -0.7 -0.8 -0.9 -1
0
0.25
0.5
0.75
1
x Figure 6. Free-surface mesh deformed with the Laplacian operator over ten intermediate steps.
c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
MESH DEFORMATION USING THE BIHARMONIC OPERATOR
17
0.2 0.1 0 -0.1 -0.2
y
-0.3 -0.4 -0.5 -0.6 -0.7 -0.8 -0.9 -1
0
0.25
0.5
0.75
1
x Figure 7. Free-surface mesh deformed with the biharmonic operator over ten intermediate steps.
c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
18
Maximum Correction to Vertex Position
B. T. HELENBROOK
10
-1
10
-2
10-3
Biharmonic 10-4 10
-5
10
-6
10-7
Laplacian 10-8 0
20
40
60
80
w-cycles Figure 8. Multigrid convergence rates for the Laplacian and Biharmonic methods.
c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
MESH DEFORMATION USING THE BIHARMONIC OPERATOR
19
1
y
0.75
0.5
0.25
0 0
0.25
0.5
0.75
1
x Figure 9. Laplacian method for deforming a mesh with a rising interface
c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
20
B. T. HELENBROOK
1
y
0.75
0.5
0.25
0 0
0.25
0.5
0.75
1
x Figure 10. Biharmonic method for deforming a mesh with a rising interface
c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
MESH DEFORMATION USING THE BIHARMONIC OPERATOR
21
1
0.75
y
0.5
0.25
0
-0.25
-0.5
0
0.5
x Figure 11. Interface deformation using the biharmonic method. Left: upper and lower meshes are treated separately. Right: mesh continuity enforced.
c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
22
B. T. HELENBROOK
y
1 0
-1
0
2
x 4
6
8
Figure 12. Undeformed mesh for cylinder movement in a channel
c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
23
MESH DEFORMATION USING THE BIHARMONIC OPERATOR
y
1 0
-1
0
2
x 4
6
8
Figure 13. Deformed mesh calculated with the Laplacian method for cylinder movement in a channel
c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
24
B. T. HELENBROOK
y
1 0
-1
0
2
x 4
6
8
Figure 14. Deformed mesh calculated with the biharmonic method for cylinder movement in a channel
c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
25
MESH DEFORMATION USING THE BIHARMONIC OPERATOR
y
10
0
-10
-10
0
x
10
20
Figure 15. Undeformed mesh for cylinder movement with a boundary layer. The inset shows the mesh around the cylinder
c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
26
B. T. HELENBROOK
y
10
0
-10
-10
0
x
10
20
Figure 16. Deformed mesh calculated using the biharmonic method for cylinder movement with a boundary layer. The inset shows the mesh around the cylinder
c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
27
MESH DEFORMATION USING THE BIHARMONIC OPERATOR
y
10
0
-10
-10
0
x
10
20
Figure 17. Deformed mesh calculated using the spring method for cylinder movement with a boundary layer. The inset shows the mesh around the cylinder
c 2001 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls
Int. J. Numer. Meth. Engng 2001; 00:1–30
