- Email: [email protected]
Finite element mesh update methods for fluid–structure interaction simulations
Available online at www.sciencedirect.com
Finite Elements in Analysis and Design 40 (2004) 1259 – 1269 www.elsevier.com/locate/!nel
Finite element mesh update methods for &uid–structure interaction simulations Zhenlong Xu, Michael Accorsi∗ Department of Civil and Environmental Engineering, University of Connecticut, Storrs, CT 06269-2037, USA Received 14 November 2002; received in revised form 24 April 2003; accepted 19 May 2003
Abstract Mesh updating is a problem commonly encountered in &uid–structure interaction simulations when structures undergo a large displacement. In this paper, a set of mesh updating strategies based on a “pseudo-solid” model is proposed. Generally, a homogeneous analysis is !rst performed on the pseudo-solid with prescribed boundary displacements. Di5erent element properties and external forces are then introduced based on strain results from the !rst analysis. A second analysis is performed on this non-homogeneous model with the goal to maintain element aspect ratio while preserving element volumes. Several numerical simulations are conducted using di5erent combinations of strategies proposed in this paper. In general, the proposed strategies signi!cantly enhance the performance of the pseudo-solid mesh updating method. ? 2003 Elsevier B.V. All rights reserved. Keywords: Fluid–structure interaction; Mesh update; Finite element method
1. Introduction In &uid–structure interaction (FSI) simulations, it is quite common for the structure to undergo large displacements. To maintain compatibility between the structural and &uid meshes, re-meshing is usually performed on the &uid mesh to accommodate this change. However, for large-scale simulations, this process is computationally expensive. An alternative approach to achieve this goal is mesh moving, in which the current &uid mesh is moved along with the deforming structure so as to preserve the compatibility between these two meshes. The most common approach for mesh moving is the “pseudo-solid” method. In this method, the &uid mesh is viewed as a linear elastic pseudo-solid, thereby turning the mesh-moving problem ∗
Corresponding author. Tel.: +1-860-486-5642; fax: +1-860-486-2298. E-mail address: [email protected] (M. Accorsi).
0168-874X/$ - see front matter ? 2003 Elsevier B.V. All rights reserved. doi:10.1016/j.!nel.2003.05.001
1260
Z. Xu, M. Accorsi / Finite Elements in Analysis and Design 40 (2004) 1259 – 1269
into a linear solid mechanics problem, with the displacement of the actual structural mesh acting as boundary conditions. Because the sole purpose of the pseudo-solid model is to move the &uid mesh, the solution to this problem itself has no physical importance. Therefore, di5erent methods can be employed to manipulate the properties of the pseudo-solid so as to achieve a mesh with minimum element distortion and volume change. Strategies based on the pseudo-solid method generally fall into two categories, one-step and two-step methods. In one-step methods, sti5ness is assigned to di5erent elements according to their initial geometric properties. Typical criteria include element sizes [1–3] or their distance to a boundary of prescribed displacements. Two-step methods are based on a predictor–corrector philosophy. In the !rst step (predictor), a linear structural analysis is conducted on the pseudo-solid mesh based on a homogeneous material model subjected to the prescribed boundary displacements. Element strains obtained from this step are then utilized as the basis for assigning sti5ness to di5erent elements. A second analysis (corrector) is then conducted on this non-homogeneous solid model and results of this analysis will be used as the desired motion of the &uid mesh. The success of two-step methods depends on the selection of element sti5ness based on the !rst analysis. A study examining various two-step methods has been presented by Chiandussi et al. [4]. In this paper, several new strategies for two-step methods are proposed based upon identi!cation of key reasons for element distortion. These methods can be either implemented independently, or combined with other approaches. Numerical simulations are then conducted on a 2-D mesh to compare the e5ectiveness of the di5erent strategies. Generally, the new methods proposed in this research are shown to be highly e5ective in preventing extremely distorted elements while maintaining the overall quality of the &uid mesh. The size of the required mesh motion corresponds to the motion of the structure during the FSI simulation, which is a function of many parameters. The goal of a mesh moving procedure is to control element distortion and maintain a high-quality mesh for as large a motion as possible. For large motions, it is anticipated that mesh moving will become less e5ective and re-meshing will be required. Therefore, in the numerical examples, the e5ectiveness of the mesh moving procedures is evaluated over a large range of mesh motion sizes. 2. Methodology 2.1. De1nition of mesh quality Since the ultimate goal of mesh moving is to achieve a high-quality !nite element mesh subjected to prescribed structural motions, a general measure of mesh quality should be de!ned. The most important factor that de!nes mesh quality is the aspect ratio of an element. It is generally acknowledged that elements with large aspect ratios perform poorly in !nite elements analyses. Therefore, a successful mesh moving strategy should not only maintain the overall quality of the mesh, but also be able to limit, if not avoid totally, the number of elements with extremely large aspect ratios. In this paper, for a typical 3-D tetrahedral element, we adopt the r=R value as the indicator of the element aspect ratio, in which r and R are the radii of the inner- and circum-sphere of the element, respectively. For 2-D analysis, r and R are the radii of the inner- and circum-circle of a triangular element [5].
Z. Xu, M. Accorsi / Finite Elements in Analysis and Design 40 (2004) 1259 – 1269
1261
2.2. Pseudo-solid formulation The basic !nite element formulation for the pseudo-solid system is K1 d1 = 0 d 1 = d0
(); ( = s )
(1)
in which K1 is the sti5ness matrix, d1 the mesh motion to be determined, and d0 the boundary displacements prescribed by the motion of the actual structural mesh. If this analysis serves as the !rst step in a two-step method, K1 is usually based on a homogeneous material model for the entire mesh. If this formulation is for a one-step method, K1 is generally non-homogeneous in the sense that di5erent elements may have di5erent material models (i.e. variable sti5ness) to control their deformation. This will be further explained in the next few sections. For a typical two-step method, a second analysis is also performed on the pseudo-solid mesh with designated element sti5ness, which is assigned to each element according to their strain obtained from the !rst step. The formulation for this step is K2 (d1 )d2 = F2 (d1 ) d 2 = d0
( = s )
(); (2)
in which d2 is the motion to be determined. The sti5ness matrix is written as K2 (d1 ) to re&ect its dependence on the !rst analysis. Unlike the homogeneous analysis, external forces could also be introduced into the second analysis to eliminate local deterioration of mesh quality so long as the boundary condition is satis!ed. Since the external force is also related to results from the !rst analysis, it is written as F2 (d1 ) in the above equation to re&ect this in&uence. Several strategies can be utilized in selecting K2 (d1 ) and F2 (d1 ). They will be explained in the following sections. 2.2.1. Poisson’s ratio For isotropic solids, the value of Poisson’s ratio re&ects the in&uence of the motion in one direction on the motion in the other two perpendicular directions. The range for its value is −1 ¡ ¡ 0:5 based on the requirement that the material has a positive bulk modulus and a positive shear modulus. Poisson’s ratio values for typical solids are positive. For these materials, compression in one direction will produce elongation in the other two directions, and vice versa. However, this kind of behavior has a negative e5ect on the geometry of elements in that elongation in one direction and shortening in the other two leads to degradation of the element aspect ratios. As mentioned previously, results obtained from the pseudo-solid model have no physical significance. Therefore, material properties can be chosen solely for the purpose of conserving element geometry without concern about whether such materials exist in the physical world or not. For this reason, use of a negative Poisson’s ratio is proposed in this research. The e5ect of this adoption is apparent. For example, with such a property, an element subjected to compression in one direction would tend to shorten in the other two perpendicular directions, thereby preserving the initial aspect ratio. However, this approach does have certain disadvantages in preserving element volume. Using the same example, it is apparent that a negative Poisson’s ratio may cause elements to contract too much under compression or elongate too much under tension.
1262
Z. Xu, M. Accorsi / Finite Elements in Analysis and Design 40 (2004) 1259 – 1269
Excessive change in element volume can be prevented by combining this method with other element sti5ening methods, as described in the following sections. 2.2.2. Element sti6ening Non-homogeneous element sti5ening is the primary idea underlining almost all mesh-moving methods based on pseudo-solid models. The di5erence between di5erent methods lies only in the selection of sti5ening functions, which usually relate to element geometry for one-step methods and element deformation for two-step methods. A typical one-step and two-step sti5ening scheme will be brie&y reviewed along with a new sti5ening scheme proposed in this research. The two former sti5ening schemes are given here primarily for the purpose of comparing results. The element sti5ness matrix for the pseudo-solid elements can be expressed as Ke = BT DB |J |e e d; (3)
where B is the derivative matrix of shape functions, D the constitutive matrix, |J |e the Jacobian of the element, and e the element sti5ening factor for mesh moving. Masud’s one-step sti6ening scheme: As a typical method in the one-step category, Masud [1] introduced a sti5ening scheme based on element sizes, i.e. choosing e as |J |max − |J |min ; (4) e = 1 + |J |e where |J |max and |J |min are the maximum and minimum Jacobian among all elements, respectively. This method essentially increases the sti5ness of small elements to prevent them from being distorted by large motions. Chiandussi’s two-step sti6ening scheme: Chiandussi et al. [4] proposed several two-step sti5ening schemes based on di5erent measures of element deformation, which were shown to have similar performance. Due to this reason, only one of those methods is presented here. It is a sti5ening scheme based on the square norm of element strains, i.e. 1 12 + 22 + 32 (5) N 2 3 in which 1 ; 2 and 3 are the three principal strains for an element based on results calculated from the !rst homogeneous analysis, and N is a user speci!ed strain constant for that element. Since, in practice, the same N is typically used for all the elements, this value has no e5ect on the !nal results for a linear analysis. Two-step sti6ening scheme based on element distortion: In this section, a sti5ening method is introduced based on the observations of the causes of element distortion. The key factors that a5ect element geometry are the maximum and minimum principal strains. The reason for this lies in the notion that the quality of an element’s geometry (i.e. aspect ratio) is essentially decided by its largest and smallest dimensions. Therefore, for an element with a near-optimal aspect ratio, which should be the case for most of the elements in a well-constructed initial mesh, the deformed shape should be controlled by the two extreme principal strains. While the in&uence of the middle principal strain on element geometry exists, its e5ect is not as essential as the other two. This criterion is noticeably di5erent than those introduced by Chiandussi et al., where all three principal strains are considered to have the same e5ect on element geometry. e =
Z. Xu, M. Accorsi / Finite Elements in Analysis and Design 40 (2004) 1259 – 1269
1263
Based on this notion, the element-sti5ening factors are chosen to be functions of (1 − 3 ), which can also be interpreted as the diameter of the largest Mohr’s circle for a 3-D strain state. Obviously, larger values of (1 − 3 ) indicate more severe distortion of elements, which, in turn, require higher sti5ening factors. In practice, many di5erent functions could be selected, ranging from linear, quadratic to exponential, all with (1 −3 ) acting as the only variable. For this research, the following exponential formulation is chosen: e = a(1 −3 ) ;
(6)
where a is a constant to be decided. In9uence of tension and compression: A very important but largely neglected factor that a5ects element geometry is the di5erent e5ects of tension and compression. In this research, it is recognized that compression generally tends to have a more severe e5ect on the degradation of element quality than tension. This can be demonstrated through a very simple example. For a 2-D element with an optimal initial aspect ratio of , assume the two principal strains are 0 and 0 after a certain deformation. If 0 is tensile (i.e. positive), the new aspect ratio for this element will be roughly (1 + 0 ). If 0 is compressive (i.e. negative), then the new aspect ratio will be approximately =(1 + 0 ). For a mesh-moving problem, it is not uncommon to have very large strains caused by large motions of the structural mesh. Given a moderately large strain value of ±0:5, the !nal aspect ratio would be 1:5 for the tension case (0 = 0:5) and 2 for the compression case (0 = −0:5). Clearly, compression has caused signi!cantly more damage to the aspect ratio than tension. This e5ect would be even more pronounced for extremely large strains. Therefore, any sti5ening technique that neglects this di5erence would be discriminating against elements under compression and potentially lead to a mesh with more distortion. In this research, the in&uence of tension/compression is accounted for by introducing a modifying factor 1=(1 + 3 ) to the element sti5ener e , which now becomes a(1 −3 ) e = : (7) 1 + 3 Deviatoric stress removal: As previously expressed in the basic formulation, external forces can also be introduced to eliminate element distortion. Generally, this is not implemented in the !rst step due to the uncertainty in guessing external forces that are bene!cial. However, this technique can be utilized in the second step based on results obtained from the !rst step. It is well known that a stress state can be decomposed into two parts: the hydrostatic and deviatoric stress. While the hydrostatic stress is responsible for a change in the element volume, the deviatoric stress is related to the element distortion. For a linear element with no displacement constraints, the release of the deviatoric stress results in a purely hydrostatic stress state with no element distortion. The consistent element force needed to remove the hydrostatic stress is F2 (d1 ) = − BT D (d1 )|J | d; (8)
where D (d1 ) is the deviatoric stress based on the results from the !rst analysis. As previously identi!ed by Stein [6], the existence of boundary constraints prevents elements from achieving a distortion-free state. This is also true for the proposed method of stress removal. In the second part of the two-step method, a complete removal of deviatoric stresses cannot be achieved
1264
Z. Xu, M. Accorsi / Finite Elements in Analysis and Design 40 (2004) 1259 – 1269
due to this in&uence. However, our assumption is that even a partial removal of these stresses will be bene!cial in eliminating local element distortion. It should also be noted that this approach is based solely on linear mechanics, whereas the deformation of the pseudo-solid, in general, can be quite large and possibly nonlinear. This is expected to also limit the removal of deviatoric stress in the elements. 3. Numerical results The strategies discussed above have been applied to a 2-D example, which is the same as used by Stein [6]. In this example, the mesh occupies the region de!ned by |x| 6 1:0 and |y| 6 1:0, as shown in Fig. 1. An interior surface is located in the middle of the mesh, speci!ed as |x| 6 0:5 and y = 0:0. In general, the strategies can be applied to 3-D models without any additional diQculties. Three types of motion of the inner surface were studied in this paper: (1) a uniform vertical translation in the y direction, (2) a rotation with respect to the center of the mesh, and (3) a deformation of the inner surface from a line to a circular arc with no stretching. For each motion, the three di5erent sizes given in Table 1 were used to study the relative e5ectiveness of the methods
Fig. 1. (a) Initial mesh, (b) initial mesh (enlarged).
Table 1 Motion types and sizes Relative scale Motion I (translation) Motion II (rotation) Motion III (deformation)
0.1 0.1 =20 =5
0.2 0.2 =10 2=5
0.5 0.5 =4
Z. Xu, M. Accorsi / Finite Elements in Analysis and Design 40 (2004) 1259 – 1269
1265
Table 2 Six strategies for mesh moving Strategy Strategy Strategy Strategy Strategy Strategy
1 2 3 4 5 6
Masud’s method (Eq. (4)) Chiandussi et al. two-step sti5ener (Eq. (5)) Two-step sti5ener with tension/compression in&uence (Eq. (7)) Strategy 3 + deviatoric stress removal (Eq. (8)) Strategy 4 + negative Poisson’s ratio of −0:5 Masud’s method with two half-steps
in handling small, medium, and large motions. For uniformity, all motion sizes are reported on a relative scale ranging from 0 to 0.5. The strategies used in this study are listed in Table 2. Considering the fact that a two-step method is twice as expensive, computationally, as a one-step method for the same motion size, strategy 6 is included for comparison where Masud’s one-step method is used with two half-steps rather than just one. Except for strategy 5, the Poisson’s ratio used in all other strategies has a value of 0.3. Young’s modulus of the pseudo-solid can be chosen arbitrarily in all cases. Also, a = 2 is used for the parameter in Eq. (7). For strategies 3–5 a homogeneous solid is used for the !rst step, and the modi!cations listed are applied in the second step. As described in the methodology, the ratio between the radii of the circum- and inner-circle is used as the criterion for mesh quality. For each motion, the average ratio for the entire mesh is used to evaluate the overall quality of the mesh, and the worst case (smallest ratio of all elements) is also presented to evaluate the worst degradation of mesh quality. For the initial mesh, the average ratio is 0.4401 and the worst ratio is 0.2629. The computer program used for this study is based on matrix-free techniques for formulating the system equations and an iterative solver GMRES [7] was adopted to solve those equations. All the results calculated in this research were obtained by using 100 iterations in the inner iterative loop of GMRES. 3.1. Motion I—vertical translation Fig. 2 shows a successful mesh motion for a prescribed vertical displacement of 0.5. In Fig. 3, the average and worst ratios for di5erent mesh movement sizes using di5erent methods are presented. Obviously, the overall mesh quality deteriorates as the mesh motion becomes larger. Although the di5erent methods appear to have similar performance in this aspect, the di5erence between them tends to grow, as the prescribed displacements get larger. Because the displacement increment used in strategy 6 is actually only half the size used in the other strategies, strategy 6 is able to achieve a better overall quality than all other methods. Under large displacements, the worst ratio value drops dramatically, from around 0.23 to about 0.1. The three strategies proposed in this research (i.e. 3–5) perform noticeably better than other methods for large displacement. With the help of deviatoric stress removal, strategy 4 is slightly better than strategy 3, and strategy 5 outperforms strategy 4 by additionally utilizing a negative Poisson’s ratio. In practice, controlling the worst case is most critical because failure of one element (i.e. element tangling) is the limiting factor for the mesh-moving process.
1266
Z. Xu, M. Accorsi / Finite Elements in Analysis and Design 40 (2004) 1259 – 1269
Fig. 2. (a) Motion I—vertical translation, (b) motion I—vertical translation (enlarged).
Worst r/R Ratio -- Motion I
Average r/R Ratio -- Motion I
0.3
0.45 Strategy 1 Strategy 2 Strategy 3 Strategy 4 Strategy 5 Strategy 6
0.43 0.41 0.39 0.37
0.2 0.15 0.1 0.05 0
0.35 0
(a)
Strategy 1 Strategy 2 Strategy 3 Strategy 4 Strategy 5 Strategy 6
0.25
0.2
0.4
0.6
0
0.2
0.4
0.6
(b)
Fig. 3. (a) Mesh quality for motion I—average r=R ratio, (b) mesh quality for motion I—worst r=R ratio.
A very important property for this translational motion is that it roughly divides the pseudo-solid mesh into two regions: the compression and tension part, with a small fraction of elements lying between these two areas. From Fig. 2, it is quite obvious that elements in the compression region deform much more severely than those in the tension region. This observation validates the reasoning for using a sti5ening scheme that accounts for the e5ect of tension and compression. In Fig. 3, strategy 2 has surprisingly low values for the worst-case ratio. This can be explained by a more detailed look at the deformed shape, which is not shown here. The worst-case element is not located near the tips of the structural surface, as previously assumed, but rather on the border of the mesh somewhere between the tension and compression region. A plausible explanation is that the sti5ening formulation in this strategy did not adequately protect elements with negative
Z. Xu, M. Accorsi / Finite Elements in Analysis and Design 40 (2004) 1259 – 1269
1267
Fig. 4. (a) Motion II—rotation, (b) motion II—rotation (enlarged).
Worst r/R Ratio -- Motion II
Average r/R Ratio -- Motion II
0.3
0.45 Strategy 1 Strategy 2 Strategy 3 Strategy 4 Strategy 5 Strategy 6
0.43 0.41 0.39 0.37
0.25
Strategy 1 Strategy 2 Strategy 3 Strategy 4 Strategy 5 Strategy 6
0.2 0.15 0.1 0.05
0.35 0
0.2
(a)
0.4
0.6
0 0
0.2
0.4
0.6
(b)
Fig. 5. (a) Mesh quality for motion II—average r=R ratio, (b) mesh quality for motion II—worst r=R ratio.
principal strains since it does not consider the stronger in&uence on element geometry caused by compression. 3.2. Motion II—rotation Fig. 4 shows a successful mesh motion for a prescribed rotation of =4. In Fig. 5, the average ratio and worst ratio for di5erent mesh movement sizes using di5erent methods are presented. The trend in average and worst values agrees with results listed in the previous section, i.e. mesh quality degrades with the increase in magnitude of the prescribed displacement. One noticeable di5erence between the deformed mesh for this motion and that under translation is that most elements in this case did not undergo a large change in volume, but rather a change in shape. This is apparent by comparing Figs. 2 and 4.
1268
Z. Xu, M. Accorsi / Finite Elements in Analysis and Design 40 (2004) 1259 – 1269
Fig. 6. (a) Motion III—deformation, (b) motion III—deformation (enlarged).
For this motion, strategy 6 outperforms the other methods on both measures. One possible explanation for this is that this strategy is based on the idea of preserving element volume, which happens to agree with the behavior of the deformed mesh under this motion. The methods proposed in this research performed fairly well, although not as good as strategy 6. Strategy 2 once again was shown to have the lowest value in worst ratio. 3.3. Motion III—deformation This case is the most demanding one among all the three motions considered, and strategies 1 and 2 actually fail for the largest displacement. A typical deformed shape of the mesh for this motion under a prescribed deformation of is shown in Fig. 6. Additional simulations were conducted for strategies 1 and 2 to decide the maximal size of motion that they can perform without causing mesh tangling. This was found to be approximately 3=4 for both strategies. Although strategy 6 survived in the large displacement case, it resulted in the lowest value in the worst ratio. The extremely low value (about 0.01) indicates a severely deformed element. Strategies 3 and 4 performed similar to strategy 6 on this account. In contrast, strategy 5 was once again shown to have an advantage over other methods in preventing extremely low aspect ratios. One anomaly of the worst ratio plot in Fig. 7 is the low value for strategy 2 corresponding to a small displacement of just =5. This was found to be the result of a non-converged iteration in the iterative solver (GMRES). Attempts were made to reach a converged result for this case by increasing the number of inner iterations in GMRES to 200. Ultimately, this did not succeed either. Since a failure in the iterative solver usually indicates a lose of diagonal dominance in the matrix, this anomaly suggests that strategy 2 can potentially lead to ill-conditioning in the sti5ness matrix.
Z. Xu, M. Accorsi / Finite Elements in Analysis and Design 40 (2004) 1259 – 1269
1269
Worst r/R ratio -- Motion III
Average r/R Ratio -- Motion III
0.3
0.45 Strategy 1
0.43
Strategy 2
0.25 Strategy 1
0.2
0.41
Strategy 3
0.39
Strategy 4
0.15
Strategy 3
Strategy 5
0.1
Strategy 4
0.37
Strategy 6
Strategy 2
Strategy 5
0.05
Strategy 6
0.35 0
0.2
0.4
0.6
0
(a) (b)
0
0.2
0.4
0.6
Fig. 7. (a) Mesh quality for motion III—average r=R ratio, (b) mesh quality for motion III—worst r=R ratio.
4. Conclusions The mesh updating methods proposed in this research have been shown to be able to eliminate extremely deformed elements while maintaining overall mesh quality in mesh-moving problems for FSI simulations. Although di5erent methods perform similarly for small displacements, the approaches presented here exhibit an advantage over other methods in avoiding low element aspect ratios for prescribed large displacements of the structural surface. While the sti5ening method based on element distortion alone is very robust, the removal of deviatoric stresses and a negative Poisson’s ratio are clearly very bene!cial. Acknowledgements This paper is based upon work supported by, or in part by, the US Army Research OQce under Grant number DAAD19-99-1-0235. The authors also acknowledge support under a DoD HPC Challenge Project. References [1] A. Masud, A space–time !nite element method for &uid–structure interaction, Ph.D. Dissertation, Stanford University, 1993. [2] A. Johnson, T. Tezduyar, Mesh update strategies in parallel !nite element computation of &ow problems with moving boundaries and interfaces, Comput. Methods Appl. Mech. Eng. 119 (1994) 73–94. [3] T. Belytschko, W. Liu, B. Mora, Nonlinear Finite Elements for Continua and Structures, Wiley, New York, 2000. [4] G. Chiandussi, G. Bugeda, E. Onate, A simple method for automatic update of !nite element meshes, Commun. Numer. Methods Eng. 16 (2000) 1–19. [5] P. Fleischmann, Mesh generation for technology CAD in three dimensions, Ph.D. Dissertation, Technical University, Vienna, Austria, 2000. [6] K. Stein, Simulation and modeling techniques for parachute &uid–structure interactions, Ph.D. Dissertation, University of Minnesota, 1999. [7] Y. Saad, M. Schultz, GMRES: a generalized minimum algorithm for solving nonsymmetric linear systems, SIAM J. Sci. Statist. Comput. 7 (1986) 856–869.
Finite Elements in Analysis and Design 40 (2004) 1259 – 1269 www.elsevier.com/locate/!nel
Finite element mesh update methods for &uid–structure interaction simulations Zhenlong Xu, Michael Accorsi∗ Department of Civil and Environmental Engineering, University of Connecticut, Storrs, CT 06269-2037, USA Received 14 November 2002; received in revised form 24 April 2003; accepted 19 May 2003
Abstract Mesh updating is a problem commonly encountered in &uid–structure interaction simulations when structures undergo a large displacement. In this paper, a set of mesh updating strategies based on a “pseudo-solid” model is proposed. Generally, a homogeneous analysis is !rst performed on the pseudo-solid with prescribed boundary displacements. Di5erent element properties and external forces are then introduced based on strain results from the !rst analysis. A second analysis is performed on this non-homogeneous model with the goal to maintain element aspect ratio while preserving element volumes. Several numerical simulations are conducted using di5erent combinations of strategies proposed in this paper. In general, the proposed strategies signi!cantly enhance the performance of the pseudo-solid mesh updating method. ? 2003 Elsevier B.V. All rights reserved. Keywords: Fluid–structure interaction; Mesh update; Finite element method
1. Introduction In &uid–structure interaction (FSI) simulations, it is quite common for the structure to undergo large displacements. To maintain compatibility between the structural and &uid meshes, re-meshing is usually performed on the &uid mesh to accommodate this change. However, for large-scale simulations, this process is computationally expensive. An alternative approach to achieve this goal is mesh moving, in which the current &uid mesh is moved along with the deforming structure so as to preserve the compatibility between these two meshes. The most common approach for mesh moving is the “pseudo-solid” method. In this method, the &uid mesh is viewed as a linear elastic pseudo-solid, thereby turning the mesh-moving problem ∗
Corresponding author. Tel.: +1-860-486-5642; fax: +1-860-486-2298. E-mail address: [email protected] (M. Accorsi).
0168-874X/$ - see front matter ? 2003 Elsevier B.V. All rights reserved. doi:10.1016/j.!nel.2003.05.001
1260
Z. Xu, M. Accorsi / Finite Elements in Analysis and Design 40 (2004) 1259 – 1269
into a linear solid mechanics problem, with the displacement of the actual structural mesh acting as boundary conditions. Because the sole purpose of the pseudo-solid model is to move the &uid mesh, the solution to this problem itself has no physical importance. Therefore, di5erent methods can be employed to manipulate the properties of the pseudo-solid so as to achieve a mesh with minimum element distortion and volume change. Strategies based on the pseudo-solid method generally fall into two categories, one-step and two-step methods. In one-step methods, sti5ness is assigned to di5erent elements according to their initial geometric properties. Typical criteria include element sizes [1–3] or their distance to a boundary of prescribed displacements. Two-step methods are based on a predictor–corrector philosophy. In the !rst step (predictor), a linear structural analysis is conducted on the pseudo-solid mesh based on a homogeneous material model subjected to the prescribed boundary displacements. Element strains obtained from this step are then utilized as the basis for assigning sti5ness to di5erent elements. A second analysis (corrector) is then conducted on this non-homogeneous solid model and results of this analysis will be used as the desired motion of the &uid mesh. The success of two-step methods depends on the selection of element sti5ness based on the !rst analysis. A study examining various two-step methods has been presented by Chiandussi et al. [4]. In this paper, several new strategies for two-step methods are proposed based upon identi!cation of key reasons for element distortion. These methods can be either implemented independently, or combined with other approaches. Numerical simulations are then conducted on a 2-D mesh to compare the e5ectiveness of the di5erent strategies. Generally, the new methods proposed in this research are shown to be highly e5ective in preventing extremely distorted elements while maintaining the overall quality of the &uid mesh. The size of the required mesh motion corresponds to the motion of the structure during the FSI simulation, which is a function of many parameters. The goal of a mesh moving procedure is to control element distortion and maintain a high-quality mesh for as large a motion as possible. For large motions, it is anticipated that mesh moving will become less e5ective and re-meshing will be required. Therefore, in the numerical examples, the e5ectiveness of the mesh moving procedures is evaluated over a large range of mesh motion sizes. 2. Methodology 2.1. De1nition of mesh quality Since the ultimate goal of mesh moving is to achieve a high-quality !nite element mesh subjected to prescribed structural motions, a general measure of mesh quality should be de!ned. The most important factor that de!nes mesh quality is the aspect ratio of an element. It is generally acknowledged that elements with large aspect ratios perform poorly in !nite elements analyses. Therefore, a successful mesh moving strategy should not only maintain the overall quality of the mesh, but also be able to limit, if not avoid totally, the number of elements with extremely large aspect ratios. In this paper, for a typical 3-D tetrahedral element, we adopt the r=R value as the indicator of the element aspect ratio, in which r and R are the radii of the inner- and circum-sphere of the element, respectively. For 2-D analysis, r and R are the radii of the inner- and circum-circle of a triangular element [5].
Z. Xu, M. Accorsi / Finite Elements in Analysis and Design 40 (2004) 1259 – 1269
1261
2.2. Pseudo-solid formulation The basic !nite element formulation for the pseudo-solid system is K1 d1 = 0 d 1 = d0
(); ( = s )
(1)
in which K1 is the sti5ness matrix, d1 the mesh motion to be determined, and d0 the boundary displacements prescribed by the motion of the actual structural mesh. If this analysis serves as the !rst step in a two-step method, K1 is usually based on a homogeneous material model for the entire mesh. If this formulation is for a one-step method, K1 is generally non-homogeneous in the sense that di5erent elements may have di5erent material models (i.e. variable sti5ness) to control their deformation. This will be further explained in the next few sections. For a typical two-step method, a second analysis is also performed on the pseudo-solid mesh with designated element sti5ness, which is assigned to each element according to their strain obtained from the !rst step. The formulation for this step is K2 (d1 )d2 = F2 (d1 ) d 2 = d0
( = s )
(); (2)
in which d2 is the motion to be determined. The sti5ness matrix is written as K2 (d1 ) to re&ect its dependence on the !rst analysis. Unlike the homogeneous analysis, external forces could also be introduced into the second analysis to eliminate local deterioration of mesh quality so long as the boundary condition is satis!ed. Since the external force is also related to results from the !rst analysis, it is written as F2 (d1 ) in the above equation to re&ect this in&uence. Several strategies can be utilized in selecting K2 (d1 ) and F2 (d1 ). They will be explained in the following sections. 2.2.1. Poisson’s ratio For isotropic solids, the value of Poisson’s ratio re&ects the in&uence of the motion in one direction on the motion in the other two perpendicular directions. The range for its value is −1 ¡ ¡ 0:5 based on the requirement that the material has a positive bulk modulus and a positive shear modulus. Poisson’s ratio values for typical solids are positive. For these materials, compression in one direction will produce elongation in the other two directions, and vice versa. However, this kind of behavior has a negative e5ect on the geometry of elements in that elongation in one direction and shortening in the other two leads to degradation of the element aspect ratios. As mentioned previously, results obtained from the pseudo-solid model have no physical significance. Therefore, material properties can be chosen solely for the purpose of conserving element geometry without concern about whether such materials exist in the physical world or not. For this reason, use of a negative Poisson’s ratio is proposed in this research. The e5ect of this adoption is apparent. For example, with such a property, an element subjected to compression in one direction would tend to shorten in the other two perpendicular directions, thereby preserving the initial aspect ratio. However, this approach does have certain disadvantages in preserving element volume. Using the same example, it is apparent that a negative Poisson’s ratio may cause elements to contract too much under compression or elongate too much under tension.
1262
Z. Xu, M. Accorsi / Finite Elements in Analysis and Design 40 (2004) 1259 – 1269
Excessive change in element volume can be prevented by combining this method with other element sti5ening methods, as described in the following sections. 2.2.2. Element sti6ening Non-homogeneous element sti5ening is the primary idea underlining almost all mesh-moving methods based on pseudo-solid models. The di5erence between di5erent methods lies only in the selection of sti5ening functions, which usually relate to element geometry for one-step methods and element deformation for two-step methods. A typical one-step and two-step sti5ening scheme will be brie&y reviewed along with a new sti5ening scheme proposed in this research. The two former sti5ening schemes are given here primarily for the purpose of comparing results. The element sti5ness matrix for the pseudo-solid elements can be expressed as Ke = BT DB |J |e e d; (3)
where B is the derivative matrix of shape functions, D the constitutive matrix, |J |e the Jacobian of the element, and e the element sti5ening factor for mesh moving. Masud’s one-step sti6ening scheme: As a typical method in the one-step category, Masud [1] introduced a sti5ening scheme based on element sizes, i.e. choosing e as |J |max − |J |min ; (4) e = 1 + |J |e where |J |max and |J |min are the maximum and minimum Jacobian among all elements, respectively. This method essentially increases the sti5ness of small elements to prevent them from being distorted by large motions. Chiandussi’s two-step sti6ening scheme: Chiandussi et al. [4] proposed several two-step sti5ening schemes based on di5erent measures of element deformation, which were shown to have similar performance. Due to this reason, only one of those methods is presented here. It is a sti5ening scheme based on the square norm of element strains, i.e. 1 12 + 22 + 32 (5) N 2 3 in which 1 ; 2 and 3 are the three principal strains for an element based on results calculated from the !rst homogeneous analysis, and N is a user speci!ed strain constant for that element. Since, in practice, the same N is typically used for all the elements, this value has no e5ect on the !nal results for a linear analysis. Two-step sti6ening scheme based on element distortion: In this section, a sti5ening method is introduced based on the observations of the causes of element distortion. The key factors that a5ect element geometry are the maximum and minimum principal strains. The reason for this lies in the notion that the quality of an element’s geometry (i.e. aspect ratio) is essentially decided by its largest and smallest dimensions. Therefore, for an element with a near-optimal aspect ratio, which should be the case for most of the elements in a well-constructed initial mesh, the deformed shape should be controlled by the two extreme principal strains. While the in&uence of the middle principal strain on element geometry exists, its e5ect is not as essential as the other two. This criterion is noticeably di5erent than those introduced by Chiandussi et al., where all three principal strains are considered to have the same e5ect on element geometry. e =
Z. Xu, M. Accorsi / Finite Elements in Analysis and Design 40 (2004) 1259 – 1269
1263
Based on this notion, the element-sti5ening factors are chosen to be functions of (1 − 3 ), which can also be interpreted as the diameter of the largest Mohr’s circle for a 3-D strain state. Obviously, larger values of (1 − 3 ) indicate more severe distortion of elements, which, in turn, require higher sti5ening factors. In practice, many di5erent functions could be selected, ranging from linear, quadratic to exponential, all with (1 −3 ) acting as the only variable. For this research, the following exponential formulation is chosen: e = a(1 −3 ) ;
(6)
where a is a constant to be decided. In9uence of tension and compression: A very important but largely neglected factor that a5ects element geometry is the di5erent e5ects of tension and compression. In this research, it is recognized that compression generally tends to have a more severe e5ect on the degradation of element quality than tension. This can be demonstrated through a very simple example. For a 2-D element with an optimal initial aspect ratio of , assume the two principal strains are 0 and 0 after a certain deformation. If 0 is tensile (i.e. positive), the new aspect ratio for this element will be roughly (1 + 0 ). If 0 is compressive (i.e. negative), then the new aspect ratio will be approximately =(1 + 0 ). For a mesh-moving problem, it is not uncommon to have very large strains caused by large motions of the structural mesh. Given a moderately large strain value of ±0:5, the !nal aspect ratio would be 1:5 for the tension case (0 = 0:5) and 2 for the compression case (0 = −0:5). Clearly, compression has caused signi!cantly more damage to the aspect ratio than tension. This e5ect would be even more pronounced for extremely large strains. Therefore, any sti5ening technique that neglects this di5erence would be discriminating against elements under compression and potentially lead to a mesh with more distortion. In this research, the in&uence of tension/compression is accounted for by introducing a modifying factor 1=(1 + 3 ) to the element sti5ener e , which now becomes a(1 −3 ) e = : (7) 1 + 3 Deviatoric stress removal: As previously expressed in the basic formulation, external forces can also be introduced to eliminate element distortion. Generally, this is not implemented in the !rst step due to the uncertainty in guessing external forces that are bene!cial. However, this technique can be utilized in the second step based on results obtained from the !rst step. It is well known that a stress state can be decomposed into two parts: the hydrostatic and deviatoric stress. While the hydrostatic stress is responsible for a change in the element volume, the deviatoric stress is related to the element distortion. For a linear element with no displacement constraints, the release of the deviatoric stress results in a purely hydrostatic stress state with no element distortion. The consistent element force needed to remove the hydrostatic stress is F2 (d1 ) = − BT D (d1 )|J | d; (8)
where D (d1 ) is the deviatoric stress based on the results from the !rst analysis. As previously identi!ed by Stein [6], the existence of boundary constraints prevents elements from achieving a distortion-free state. This is also true for the proposed method of stress removal. In the second part of the two-step method, a complete removal of deviatoric stresses cannot be achieved
1264
Z. Xu, M. Accorsi / Finite Elements in Analysis and Design 40 (2004) 1259 – 1269
due to this in&uence. However, our assumption is that even a partial removal of these stresses will be bene!cial in eliminating local element distortion. It should also be noted that this approach is based solely on linear mechanics, whereas the deformation of the pseudo-solid, in general, can be quite large and possibly nonlinear. This is expected to also limit the removal of deviatoric stress in the elements. 3. Numerical results The strategies discussed above have been applied to a 2-D example, which is the same as used by Stein [6]. In this example, the mesh occupies the region de!ned by |x| 6 1:0 and |y| 6 1:0, as shown in Fig. 1. An interior surface is located in the middle of the mesh, speci!ed as |x| 6 0:5 and y = 0:0. In general, the strategies can be applied to 3-D models without any additional diQculties. Three types of motion of the inner surface were studied in this paper: (1) a uniform vertical translation in the y direction, (2) a rotation with respect to the center of the mesh, and (3) a deformation of the inner surface from a line to a circular arc with no stretching. For each motion, the three di5erent sizes given in Table 1 were used to study the relative e5ectiveness of the methods
Fig. 1. (a) Initial mesh, (b) initial mesh (enlarged).
Table 1 Motion types and sizes Relative scale Motion I (translation) Motion II (rotation) Motion III (deformation)
0.1 0.1 =20 =5
0.2 0.2 =10 2=5
0.5 0.5 =4
Z. Xu, M. Accorsi / Finite Elements in Analysis and Design 40 (2004) 1259 – 1269
1265
Table 2 Six strategies for mesh moving Strategy Strategy Strategy Strategy Strategy Strategy
1 2 3 4 5 6
Masud’s method (Eq. (4)) Chiandussi et al. two-step sti5ener (Eq. (5)) Two-step sti5ener with tension/compression in&uence (Eq. (7)) Strategy 3 + deviatoric stress removal (Eq. (8)) Strategy 4 + negative Poisson’s ratio of −0:5 Masud’s method with two half-steps
in handling small, medium, and large motions. For uniformity, all motion sizes are reported on a relative scale ranging from 0 to 0.5. The strategies used in this study are listed in Table 2. Considering the fact that a two-step method is twice as expensive, computationally, as a one-step method for the same motion size, strategy 6 is included for comparison where Masud’s one-step method is used with two half-steps rather than just one. Except for strategy 5, the Poisson’s ratio used in all other strategies has a value of 0.3. Young’s modulus of the pseudo-solid can be chosen arbitrarily in all cases. Also, a = 2 is used for the parameter in Eq. (7). For strategies 3–5 a homogeneous solid is used for the !rst step, and the modi!cations listed are applied in the second step. As described in the methodology, the ratio between the radii of the circum- and inner-circle is used as the criterion for mesh quality. For each motion, the average ratio for the entire mesh is used to evaluate the overall quality of the mesh, and the worst case (smallest ratio of all elements) is also presented to evaluate the worst degradation of mesh quality. For the initial mesh, the average ratio is 0.4401 and the worst ratio is 0.2629. The computer program used for this study is based on matrix-free techniques for formulating the system equations and an iterative solver GMRES [7] was adopted to solve those equations. All the results calculated in this research were obtained by using 100 iterations in the inner iterative loop of GMRES. 3.1. Motion I—vertical translation Fig. 2 shows a successful mesh motion for a prescribed vertical displacement of 0.5. In Fig. 3, the average and worst ratios for di5erent mesh movement sizes using di5erent methods are presented. Obviously, the overall mesh quality deteriorates as the mesh motion becomes larger. Although the di5erent methods appear to have similar performance in this aspect, the di5erence between them tends to grow, as the prescribed displacements get larger. Because the displacement increment used in strategy 6 is actually only half the size used in the other strategies, strategy 6 is able to achieve a better overall quality than all other methods. Under large displacements, the worst ratio value drops dramatically, from around 0.23 to about 0.1. The three strategies proposed in this research (i.e. 3–5) perform noticeably better than other methods for large displacement. With the help of deviatoric stress removal, strategy 4 is slightly better than strategy 3, and strategy 5 outperforms strategy 4 by additionally utilizing a negative Poisson’s ratio. In practice, controlling the worst case is most critical because failure of one element (i.e. element tangling) is the limiting factor for the mesh-moving process.
1266
Z. Xu, M. Accorsi / Finite Elements in Analysis and Design 40 (2004) 1259 – 1269
Fig. 2. (a) Motion I—vertical translation, (b) motion I—vertical translation (enlarged).
Worst r/R Ratio -- Motion I
Average r/R Ratio -- Motion I
0.3
0.45 Strategy 1 Strategy 2 Strategy 3 Strategy 4 Strategy 5 Strategy 6
0.43 0.41 0.39 0.37
0.2 0.15 0.1 0.05 0
0.35 0
(a)
Strategy 1 Strategy 2 Strategy 3 Strategy 4 Strategy 5 Strategy 6
0.25
0.2
0.4
0.6
0
0.2
0.4
0.6
(b)
Fig. 3. (a) Mesh quality for motion I—average r=R ratio, (b) mesh quality for motion I—worst r=R ratio.
A very important property for this translational motion is that it roughly divides the pseudo-solid mesh into two regions: the compression and tension part, with a small fraction of elements lying between these two areas. From Fig. 2, it is quite obvious that elements in the compression region deform much more severely than those in the tension region. This observation validates the reasoning for using a sti5ening scheme that accounts for the e5ect of tension and compression. In Fig. 3, strategy 2 has surprisingly low values for the worst-case ratio. This can be explained by a more detailed look at the deformed shape, which is not shown here. The worst-case element is not located near the tips of the structural surface, as previously assumed, but rather on the border of the mesh somewhere between the tension and compression region. A plausible explanation is that the sti5ening formulation in this strategy did not adequately protect elements with negative
Z. Xu, M. Accorsi / Finite Elements in Analysis and Design 40 (2004) 1259 – 1269
1267
Fig. 4. (a) Motion II—rotation, (b) motion II—rotation (enlarged).
Worst r/R Ratio -- Motion II
Average r/R Ratio -- Motion II
0.3
0.45 Strategy 1 Strategy 2 Strategy 3 Strategy 4 Strategy 5 Strategy 6
0.43 0.41 0.39 0.37
0.25
Strategy 1 Strategy 2 Strategy 3 Strategy 4 Strategy 5 Strategy 6
0.2 0.15 0.1 0.05
0.35 0
0.2
(a)
0.4
0.6
0 0
0.2
0.4
0.6
(b)
Fig. 5. (a) Mesh quality for motion II—average r=R ratio, (b) mesh quality for motion II—worst r=R ratio.
principal strains since it does not consider the stronger in&uence on element geometry caused by compression. 3.2. Motion II—rotation Fig. 4 shows a successful mesh motion for a prescribed rotation of =4. In Fig. 5, the average ratio and worst ratio for di5erent mesh movement sizes using di5erent methods are presented. The trend in average and worst values agrees with results listed in the previous section, i.e. mesh quality degrades with the increase in magnitude of the prescribed displacement. One noticeable di5erence between the deformed mesh for this motion and that under translation is that most elements in this case did not undergo a large change in volume, but rather a change in shape. This is apparent by comparing Figs. 2 and 4.
1268
Z. Xu, M. Accorsi / Finite Elements in Analysis and Design 40 (2004) 1259 – 1269
Fig. 6. (a) Motion III—deformation, (b) motion III—deformation (enlarged).
For this motion, strategy 6 outperforms the other methods on both measures. One possible explanation for this is that this strategy is based on the idea of preserving element volume, which happens to agree with the behavior of the deformed mesh under this motion. The methods proposed in this research performed fairly well, although not as good as strategy 6. Strategy 2 once again was shown to have the lowest value in worst ratio. 3.3. Motion III—deformation This case is the most demanding one among all the three motions considered, and strategies 1 and 2 actually fail for the largest displacement. A typical deformed shape of the mesh for this motion under a prescribed deformation of is shown in Fig. 6. Additional simulations were conducted for strategies 1 and 2 to decide the maximal size of motion that they can perform without causing mesh tangling. This was found to be approximately 3=4 for both strategies. Although strategy 6 survived in the large displacement case, it resulted in the lowest value in the worst ratio. The extremely low value (about 0.01) indicates a severely deformed element. Strategies 3 and 4 performed similar to strategy 6 on this account. In contrast, strategy 5 was once again shown to have an advantage over other methods in preventing extremely low aspect ratios. One anomaly of the worst ratio plot in Fig. 7 is the low value for strategy 2 corresponding to a small displacement of just =5. This was found to be the result of a non-converged iteration in the iterative solver (GMRES). Attempts were made to reach a converged result for this case by increasing the number of inner iterations in GMRES to 200. Ultimately, this did not succeed either. Since a failure in the iterative solver usually indicates a lose of diagonal dominance in the matrix, this anomaly suggests that strategy 2 can potentially lead to ill-conditioning in the sti5ness matrix.
Z. Xu, M. Accorsi / Finite Elements in Analysis and Design 40 (2004) 1259 – 1269
1269
Worst r/R ratio -- Motion III
Average r/R Ratio -- Motion III
0.3
0.45 Strategy 1
0.43
Strategy 2
0.25 Strategy 1
0.2
0.41
Strategy 3
0.39
Strategy 4
0.15
Strategy 3
Strategy 5
0.1
Strategy 4
0.37
Strategy 6
Strategy 2
Strategy 5
0.05
Strategy 6
0.35 0
0.2
0.4
0.6
0
(a) (b)
0
0.2
0.4
0.6
Fig. 7. (a) Mesh quality for motion III—average r=R ratio, (b) mesh quality for motion III—worst r=R ratio.
4. Conclusions The mesh updating methods proposed in this research have been shown to be able to eliminate extremely deformed elements while maintaining overall mesh quality in mesh-moving problems for FSI simulations. Although di5erent methods perform similarly for small displacements, the approaches presented here exhibit an advantage over other methods in avoiding low element aspect ratios for prescribed large displacements of the structural surface. While the sti5ening method based on element distortion alone is very robust, the removal of deviatoric stresses and a negative Poisson’s ratio are clearly very bene!cial. Acknowledgements This paper is based upon work supported by, or in part by, the US Army Research OQce under Grant number DAAD19-99-1-0235. The authors also acknowledge support under a DoD HPC Challenge Project. References [1] A. Masud, A space–time !nite element method for &uid–structure interaction, Ph.D. Dissertation, Stanford University, 1993. [2] A. Johnson, T. Tezduyar, Mesh update strategies in parallel !nite element computation of &ow problems with moving boundaries and interfaces, Comput. Methods Appl. Mech. Eng. 119 (1994) 73–94. [3] T. Belytschko, W. Liu, B. Mora, Nonlinear Finite Elements for Continua and Structures, Wiley, New York, 2000. [4] G. Chiandussi, G. Bugeda, E. Onate, A simple method for automatic update of !nite element meshes, Commun. Numer. Methods Eng. 16 (2000) 1–19. [5] P. Fleischmann, Mesh generation for technology CAD in three dimensions, Ph.D. Dissertation, Technical University, Vienna, Austria, 2000. [6] K. Stein, Simulation and modeling techniques for parachute &uid–structure interactions, Ph.D. Dissertation, University of Minnesota, 1999. [7] Y. Saad, M. Schultz, GMRES: a generalized minimum algorithm for solving nonsymmetric linear systems, SIAM J. Sci. Statist. Comput. 7 (1986) 856–869.