c
TU Delft, The Netherlands, 2006
ON THE SOLUTION OF THE DIE-SWELL PROBLEM
USING AN UNSTRUCTURED MESH
Jukka I. Toivanen∗ and Kari T. K¨arkk¨ainen†
∗University of Jyv¨askyl¨a, Department of Mathematical Information Technology P.O. Box 35 (Agora), FI-40014 University of Jyv¨askyl¨a, Finland
e-mail: jitoivan@mit.jyu.fi
†University of Jyv¨askyl¨a, Information Technology Research Institute e-mail: kari.karkkainen@titu.jyu.fi
Key words: Die-swell, Free Boundary Problems, Mesh Sensitivities
Abstract. In this work we consider the problem of determining the free surface for a stationary flow of viscous fluid. We introduce a solution algorithm based on interface-tracking of a geometry modeled by an unstructured mesh. We shall study the influence of the area of the geometry movement to the convergence of the Newton-type solution method. Detailed implementation for a model free boundary problem is presented along with the numerical results.
1 INTRODUCTION
Free boundaries are of crucial interest in many fluid flow applications. For example, an industrial process involving bubbles, jets or waves requires handling of a free boundary, which can play an important role in the quality or cost of a product. Therefore it is important to be able to simulate accurately free boundary problems of fluid flows.
The simulation of a free boundary problem requires the solution of state variables (e.g. pressure, velocity) in a geometry that is not known in advance. The solution geometry can be sought either by interface tracking or interface capturing methods11. In
the interface capturing method the solution geometry is constructed by using calculated fluid properties, e.g. density of the fluid. In the interface tracking method the geometry of the free fluid surface is solved by moving the nodal points presenting the free surface along with the free boundary.
In the following we shall sketch the proposed solution algorithm. First we calculate the dependence of each internal mesh nodal point to each nodal point at the free boundary. The dependencies are calculated efficiently using a Laplace type equation that controls the influence of each nodal point to the internal nodes. Second we solve the fluid flow problem iteratively in a sequence of geometries using a Newton type solution algorithm. The update to the geometry is solved at each Newton iteration step along with the changes to the state solution. The numerical results reveal that an accurate solution of the free fluid surface problem is achieved by 4–6 iteration steps.
An essential part in solving the free boundary problem is the calculation of the sensi-tivities of the residuals with respect to the geometry. Here we shall use two approaches. The first approach deals explicitly with the exact movement of the nodes of the mesh and the derivatives required for the Newton steps are calculated accurately. In the second ap-proach we calculate the sensitivities using a smaller area of influence of the free boundary movement, i.e. the contribution of the free boundary node movement is restricted in the neighborhood of the free boundary nodes but the geometry is transferred using a larger area of influence. Later we shall see that we can obtain a rapidly converging algorithm even though the derivatives are not calculated exactly.
In this paper we shall further consider the optimal solution strategies for free fluid surface problems. From the theory developed to shape optimization12 we know that
from the sensitivity calculation of the state equation the derivative with respect to the geometry is restricted on the boundary of the geometry. In references7,8 the automation
of the calculation of shape sensitivities is applied to solve the free boundary problems. In principle the numerical solution of a boundary value problem is sensitive also for internal nodal movement. This means that when the internal nodal points are moved along with the free boundary the state solution also changes slightly. If the internal nodal points are not taken into account in the calculation of the derivatives of the residuals for the Newton algorithm we can only expect linear convergence rate. However, the state solution boundary conditions can be chosen in a special way to guarantee that the state solution does not change when moving the free boundary in the normal direction. This was suggested by Garabedian3 for a model free boundary problem and later studied by
Cryer1 and Tiihonen14.
In the Newton algorithm the geometry of the free boundary is solved simultaneously with the state variables. FEM discretization of the coupled system yields a residual vector, which can be differentiated exactly for example using automatic differentiation4,5, and the
Newton iteration can be used to solve the coupled system efficiently.
In this approach, the sensitivities of the mesh node co-ordinates with respect to the boundary node locations are needed to differentiate the residual. Such sensitivities are not trivial to obtain if we use an unstructured mesh generator, since they often take iterative improvement steps13 which can make the relationship between the inner mesh node and
the boundary node locations highly nonlinear or even discontinuous.
pro-cedure that relocates the inner mesh nodes after the boundary displacement is known, and differentiate that procedure. For example the Laplace equation is widely used to propagate the boundary displacement into the domain, but with this approach, unless explicitly restricted, all the inner mesh nodes will have a non-zero sensitivity with respect to all the nodes on the free boundary. This significantly increases the time needed to solve the coupled system, since certain blocks of the Jacobian of the residual become dense.
In the reference10 the Poisson equation is used to estimate the distance from the
bound-ary. Similarly, we use the anisotropic Laplace equation to find a neighborhood of each boundary node. The displacement of the boundary node is then propagated into its neigh-borhood using the same equation, this time setting a homogeneous Dirichlet condition for the nodes that are ’far’ and eliminating them from the system. This results in a different subproblem for each node on the free boundary, but the subproblems can be solved effi-ciently since they are small and the related matrices representing the Laplace operator are sparse and positive definite. Moreover, it is sufficient to assemble the matrix representing the Laplace equation once.
2 Model free boundary problem
In this paper we shall apply the suggested solution algorithm to a classical benchmark-ing free boundary problem, so called 2D die-swell problem. It models the extrusion of a viscous incompressible jet from a die. Unknowns are the fluid velocity ~u, the hydrostatic pressure p and the free jet surface.
Problem geometry is presented in figure 1. The flow is assumed to be symmetric with respect to the x1-axis, which is denoted with Γs in the figure. Other boundary sections
are the free surface Γf, the wall of the die Γw, the inlet Γi and the outlet Γo. Dimensions
used in this example were b = 2, L1 = 3.5 and L2 = 20.0.
b/2 i Γ Γw Γf Γo Γs x1= L 2 x1= − L 1 x1= 0 P
The fluid flow is governed by the Navier-Stokes equations −∇ · σ + ρ~u · ∇~u = 0 in Ω
∇ · ~u = 0 in Ω (1)
where
σ = µ(∇~u + ∇~uT) − pI.
Boundary conditions are such that
u1 = 32U (1 − 4x22/b2), u2 = 0 on Γi ~u = 0 on Γw ~n · σ · ~n = 0, ~u · ~t = 0 on Γo ~t · σ · ~n = 0, ~u · ~n = 0 on Γs ~n · σ · ~n = λκ, ~t · σ · ~n = 0, ~u · ~n = 0 on Γf (2)
where λ is the surface tension coefficient, κ is the curvature of the free boundary, U is the mean velocity of the flow at the inlet, b is the inlet height and ~n and ~t are the surface normal and tangent vectors.
The equations were discretized using the finite element method. Quadratic triangular elements were used for the velocity and linear triangular elements for the pressure. The variational formulation reads9
r1 = Z Ω σ1· ∇ϕ dx + Z Ω ρu · ∇u1ϕ dx + λ Z Γf t1 dϕ ds − λϕ(P ) = 0 (3) r2 = Z Ω σ2· ∇ϕ dx + Z Ω ρ~u · ∇u2ϕ dx + λ Z Γf t2 dϕ ds − λϕ(P ) = 0 (4) r3 = Z Ω (∇ · ~u) φ dx = 0 (5) r4 = Z Γf (~u · ~n) φ dx = 0 (6)
where ϕ and φ are the quadratic and linear basis functions respectively and σi is the i:th
row of the tensor. The problem is characterized by the Reynolds number Re = bU ρµ and inverse of the capillary number 1/Ca = λ
µU.
Similarly to6, we solved the fluid flow and the location of the free boundary
simulta-neously. The coupled system had thus four-by-four blocks, where the unknowns were the velocity components q1 and q2, the pressure q3 and the x2 co-ordinates of the nodes on
the free boundary q4.
We used the Newton linearization to solve the coupled system. That is, in each iteration the solution q is updated by setting q := q − δq, where δq is solved from
∂r ∂q
Here r is the residual vector and ∂r/∂q is its Jacobian.
Part of the Jacobian now contains the sensitivities of the residual with respect to the x2 locations of the nodes on the free boundary, which according to the chain rule can be
written as ∂r ∂q4 = ∂r ∂X ∂X ∂q4 . (8)
Since the matrix ∂r/∂X is constant, the sparsity of the Jacobian blocks ∂r/∂q4 is
deter-mined solely by the structure of the nodal sensitivity field ∂X/∂q4.
3 UTILIZING UNSTRUCTURED MESH
In general, unstructured meshes have some significant advantages over structured meshes. For example, the mesh generation phase is easier to automate, and the mesh can be refined locally near the difficult locations regarding the flow solution without significantly increasing the number of degrees of freedom. In figure 2 we can see an un-structured mesh used in the solution of the die-swell problem. The mesh is refined near the exit of the die because the solution of the state equation has a singularity at that location. In figure 4 there is a structured mesh with a similar number of nodes and re-finement at the same location. The element size at that area is nonetheless much larger than in the unstructured mesh.
With unstructured meshes, the mesh sensitivities ∂X/∂q4 needed in (8) are not trivial
to obtain, since unstructured mesh generators often take iterative improvement steps13,
which can make the relationship between the inner mesh node and the boundary node locations non-smooth or even discontinuous.
Instead of differentiating the mesh generator, we can use some mesh deformation pro-cedure that relocates the inner mesh nodes after the boundary displacement is known, and differentiate that procedure. For example, we can assume that the x1 co-ordinates
of the nodes are fixed, and the mesh node displacements v in x2 direction are given as a
solution to the anisotropic diffusion equation c1 ∂2v ∂x2 1 + c2 ∂2v ∂x2 2 = 0 (9)
with Dirichlet boundary conditions. In this work, we used the values c1 = 0.1 and c2 = 1.0
for the diffusion coefficients.
Being a linear equation, (9) defines a linear mapping L from the boundary node lo-cations q4 to the inner mesh node locations X. The mesh node locations are thus given
by
X = X0+ L(qn
4 − q04) (10)
where X0 are the initial co-ordinates and qi
4 are the locations of the boundary nodes at
step i. The mesh sensitivities are then ∂X ∂q4
In this approach, unless explicitly restricted, all the inner mesh nodes will have a non-zero sensitivity with respect to all the nodes on the free boundary. This significantly increases the time needed to solve the linearized equations (7), since the residual blocks ∂r/∂q4 become dense.
In10 the Poisson equation is used to estimate the distance from the boundary. We can
do the same thing using the equation (9). We assemble the matrix A representing the discretized equation (9), and for each node i on the free boundary we solve the system
Adi = ei
where ei is the i:th Euclidean basis vector.
As the solution to these equations we get the vectors di, which can be used as distance
estimates from the corresponding free boundary nodes. We then define the index sets Ii =j|dji < tol , which represent the nodes in the neighborhood of each boundary node.
Here tol is a tolerance parameter. In figure 2 we can see the neighborhood of one boundary node estimated in this way using the tolerance parameter tol = 0.05.
In the next step, we solve the condensed systems Aivi = ei,
where only those rows and columns of A and entries of ei that are indexed by Ii are
present.
For later use, we define two linear mappings from the nodes on the free boundary to the whole domain. Columns of the first mapping Ld are composed of the vectors di, and
the matrix Ld is thus dense. The other is a sparse mapping Ls, whose sparsity pattern is
given by the sets Ii and the non-zero entries are given by the vectors vi.
Figure 2: Part of an unstructured mesh and the estimated neighborhood of one boundary node (tol = 0.05).
The previous computations are performed using linear elements. Since we always require that the extra nodes in the quadratic elements are located in the midpoint of two element vertices, those entries in the matrices Ld and Ls that correspond to the
4 NUMERICAL RESULTS
The amount of non-zero entries in the sensitivity field ∂X/∂q4 effects the time needed
to perform one solution iteration in two ways. The Jacobian ∂r/∂q was computed using the so called forward mode of automatic differentiation, and the computational cost is therefore directly proportional to the number of partial derivatives that has to be com-puted. The structure of the sensitivity field also has a direct effect on the sparsity of the Jacobian, which affects the time needed to solve the linear system (7). The package SuperLU2 was used as the linear solver. It implements LU decomposition and reordering
of the unknowns for sparsity.
4.1 Comparison of different approaches
We tested different approaches to solve the die-swell problem. The results are illus-trated in figure 3, where the convergence of the residual versus the computation time using different methods is plotted. In this test we used an unstructured mesh with 6297 nodes. The non-dimensional parameters had the values Re = 75 and 1/Ca = 2.0.
For the previously mentioned reasons, we should keep the sensitivity field ∂X/∂q4 as
sparse as possible to minimize the computation time. However, the mapping used to deform the mesh according to the equation (10) can not be too sparse, or the mesh will become distorted. Therefore, we should basically use as large a tolerance parameter tol as possible such that the mapping Ls is still feasible in terms of mesh deformation. In
our tests the tolerance parameter 0.05 proved to be sufficient, but the required tolerance obviously depends on many things, such as the maximum deformation and the density of the mesh. This approach is called case A in figure 3.
Due to the simplicity of the geometry, we can form the mapping Ls analytically in the
following way. If for node k with co-ordinates (xk
1, xk2) there exists two consecutive free
boundary nodes a and b with x1 co-ordinates xa1 and xb1 such that xa1 ≤ xk1 ≤ xb1 we set
Ls k,a = xk2 xb 1− xk1 xb 1− xa1 and Ls k,b= xk2 xk 1 − xa1 xb 1− xa1 .
Computation with this kind of mapping is called case B in the figure. In case of more complex geometries finding such relationships may not be possible.
Another possibility is to only compute gradients approximately. That is, we deform the mesh with the dense mapping Ld, but to compute the Jacobian we set ∂X/∂q
4 = Ls.
Now the tolerance parameter can be kept relatively large. A run where the mapping Ls
was computed using a tolerance parameter 0.2 is plotted as case C in the figure. The drawback of this kind of approach is that since the gradients are only approximate, a larger number of iterations is needed for convergence.
kind of field is called case D in the figure 3. We can see, that this approach is very competitive for this problem. Although we need slightly more iterations than with the previous methods, the convergence is still quite fast. Also, the computation time per iteration is very small.
As a reference we also tested a field where only the nodes on the boundary have a non-zero sensitivity. This is called case E in the figure. A clear difference to case D can be seen there. If we restrict the field strictly on the boundary, the convergence rate becomes quite poor and the total computation time increases.
The obtained results vary a little depending on the mesh used, but the overall perfor-mance of the methods is similar. It can clearly be seen, that the cases where the Jacobian is computed exactly (cases A and B) converge in a fewer number of iterations than the approximate methods. However, in terms of the total computation time the difference is not that clear anymore.
−14 −12 −10 −8 −6 −4 −2 0 2 0 20 40 60 80 100 120 CASE A CASE B CASE C CASE D CASE E
Figure 3: Convergence of the residual versus the computation time using different methods.
4.2 Comparison of structured and unstructured meshes
The traditional approach to solve the die-swell problem is to use a structured mesh, for which the nodal sensitivities are easy to obtain and the sensitivity field is inherently sparse. An example of a structured mesh and the nodes that are affected by the movement of one boundary node are illustrated in figure 4.
Figure 4: Part of a structured mesh and the nodes affected by the movement of one boundary node. equation has a singularity. We performed a series of computations on different structured and unstructured meshes. The number of nodes varied from 511 to 23805 on structured and from 1161 to 17657 on unstructured meshes.
For the structured mesh the relationship between the domain and the boundary nodes, and thus also the nodal sensitivities, are trivial to obtain. With unstructured meshes we used the previously described approach, where the mapping Ls was computed using
the anisotropic Laplace equation (tolerance parameter 0.05) and the Jacobian matrix was computed exactly.
The non-dimensional parameters had the values Re = 75 and 1/Ca = 2.0. We were interested in the so called die-swell percentage defined as 100(h(Γo) − h(Γi))/h(Γi), where
h(Γi) and h(Γo) are the heights of the inlet and the outlet respectively. The obtained
die-swell percentage versus the total computation time, including the formation of the mapping Ls, on different meshed is plotted in figure 5. Our computations with dense
structured and unstructured meshes suggest the value −11.13 for the die-swell percentage. Other values appearing in the literature for these parameters are for example -11.21 in reference9 and -11.06 in reference6. We can conclude that realistic results can be obtained
with low computational cost using an unstructured mesh with refinement near the exit of the die. -11.5 -11.4 -11.3 -11.2 -11.1 0 50 100 150 200 250 Structured Unstructured
Part of the final deformed mesh is shown in figure 6 along with the pressure contours of the solution. We can see, that the quality of the mesh is still relatively good.
Figure 6: The deformed mesh and the pressure contours near the exit of the die.
5 CONCLUSIONS
In this paper we studied an implementation of a quasi-Newton type solution algorithm to solve a stationary viscous incompressible free fluid surface problem. We presented a generic method to extend the movement of the free boundary to the interior nodal points by solving an anisotropic Laplace equation. This enabled an efficient and generic calculation of shape sensitivities using automatic differentiation. The algorithm was tested for a simple geometry but the intended use of the method is for complex geometries. Numerical results indicated that the method is competitive against analytically extended linear geometry transformation and thus the suggested approach gives the possibility to study more complex free boundary problems modeled on an unstructured mesh easily and more accurately.
REFERENCES
[1] C. W. Cryer. On the approximate solution of free boundary problems using finite differences. Journal of the Association for Computing Machinery, 17(3):397–411, July 1970.
[3] P.R. Garabedian. The mathematical theory of three dimensional cavities and jets. Bulletin of the American Mathematical Society, 62:219–235, 1956.
[4] A. Griewank. Evaluating Derivatives: Principles and Techniques of Algorithmic Dif-ferentiation. SIAM, 2000.
[5] J. Haslinger and R. A. E. M¨akinen. Introduction to Shape Optimization: Theory, Approximation, and Computation. Society for Industrial and Applied Mathematics, Philadelphia, 2003.
[6] K. Hiltunen, M. Laitinen, A. Niemist¨o, and P. Tarvainen. Using mathematical con-cepts in software design of computational mechanics. In Proceedings of the VIII Finnish Mechanics Days, 2003.
[7] K. K¨arkk¨ainen and T. Tiihonen. Free surfaces: shape sensitivity analysis and numerical methods. International Journal for Numerical Methods in Engineering, 44(8):1079–1098, 1999.
[8] Kari K¨arkk¨ainen. Shape Sensitivity Analysis for Numerical Solution of Free Boundary Problems. PhD thesis, University of Jyv¨askyl¨a, 2005.
[9] N.P. Kruyt, C. Cuvelier, A. Segal, and J. van der Zanden. A total linearization method for solving viscous free boundary flow problems by the finite element method. International Journal for Numerical Methods in Fluids, 8:351–363, 1988.
[10] R. L¨ohner and C. Yang. Improved ale mesh velocities for moving bodies. Communi-cations in Numerical Methods in Engineering, 12:599–608, 1996.
[11] Anton Smolianski. Numerical modeling of two-fluid interfacial flows. PhD thesis, University of Jyv¨askyl¨a, 2001.
[12] Jan Sokolowski and J.P. Zol´esio. Introduction to shape optimization. Springer-Verlag, 1992.
[13] Joe F. Thompson, Bharat K. Soni, and Nigel P. Weatherhill, editors. Handbook of Grid Generation. CRC Press, 1999.
