- Email: [email protected]
A remeshing algorithm based on bubble packing method and its application to large deformation problems
Finite Elements in Analysis and Design 39 (2003) 301 – 324 www.elsevier.com/locate/ nel
A remeshing algorithm based on bubble packing method and its application to large deformation problems Soon Wan Chung, Seung Jo Kim ∗ Department of Aerospace Engineering, School of Mechanical and Aerospace Engineering, Seoul National University, San 56-1, Shinlim-dong, Kwanak-gu, Seoul 151-742, South Korea Received 15 April 2001; accepted 11 January 2002
Abstract The remeshing algorithm based on the bubble packing method (BPM) is developed for nite element analysis (FEA) because the optimal position of nodes can be obtained systematically when only the bubble size function is given. The triangular meshes are generated by Delaunay triangulation with the advancing front concept. To use the automatically generated mesh information in FEA, a new bandwidth minimization scheme with high e5ciency in CPU time is also developed. A re ning circle is proposed to specify the re nement area using the bubble size in the remeshing step. The remeshing algorithm is applied to the metal forming problem to replace the distorted mesh at the largely deformed area with regular and ne meshes. The numerical results show that the remeshing algorithm based on the BPM works well at the region with large error and is able to control the re nement area and the new mesh size easily through two parameters (∗max ; q). ? 2002 Elsevier Science B.V. All rights reserved. Keywords: Remeshing; Bubble packing method; Bandwidth minimization; Re ning circle; Large deformation
1. Introduction For several decades, adaptive remeshing by referencing the error estimation in nite element analysis (FEA) has been investigated actively, because mesh re nement is crucial in nonlinear problems, such as the metal forming process. Therefore, the remeshing algorithm based on bubble packing method (BPM), which is one of useful mesh generation algorithms, is developed and applied to large deformation problems in this paper. Generally, FEA including the remeshing is composed of the following three modules: automatic mesh generation, error estimation and adaptive remeshing. A number of automatic mesh generation ∗
Corresponding author. Tel.: +82-2-880-7388; fax: +82-2-887-2662. E-mail address: [email protected] (S.J. Kim).
0168-874X/02/$ - see front matter ? 2002 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 8 - 8 7 4 X ( 0 2 ) 0 0 0 7 5 - 6
302
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324
algorithms in the two-dimensional domain have been proposed such as the node connection approach [1–3] and the spatial decomposition procedure [4 – 6]. In the node connection approach, Delaunay triangulation or the advancing front method is generally used, and in the spatial decomposition method, the domain is subdivided recursively until the mesh size attains the desired level. Recently, the BPM [7] has been proposed, which has the advantages of the optimal node placement by the physically based relaxation and ease of mesh size control. In the eld of error estimation, the Zienkiewicz–Zhu error estimator [8,9] is popular due to simple implementation and cost eIectiveness. To estimate the discretization error, the recovered values are used in place of the exact solutions when the exact ones are not available in the practical problems. As the recovery technique, superconvergent patch recovery (SPR) [10,11] and recovery by equilibrium in patches (REP) [12] are representative. The adaptive remeshing is needed at the region where the solution gradient is high or when the initial mesh becomes distorted severely, for example, in metal forming problems. The remeshing criterion aIecting re nement area is generally based on error estimation [13–16] or the element distortion indicator [17]. In addition, the routine for accurate data transfer between an old mesh and a new mesh should be included. In this paper, the BPM is used in the remeshing algorithm [18] as well as mesh generation algorithm. As the BPM considers a node as the center of a bubble (a circle in two dimensions) and represents mesh size by bubble size, the re nement area and new mesh size in this area can be calculated from the bubble size function. We propose a re ning circle centered on the candidate element for re nement and calculate the new bubble size in the circle using control factor q. In the automatic mesh generation module, Delaunay triangulation with the advancing front concept is used to connect the nodes generated by the BPM. Also we adopt the Zienkiewicz–Zhu error estimator combined with SPR. As the numerical experiments, we rst investigate an elastostatic problem, the tension problem of a plate with an elliptical hole. Secondly, the metal forming process is simulated for the large deformation problem. The equilibrium equation based on the elasto-viscoplastic model and total Lagrangian formulation are used for the nonlinear nite element approximation. 2. Automatic mesh generation algorithm The automatic mesh generation in this study consists of two main parts: (1) optimal node placement through the BPM and (2) Delaunay triangulation with the advancing front concept. The algorithm of the BPM is based on Ref. [7] and is modi ed for simple implementation. The enhanced bandwidth minimization technique based on Ref. [19] is developed to reduce CPU time during the bandwidth minimization. 2.1. Bubble packing method This method determines the node placement for triangular elements in consideration of the equilibrium position by the attractive=repulsive interbubble forces of the nodes. The forces applied at the bubbles are calculated by a bubble size function. The close packing of bubbles is similar to a pattern of Voronoi polygons and therefore, Delaunay triangulation of these bubbles yields well-shaped triangles close to the regular triangles. Since the BPM compacts the bubbles directly inside the
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324
303
Fig. 1. Procedure of node distribution based on bubble packing method: (a) nodes given in input data, (b) boundary line bubble packing, (c) Surface bubble packing-1, (d) Surface bubble packing-2, (e) Bubble distribution after dynamic simulation.
domain according to the bubble size function, it is noted that BPM does not have to interpolate the background mesh [20,21] to de ne the mesh size. The mesh size in adaptive remeshing is determined by the concept of a re ning circle and will be described in Section 4. The BPM consists of the following steps: 1. 2. 3. 4.
preparation of input data for the boundary of the computational domain, boundary line packing of bubbles, surface packing of bubbles, physically based dynamic simulation.
Step 1. Preparation of input data for the boundary of the computational domain (Fig. 1(a)): In order to use in automatic mesh generation, the input data provided by a user should be concise. Hence, nodes at the vertex, lines connecting the nodes, and surfaces composed of lines are written in the input le for mesh generation. The bubble size of a node is used in the following packing procedure and the dynamic simulation. A xed node is often necessary when boundary condition must be imposed at a speci c point in FEA. For this case, we make its bubble size a negative value for distinction. Step 2. Boundary line packing of bubbles: The bubbles are packed along the boundary lines connecting nodes given in the input data (Fig. 1(b)). At this time, the bubbles are placed compactly according to the bubble size function.
304
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324
Fig. 2. Delaunay triangulation with the advancing front concept.
Step 3. Surface packing of bubbles: The center of the domain is given in the input data and then the line connecting the center and the vertex is packed in the same way as in Step 2. We call this line diagonal line for surface packing (DLSP). After each DLSP is packed, bubbles are inserted at the line between nodes on diIerent DLSPs (Fig. 1(c)). The nodes generated in a hollow domain such as a void are removed by the vertical ray shooting algorithm (Fig. 1(d)) [22]. Step 4. Physically based dynamic simulation: The above procedures (Steps 1–3) are for the initial bubble placement. The bubble distribution from the procedures may make gaps or overlaps between bubbles. To overcome this problem, bubble con guration balanced with the interbubble attractive=repulsive forces is found in this step (Fig. 1(e)). The force function is de ned as a cubic function of the distance between the centers of two bubbles. The forces applied at boundary bubbles are projected to the boundary line to constrain the bubbles on the boundary. The xed node speci ed at Step 1 applies the force to other bubbles, but it does not take any forces by other bubbles. The equations of motion of each bubble are integrated by the fourth-order Runge–Kutta method [7]. After integration, a bubble is deleted when two bubbles are too close for their bubble sizes and when the bubble exists outside the domain. 2.2. Delaunay triangulation with the advancing front concept Triangular meshes are generated by connecting nodes obtained from the BPM, and for this procedure Delaunay triangulation with the advancing front concept is incorporated. As the triangulation is started from boundary edges and advances into the interior domain until the whole domain is covered completely, we carry out Delaunay triangulation using the largest interior angle property of the Voronoi nearest neighbors [23] (Fig. 2). Newly generated edges are considered as new boundaries and the nearest node is searched by the largest interior theorem. This node connection algorithm can be as e5cient as the point insertion method [24], if nodes are organized properly [23]. Moreover, it is easy to extend the algorithm in the case that the problem domain is multiply connected because the generation of fronts from the inner boundaries is straightforward. 2.3. Additional algorithms for checking node position The nodes lying outside the domain in Step 3 or 4 of the BPM are removed using the vertical ray shooting algorithm. In the Delaunay triangulation stage, the node search algorithm (NSA) on boundary edge is incorporated keeping in mind that the candidate node should be located at only one side. If we assume that the direction vectors (vx ; vy ) of boundary lines are in a counterclockwise
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324
305
Fig. 3. Sketch for the node search algorithm.
direction as in Fig. 3, the admissible region of the outer (inner) boundary can be determined from ◦ the normal vector obtained by rotating the direction vector 90 counterclockwise. Since the normal vector of a boundary line is (−vy ; vx ), the line equation is as follows: − vy x + vx y + c = 0:
(1)
If a node (x1 ; y1 ) is on the left side of the direction vector, − vy x1 + vx y1 + c ¿ 0;
(2)
and if a node is on the right side of the direction vector, − vy x1 + vx y1 + c ¡ 0:
(3)
As the admissible region for the outer=inner boundary is on the left=right side, the node satisfying inequality (2)=(3) is allowable as the candidate node for the boundary lines. 2.4. Recursive=adaptive bubble population control To increase the regularity of a triangular mesh, adaptive and recursive bubble population controls are used. Adaptive control: The number of surrounding meshes of each interior node is calculated, and if the average value of the numbers is smaller than 5.9 or larger than 6.1, bubble addition or deletion is performed, because the number for regular meshes is 6. If the number of surrounding meshes of a node is smaller than 5, a new bubble is inserted at the midpoint of the edge with maximum interior angle (Fig. 4(a)). If the number is larger than 7, a node which is adjacent to a neighboring mesh with a smaller interior angle than the other neighboring mesh ( ¡ ) is deleted between two bubbles on the edge with minimum interior angle (Fig. 4(b)). Recursive control: When the adaptive bubble control is used, this algorithm returns to the dynamic simulation (Step 4 of BPM) to nd a new equilibrium position of the bubbles. This recursion is repeated until the adaptive control is not activated or the number of repetitions reaches 5. 2.5. Bandwidth minimization Node and element numberings from automatic mesh generation are distributed randomly, and the maximum bandwidth of global stiIness matrix may be very large. Therefore, a bandwidth minimization technique is needed if the band solver is used to solve the simultaneous equations resulted by FE approximation. The bandwidth minimization technique in this study is based on Puttonen’s
306
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324
(a)
(b)
Fig. 4. Adaptive bubble control: (a) bubble addition and (b) bubble deletion.
bandwidth reduction algorithm [19]. The basic idea is to compact the rows and columns of the symmetric node connectivity matrix. If the node connectivity matrix is a square matrix, the required memory in the problem with 10 000 nodes is 100 Mbytes, for example. However, if the condensed connectivity matrix (CCM) which contains only nonzero column numbers on every row is used for row=column change, the required memory is at most 100 kbytes. In Ref. [19], the use of the CCM is commented on, but a full description is not presented. Hence, in this study the algorithm using the CCM is explained and the enhancement of CPU time by this algorithm is shown in several examples. Since Puttonen’s algorithm uses row=column exchanges, row=column exchanges in the CCM can be classi ed into four cases according to the existence of nodes at the corresponding row=column. Fig. 5 shows four cases taking place in the CCM when the ith column and the jth column are exchanged in any kth row. In Cases 1 and 2, no changes appear in the CCM, and the changes in the CCM in other cases are as follows: Case 3: CCM (k; 1) = i; CCM (k; 2) = i + 1; CCM (k; 3) = i + 3 ⇒ CCM (k; 1) = i + 1; CCM (k; 2) = i + 3; CCM (k; 3) = j, Case 4: CCM (k; 1) = i + 1; CCM (k; 2) = i + 2; CCM (k; 3) = j ⇒ CCM (k; 1) = i; CCM (k; 2) = i + 1; CCM (k; 3) = i + 2. To show the e5ciency of the CCM, the CPU time is compared with that by using square connectivity matrix in three models (Fig. 6) through a PC with an Alpha-533 MHz CPU (Table 1). From the
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324 i
307
j
Case 1
Case 2
Case 3
Case 4
Fig. 5. Four cases of row=column exchange for the CCM.
Fig. 6. Models used for comparison of CPU time during bandwidth minimization.
Table 1 Comparison of CPU time
Model 1 Model 2 Model 3
Number of nodes
Number of elements
Square connectivity matrix (s)
CCM (s)
Change in half-bandwidth
173 399 593
304 716 1104
0.141 0.640 3.080
0.068 0.139 0.425
129 → 22 273 → 14 534 → 52
results, CPU time with the use of CCM is found to be much smaller than the case without the CCM, and the e5ciency is more prominent in the model with the larger number of nodes. As the node number is reordered through the bandwidth minimization procedure, the extra storage for relating the node numbers before and after bandwidth minimization is added. The Oowchart of the above automatic mesh generation algorithm is presented in Fig. 7.
308
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324 Input Minimal geometric data, Bubble size function, Constants(m,c)
Bubble Packing Method
Boundary line packing of bubbles
Surface packing of bubbles
Physically-based dynamic simulation
Delaunay triangulation with advancing front concept
Yes (Recursive bubble control)
Is adaptive bubble population control used? No Bandwidth minimization : CCM
Mesh generation completed!
Fig. 7. Flowchart of automatic mesh generation algorithm.
3. Nonlinear nite element analysis To use the remeshing algorithm based on the BPM in the large deformation problem, we perform the nonlinear nite element analysis. A kinematic strain measure [25] utilizing polar decomposition of the deformation gradient and constitutive equations based on the theory of material of type-N [26] are used to simulate geometrical and materially nonlinear behaviors. To realize contact conditions, the extended interior penalty method [27] is adopted. The resulting highly nonlinear equilibrium equation is reduced to the incremental weak form and approximated by the total Lagrangian formulation.
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324
309
3.1. Kinematics Assuming the unrotated current con guration Ct − and an intermediate con guration Cp− ; U can be decomposed additively as follows. [25]: U = Ue + Up − I;
(4)
where Ue and Up represent the elastic stretch and the plastic stretch tensors, respectively. The velocity gradient tensor and the velocity strain tensor take the following forms: ˙ t + R(U˙ e U−1 + U˙ p U−1 )Rt ; ˙ −1 = RR L = FF
(5)
D = L|sym = R(U˙ e U−1 |sym + U˙ p U−1 |sym )Rt :
(6)
From Eqs. (5) and (6), one can de ne the elastic and plastic strain rates as E˙ = U˙ e U−1 |sym ;
P˙ = U˙ p U−1 |sym :
(7)
3.2. Elasto-viscoplastic material model The theory of materials of type N consists of ve eld equations (conservation of mass, balance of linear and angular momenta, conservation of energy, law of entropy production) and the following two potentials: • Free energy functional: 1 1 1 2 2 {(tr E) + 2(tr E )} + h1 z + (h1 − h0 )exp(−mz) : = 0 2 m • Flow potential functional: = D0
∞ i=0
(−1)
i+1
i h2i−1 i!(2i − 1)J2i−1=2
1 where = n 3
The nal constitutive equations can be obtained as follows: S = [(tr E)I + 2E]; 0 h = [h1 + (h0 − h1 )exp(−mz)]; 0
h2 1 ˙ √ P = D0 S ; exp − J 2h J2 2 1 ˙ z˙ = S : P: h
n+1 n
(8)
:
(9)
(10) (11) (12) (13)
To integrate the above constitutive equations (12), (13) with respect to time, the predictor–corrector method for a hardening variable (z) and the backward Euler method for plastic strain (P) are used.
310
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324
4. Adaptive remeshing algorithm Whenever FEM is used to solve a diIerential equation, discretization error is introduced and the distorted mesh due to the large deformation process, in particular, may induce large approximation errors. In order to decrease this error, the remeshing procedure is performed over the whole domain and the re ned meshes are generated at the element with large error in the present study. For the re nement criterion, the maximum permissible error and the bubble size control factor are used. The mesh re nement can be implemented easily by reducing the bubble size in the BPM because the mesh size is determined by bubble size. Since node regeneration by the BPM is started at the boundary and the nodes may be added at the boundary lines during the re nement stage, the update of boundary information is important. The extra storage for boundary nodes commented in Section 2.5 is needed in this re nement stage of boundary lines in order to identify the boundary nodes among the nodes that resulted from the bandwidth minimization technique. After the boundary update is nished, the mesh generation procedure described in Section 2 is performed again. The Oowchart of this remeshing algorithm is shown in Fig. 8. 4.1. Recovery technique and error estimation To estimate a posterior error, the nodal point value is computed from the Gauss point value by the recovery technique and the recovered value obtained by interpolating the nodal point value is
Finite Element Analysis at (i) - th increment step
Error estimation
Refinement criterion
i=i+1
Is remeshing necessary?
Yes Automatic mesh generation
Mapping of state variables (stress, strain, displacement, et al)
Fig. 8. Flowchart of remeshing algorithm.
No i=i+1
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324
311
Node determined by recovery procedure Sampling point selected firstly Sampling point selected secondly(either of)
Fig. 9. Recovery at boundary nodes.
compared with the nite element solution. As the recovery technique, SPR [10] is adopted in this research. For the linear triangular element, the linear polynomial of strain is used as shown below: !1∗ = Pa;
(14) T
where P = [1; x; y] and a = [a1 ; a2 ; a3 ] . The unknown parameters a of Eq. (14) are calculated by a least-squares t of the set of superconvergent points. After minimizing the diIerence between the discontinuous strain (!) ˆ derived by nite element analysis and the linear expansion !1∗ , we can obtain the following matrix form: (15) a = A−1 b; n n T T where A = i=1 P (xi ; yi )P(xi ; yi ) and b = i=1 P (xi ; yi )!(x ˆ i ; yi ) and n is the number of sampling points. As the number of unknowns is 3, the number of sampling points must be 3 or more. To prevent the abrupt increase of recovered values at the boundary, we select the number of sampling points to be more than 3. If the number is less than 4 at the boundary nodes, the neighboring elements connected with any two nodes on the element patch are added (Fig. 9). After a is known, the nodal parameter (!S∗ ) is found out by substituting the coordinates of the node in Eq. (14), and the smoothed strain eld is interpolated by the nodal parameter. !∗ = N !S∗ :
(16)
Here, the interpolation function N is the same as the one used for the interpolation of displacements in FEM. The strain error using the recovered value is de ned by e! ≡ !∗ − !:ˆ The L2 norm of strain error e! and the predicted relative error can be written as m m
e! 2 ∗ e! = e! 2i = (!∗ − !) ˆ 2i d%; = ; !∗ 2 + e! 2 %i i=1 i=1
(17)
(18)
312
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324
Refining circle d1 2
d1+d2+d3 3
d3 2 d2 2
Candidate element for refinement Fig. 10. Re ning circle de ned by bubble size.
where m is the total number of elements. If a maximum permissible relative error (∗max ) is provided, the maximum permissible error in any element i can be de ned as shown [16]:
1=2 ∗ e! i; max = √max !∗ 2 + e! 2 : (19) m 4.2. Re9nement criterion In this research, since a re nement is provided at the element with error more than the maximum permissible error in an element ( e! i; max ), only the elements with high error are re ned as ∗max of Eq. (19) is larger. In other words, if ∗max is large, the re ned region is narrow and the increase in number of elements after re nement is also small. Therefore, ∗max can be chosen according to the error distribution of the example, and we will compare the re ned meshes according to the variation of ∗max in the linear elastostatic problem. After an element with large error is found, the re nement region and the new mesh size within the region should be systematically speci ed. A re ning circle is introduced in this algorithm to indicate the re nement region centered on the candidate element for re nement (Fig. 10). The center corresponds to the centroid of the candidate element for re nement and the radius is taken as the average bubble size of three vertices. The bubble size of a new node, which will be placed within this re ning circle, is calculated as follows: e! i; max q dnew = dave : (20) e! i Here, dave is the average bubble size of three vertices, and q is the factor to control the new bubble size dnew . When q ¿ 1, the reduction ratio of bubble size (dave −dnew )=dave becomes large, and when q ¡ 1, the ratio is relatively small. That is to say, the re nement region is determined by ∗max and the bubble size generated in the re ning circle can be controlled by q. Therefore, the new bubble
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324
313
size distribution by two parameters (∗max ; q) in addition to the updated boundary information are essential to this remeshing algorithm using BPM. 4.3. Data transfer When the remeshing stage is performed, the data transfer should be treated carefully. In the linear elastic problem, the transfer of a boundary node, at which boundary conditions or loading conditions are imposed, is only necessary because the loading is applied again at the undeformed con guration. The information at the boundary node is particularly important in this algorithm, because the automatic mesh generation is started at the boundary nodes and new nodes may be inserted at the boundary lines for two cases in the remeshing stage. One case is where a new node is inserted at the middle of the boundary line when the candidate element for re nement is adjacent to the boundary line. The other case is where new nodes are inserted during Step 2 of the BPM (boundary line packing of bubbles), when the distance between two nodes is greater than the sum of two bubble diameters due to the large deformation. If the boundary conditions of two end nodes on the boundary line agree, the same condition is applied at the new node. In the remeshing stage of a nonlinear problem, not only the information at the boundary node but also the state variable such as displacement, stress and strain should be considered. A Gauss point value such as stress or strain is transferred by the following three steps, and a nodal point value such as displacement is transferred by only Step 2 [15]. (1) Transformation of Gauss point value into nodal point value in old mesh by SPR, (2) Mapping between nodal point values in old and new meshes, (3) Transformation of nodal point value into Gauss point value in new mesh. In step 2, we rst use NSA in order to search an old element to which a new node belongs. When we consider three sides of an old element as outer boundaries, a new node exists in the element if the node is located in the admissible domain, that is to say, the node satis es Eq. (2) for all three lines. The new nodal point value is calculated by interpolating the old nodal point value using the area coordinates. In step 3, the Gauss point value is the average of three nodal point values because one point Gauss quadrature is utilized. In Ref. [17], it is addressed that these steps 2 and 3 procedures are superior to the direct transfer when the nodal value of an old mesh is transferred into the Gauss point of a new mesh. The set of contact candidate nodes is also considered in the upsetting problem when the information at the boundary node is transferred. If a new node added at the boundary line is between the contact candidate nodes, the node is included in the set of contact candidate nodes, and the contact force applied at the node is calculated.
5. Numerical experiments 5.1. Linear elastic analysis We consider a linear elastostatic problem, the tension of a square plate with an elliptical hole under the plane strain assumption. The geometry of the model is shown in Fig. 11. Young’s modulus and
314
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324
Fig. 11. Geometric description of square plate with elliptical hole.
Fig. 12. Meshes generated by uniform re nement in the tension problem: (a) mesh-1 (bubble size = 0:4), (b) mesh-2 (bubble size = 0:2), (c) mesh-3 (bubble size = 0:1).
Poisson’s ratio are E = 105 ; ( = 0:3, respectively, and the error in the equivalent strain (!eq = 2 ! : !ij ) is used for the error estimation. 3 ij At rst, three meshes are generated with uniform bubble sizes: 0.4, 0.2 and 0.1, respectively. The resulting meshes are regular as shown in Figs. 12(a) – (c) and the numbers of elements are 250, 1031 and 4264, respectively. The error ( e! i ) distributions of the three meshes are presented in Fig. 13, and the error decreases as the element size decreases. Since the error at point (A) in Fig. 11 is larger than any other region, it can be expected that more re ned meshes will be necessary at point (A). The converged equivalent strain is obtained by the extrapolation based on the following equation using the strains for the three uniform bubble sizes [28]. e 6 Chmin(p; ) ≈ CN −0:5 min(p; ) :
(21)
Here, h is the size of the element, N the number of degrees of freedom, C a positive constant depending on the problem, p the polynomial degree of shape function, and the strength of the
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324
315
Fig. 13. Error (e! i ) distribution of uniformly re ned meshes in the tension problem: (a) mesh-1, (b) mesh-2, (c) mesh-3.
singularities. Finally, the relative error in the L2 norm is calculated by 1=2 ! 2 − ! e! ˆ 2 = = : ! ! 2
(22)
For the adaptive re nement, Mesh-1 is used as the initial mesh and remeshing is continued until the relative error () becomes less than 5%. We take account of ve ∗max ’s (10%, 9%, 8%, 7%, 6%) for q = 1. The re ned meshes are generated around point (A) as shown in Fig. 14, and the mesh size is reduced as the re nement progresses. When ∗max is small, the number of elements generated during each remeshing step increases. The mesh quality de ned by the following equation is calculated and the averages of mesh qualities are summarized in Table 2. r (23) Q=2 ; R where r is the radius of the inscribed circle and R the radius of the circumscribed circle. The comparison of the convergence rates in uniform and adaptive re nements is shown in Fig. 15. Obviously, the rate of convergence in adaptive re nement is higher than that in uniform re nement,
316
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324
Fig. 14. Meshes generated by adaptive re nement in the tension problem (Mesh-r{∗max }%-{remeshing order} (number of elements after remeshing), q = 1): (a) Mesh-r10-1 (348), (b) Mesh-r10-2 (499), (c) Mesh-r10-3 (664), (d) Mesh-r10-4 (776), (e) Mesh-r8-1 (437), (f) Mesh-r8-2 (734), (g) Mesh-r8-3 (991), (h) Mesh-r6-1 (683), (i) Mesh-r6-2 (1299).
Table 2 Average quality of meshes generated by remeshing in the tension problem Re nement criterion
First remeshing
Second remeshing
Third remeshing
Fourth remeshing
Fifth remeshing
∗max = 10%; q = 1 ∗max = 9%; q = 1 ∗max = 8%; q = 1 ∗max = 7%; q = 1 ∗max = 6%; q = 1
0.912 0.912 0.900 0.903 0.910
0.910 0.891 0.903 0.909 0.911
0.907 0.908 0.905 0.887
0.911 0.897 0.918
0.900
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324
317
Relative error (η) (%)
30
10
Refinement type [Uniform and * , q)] Adaptive(ηmax Uniform 10%, 1 9%, 1 8%, 1 7%, 1 6%, 1
1 100
1000
5000
Number of degree of freedom
Fig. 15. Convergence rates of relative error(%) in the tension problem.
and the tendency is consistent with respect to ∗max . From these results, it can be said that this remeshing algorithm works well at the region with high error and the remeshing=re nement procedure using the BPM is eIective in controlling the re ning area and the mesh size. 5.2. Nonlinear elasto-viscoplastic analysis As a representative nonlinear problem, we consider the metal forming process, that is, the upsetting problem. The numerical model is an Al 2024 axisymmetric billet consisting of 16 mm radius and 40 mm height, and the boundary and loading conditions are shown in Fig. 16. If the solution does not converge within the given iteration number in an increment step, we halve the increment of die movement to search the equilibrium position under a smaller loading increment. The error in equivalent strain is used for the error calculation. For the uniform re nement, three meshes are generated with the uniform bubble sizes 4, 2 and 1 mm in the BPM (Fig. 17(a-1,-2,-3)) and the resulting numbers of elements are 100, 384 and 1512, respectively. From the deformed meshes of Mesh-1–3 (Fig. 17(b-1,-2,-3)), it is recognized that the degree of distortion increases more and more with the increase of height reduction. In particular, the deformation of the coarse mesh such as Mesh-1 is unacceptable in the physical sense, because the mesh around point (E) in Fig. 16 interferes with the lower die. In other words, the initial coarse mesh cannot be compatible with the complex geometry of die when large deformation occurs. To overcome this di5culty, Liu et al. proposed the algorithm of automatic local mesh subdivision [29]. Although the degree of distortion diminishes as a whole for the ner mesh such as Mesh-2 or Mesh-3, the error remains large as shown in Fig. 20 because of the local severe distortion. Therefore, the remeshing technique is particularly necessary in the large deformation problem. To use this remeshing technique in the upsetting problem, we add the predicted relative error (∗ ) in the re nement criterion. That is, we let remeshing to take place when the predicted relative error reaches
318
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324
Fig. 16. Boundary and loading conditions of the upsetting problem.
∗max . This is because the remeshing may occur in every increment step if the maximum permissible error in an element ( e! i; max ) is only compared with the error in each element, such as in the elastic analysis. Mesh-1 is used as the initial mesh and two re nement criteria (∗max = 8%; q = 1 and ∗max = 8%; q = 1:5) are utilized. The meshes before and after remeshing and error distribution at the moment are shown in Figs. 18 and 19 for each re nement criterion. The number of elements generated in one remeshing step for q = 1:5 is larger than that for q = 1. The averages of mesh qualities for the two cases are shown in Table 3. Since the mesh quality is reduced as height reduction increases, we now study the mesh smoothing procedure to enhance the mesh quality at the boundary mesh. In Fig. 20, e! ’s for the above meshes are compared as height reduction increases. The errors of Mesh-1–3 increase monotonically but the errors of Mesh-r8-1 and Mesh-r8-1.5 drop at the instant of every remeshing. The forces required during the upsetting process are plotted in Fig. 21 . It is seen that the required force after remeshing decreases because the internal work decreases after the mapping of the state variables. The slope of the force curve tends to approach from Mesh-1 to Mesh-3 as the re nement progresses. If the total number of elements is not changed during the remeshing, the force curve after remeshing will follow the force curve of Mesh-1 [17]. From the above results, it can be concluded that this remeshing algorithm based on the BPM generates well-shaped ne meshes controlled by two parameters in the large deformation problem, also.
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324
319
Fig. 17. Initial and deformed shapes of uniform meshes in the upsetting problem: (a-1) Mesh-1 (bubble size = 4), (a-2) Mesh-2 (bubble size = 2), (a-3) Mesh-3 (bubble size = 1), (b-1) 84.88% reduction, (b-2) 81.76% reduction, (b-3) 84.02% reduction.
6. Concluding remarks In this paper, we have proposed a new remeshing algorithm using the BPM in two-dimensional nite element analysis and applied this algorithm to the large deformation problem. The optimal node position was determined by the BPM and the advancing front concept in Delaunay triangulation was adopted to construct the triangular meshes. To use the information of nodes and elements generated by the BPM in FEA, the eIective bandwidth minimization technique using the CCM was developed. In the remeshing procedure, a re ning circle was introduced to specify the region where ne and regular meshes are generated. The new bubble size to be located within the re ning circle is determined by the maximum permissible relative error ∗max and the bubble size control factor q. The remeshing=re nement algorithm using the BPM is simple and useful because the global remeshing is carried out automatically if only ∗max and q are assigned. The capability of this algorithm was presented in the linear and nonlinear problems through several re nement criteria. From the elastic examples, it is seen that ne meshes are generated at the point where strain gradients are high, and a high rate of convergence in the adaptive re nement is obtained. Above all, we are sure that the present algorithm is also helpful in nding the converged and suitable solutions in a large deformation process such as the upsetting problem. The development of a three-dimensional remeshing algorithm based on the BPM is now in progress.
320
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324
Fig. 18. Finite element meshes before and after remeshing (left: before remeshing, middle: error distribution, right: after remeshing, ∗max = 8%; q = 1): (a) rst remeshing (1.95% reduction, number of elements after remeshing: 185), (b) second remeshing (19.38% reduction, number of elements after remeshing: 303), (c) third remeshing (38.75% reduction, number of elements after remeshing: 398).
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324
321
Fig. 19. Finite element meshes before and after remeshing (left: before remeshing, middle: error distribution, right: after remeshing, ∗max =8%; q=1:5): (a) rst remeshing (1.95% reduction, number of elements after remeshing: 291), (b) second remeshing (29.77% reduction, number of elements after remeshing: 750), (c) third remeshing (50.76% reduction, number of elements after remeshing: 1072).
322
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324
Table 3 Average quality of meshes generated by remeshing in the upsetting problem Re nement criterion
First remeshing
Second remeshing
Third remeshing
∗max = 8%; q = 1 ∗max = 8%; q = 1:5
0.883 0.862
0.880 0.847
0.869 0.844
Fig. 20. L2 norm in strain error (e! ) vs. height reduction.
Fig. 21. Required force vs. height reduction.
Acknowledgements The authors gratefully acknowledge the nancial support of the Ministry of Science and Technology under the National Research Laboratory program (Grant Number 00-N-NL-01-C-026).
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324
323
References [1] C.K. Lee, R.E. Hobbs, Automatic adaptive nite element mesh generation over arbitrary two-dimensional domain using advancing front technique, Comput. & Structures 71 (1999) 9–34. [2] S.H. Lo, Delaunay triangulation of non-convex planar domains, Internat. J. Numer. Methods Eng. 28 (1989) 2695–2707. [3] T.D. Blacker, M.B. Stephenson, Paving: a new approach to automated quadrilateral mesh generation Internat. J. Numer. Methods Eng. 32 (1991) 811–847. [4] P.L. Baehmann, S.L. Wittchen, M.S. Shephard, K.R. Grice, M.A. Yerry, Robust, geometrically based, automatic two-dimensional mesh generation, Internat. J. Numer. Methods Eng. 24 (1987) 1043–1078. [5] D.O. Potyondy, P.A. Wawrzynek, A.R. IngraIea, An algorithm to generate quadrilateral or triangular element surface meshes in arbitrary domains with applications to crack propagation, Internat. J. Numer. Methods Eng. 38 (1995) 2677–2701. [6] S.B. Petersen, J.M.C. Rodrigues, P.A.F. Martins, Automatic generation of quadrilateral meshes for the nite element analysis of metal forming processes, Finite Elem. Anal. Des. 35 (2000) 157–168. [7] K. Shimada, D.C. Gossard, Automatic triangular mesh generation of trimmed parametric surfaces for nite element analysis, Comput. Aided Geom. Design 15 (1998) 199–222. [8] O.C. Zienkiewicz, J.Z. Zhu, A simple error estimator and adaptive procedure for practical engineering analysis, Internat. J. Numer. Methods Eng. 24 (1987) 337–357. [9] H.-S. Oh, R.C. Batra, Application of Zienkiewicz–Zhu’s error estimate with superconvergent patch recovery to hierarchical p-re nement, Finite Elem. Anal. Des. 31 (1999) 273–280. [10] O.C. Zienkiewicz, J.Z. Zhu, The superconvergent patch recovery and a posteriori error estimates, Part 1: the recovery technique, Internat. J. Numer. Methods Eng. 33 (1992) 1331–1364. [11] O.C. Zienkiewicz, J.Z. Zhu, The superconvergent patch recovery and a posteriori error estimates, Part 2: error estimates and adaptivity, Internat. J. Numer. Methods Eng. 33 (1992) 1365–1382. [12] B. Boroomand, O.C. Zienkiewicz, Recovery by equilibrium in patches, Internat. J. Numer. Methods Eng. 40 (1997) 137–164. [13] H.T.Y. Yang, M. Heinstein, J.M. Shih, Adaptive 2D nite element simulation of metal forming processes, Internat. J. Numer. Methods Eng. 28 (1989) 1409–1428. [14] A.R. Khoei, R.W. Lewis, Adaptive nite element remeshing in a large deformation analysis of metal powder forming, Internat. J. Numer. Methods Eng. 45 (1999) 801–820. [15] J.-H. Cheng, Automatic adaptive remeshing for nite element simulation of forming process, Internat. J. Numer. Methods Eng. 26 (1988) 1–18. [16] G.H. Paulino, I.F.M. Menezes, J.B. Cavalcante Neto, L.F. Martha, A methodology for adaptive nite element analysis: towards an integrated computational environment, Comput. Mech. 23 (1999) 361–388. [17] Ch. Pavana Chand, R. Krishna Kumar, Remeshing issues in the nite element analysis of metal forming problems, J. Mater. Process. Technol. 75 (1998) 63–74. [18] K. Shimada, N. Mori, T. Kondo, T. Itoh, K. Kase, A. Makinouchi, Automated mesh generation for nite element analysis of sheet metal forming, Internat. J. Vehicle Design 21 (2=3) (1999) 278–291. [19] J. Puttonen, Simple and eIective bandwidth reduction algorithm, Internat. J. Numer. Methods Eng. 19 (1983) 1139–1152. [20] S.J. Owen, S. Saigal, Surface mesh sizing control, Internat. J. Numer. Methods Eng. 47 (2000) 497–511. [21] K. Mulmuley, Computational Geometry: an Introduction through Randomized Algorithms, Prentice-Hall, Englewood CliIs, NJ, 1994, pp. 288–289. [22] C. Du, A note on nding nearest neighbours and constructing Delaunay triangulation in the plane, Comm. Numer. Methods Eng. 14 (1998) 871–877. [23] D.F. Watson, Computing the n-dimensional Delaunay tesselation with application to Voronoi polytopes, Comput. J. 24 (1981) 167–172. [24] S.J. Kim, S.W. Chung, Numerical simulation of nite rotation based on strain kinematics by polar decomposition, International Conference on Computational Engineering and Science, Atlanta, GA, October 6 –9, 1998, pp. 691– 696. [25] S.J. Kim, J.T. Oden, Finite element analysis of a class of problems in nite elastoplasticity based on the thermodynamical theory of materials of type N , Comput. Methods Appl. Mech. Eng. 53 (1985) 277–307.
324
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324
[26] S.J. Kim, J.H. Kim, Finite element analysis of laminated composites with contact constraint by extended interior penalty method, Internat. J. Numer. Methods Eng. 36 (1993) 3421–3439. [27] O.C. Zienkiewicz, J.Z. Zhu, N.G. Gong, EIective and practical h-p-version adaptive analysis procedure for the nite element method, Internat. J. Numer. Methods Eng. 28 (1989) 879–891. [28] D. Liu, Z.J. Luo, M.X. Gu, The algorithm of automatic local mesh subdivision and its application to nite-element analysis of a large deformation forming process, J. Mater. Process. Technol. 83 (1998) 164–169. [29] J.Z. Zhu, O.C. Zienkiewicz, E. Hinton, J. Wu, A new approach to the development of automatic quadrilateral mesh generation, Internat. J. Numer. Methods Eng. 32 (1991) 849–866.
A remeshing algorithm based on bubble packing method and its application to large deformation problems Soon Wan Chung, Seung Jo Kim ∗ Department of Aerospace Engineering, School of Mechanical and Aerospace Engineering, Seoul National University, San 56-1, Shinlim-dong, Kwanak-gu, Seoul 151-742, South Korea Received 15 April 2001; accepted 11 January 2002
Abstract The remeshing algorithm based on the bubble packing method (BPM) is developed for nite element analysis (FEA) because the optimal position of nodes can be obtained systematically when only the bubble size function is given. The triangular meshes are generated by Delaunay triangulation with the advancing front concept. To use the automatically generated mesh information in FEA, a new bandwidth minimization scheme with high e5ciency in CPU time is also developed. A re ning circle is proposed to specify the re nement area using the bubble size in the remeshing step. The remeshing algorithm is applied to the metal forming problem to replace the distorted mesh at the largely deformed area with regular and ne meshes. The numerical results show that the remeshing algorithm based on the BPM works well at the region with large error and is able to control the re nement area and the new mesh size easily through two parameters (∗max ; q). ? 2002 Elsevier Science B.V. All rights reserved. Keywords: Remeshing; Bubble packing method; Bandwidth minimization; Re ning circle; Large deformation
1. Introduction For several decades, adaptive remeshing by referencing the error estimation in nite element analysis (FEA) has been investigated actively, because mesh re nement is crucial in nonlinear problems, such as the metal forming process. Therefore, the remeshing algorithm based on bubble packing method (BPM), which is one of useful mesh generation algorithms, is developed and applied to large deformation problems in this paper. Generally, FEA including the remeshing is composed of the following three modules: automatic mesh generation, error estimation and adaptive remeshing. A number of automatic mesh generation ∗
Corresponding author. Tel.: +82-2-880-7388; fax: +82-2-887-2662. E-mail address: [email protected] (S.J. Kim).
0168-874X/02/$ - see front matter ? 2002 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 8 - 8 7 4 X ( 0 2 ) 0 0 0 7 5 - 6
302
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324
algorithms in the two-dimensional domain have been proposed such as the node connection approach [1–3] and the spatial decomposition procedure [4 – 6]. In the node connection approach, Delaunay triangulation or the advancing front method is generally used, and in the spatial decomposition method, the domain is subdivided recursively until the mesh size attains the desired level. Recently, the BPM [7] has been proposed, which has the advantages of the optimal node placement by the physically based relaxation and ease of mesh size control. In the eld of error estimation, the Zienkiewicz–Zhu error estimator [8,9] is popular due to simple implementation and cost eIectiveness. To estimate the discretization error, the recovered values are used in place of the exact solutions when the exact ones are not available in the practical problems. As the recovery technique, superconvergent patch recovery (SPR) [10,11] and recovery by equilibrium in patches (REP) [12] are representative. The adaptive remeshing is needed at the region where the solution gradient is high or when the initial mesh becomes distorted severely, for example, in metal forming problems. The remeshing criterion aIecting re nement area is generally based on error estimation [13–16] or the element distortion indicator [17]. In addition, the routine for accurate data transfer between an old mesh and a new mesh should be included. In this paper, the BPM is used in the remeshing algorithm [18] as well as mesh generation algorithm. As the BPM considers a node as the center of a bubble (a circle in two dimensions) and represents mesh size by bubble size, the re nement area and new mesh size in this area can be calculated from the bubble size function. We propose a re ning circle centered on the candidate element for re nement and calculate the new bubble size in the circle using control factor q. In the automatic mesh generation module, Delaunay triangulation with the advancing front concept is used to connect the nodes generated by the BPM. Also we adopt the Zienkiewicz–Zhu error estimator combined with SPR. As the numerical experiments, we rst investigate an elastostatic problem, the tension problem of a plate with an elliptical hole. Secondly, the metal forming process is simulated for the large deformation problem. The equilibrium equation based on the elasto-viscoplastic model and total Lagrangian formulation are used for the nonlinear nite element approximation. 2. Automatic mesh generation algorithm The automatic mesh generation in this study consists of two main parts: (1) optimal node placement through the BPM and (2) Delaunay triangulation with the advancing front concept. The algorithm of the BPM is based on Ref. [7] and is modi ed for simple implementation. The enhanced bandwidth minimization technique based on Ref. [19] is developed to reduce CPU time during the bandwidth minimization. 2.1. Bubble packing method This method determines the node placement for triangular elements in consideration of the equilibrium position by the attractive=repulsive interbubble forces of the nodes. The forces applied at the bubbles are calculated by a bubble size function. The close packing of bubbles is similar to a pattern of Voronoi polygons and therefore, Delaunay triangulation of these bubbles yields well-shaped triangles close to the regular triangles. Since the BPM compacts the bubbles directly inside the
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324
303
Fig. 1. Procedure of node distribution based on bubble packing method: (a) nodes given in input data, (b) boundary line bubble packing, (c) Surface bubble packing-1, (d) Surface bubble packing-2, (e) Bubble distribution after dynamic simulation.
domain according to the bubble size function, it is noted that BPM does not have to interpolate the background mesh [20,21] to de ne the mesh size. The mesh size in adaptive remeshing is determined by the concept of a re ning circle and will be described in Section 4. The BPM consists of the following steps: 1. 2. 3. 4.
preparation of input data for the boundary of the computational domain, boundary line packing of bubbles, surface packing of bubbles, physically based dynamic simulation.
Step 1. Preparation of input data for the boundary of the computational domain (Fig. 1(a)): In order to use in automatic mesh generation, the input data provided by a user should be concise. Hence, nodes at the vertex, lines connecting the nodes, and surfaces composed of lines are written in the input le for mesh generation. The bubble size of a node is used in the following packing procedure and the dynamic simulation. A xed node is often necessary when boundary condition must be imposed at a speci c point in FEA. For this case, we make its bubble size a negative value for distinction. Step 2. Boundary line packing of bubbles: The bubbles are packed along the boundary lines connecting nodes given in the input data (Fig. 1(b)). At this time, the bubbles are placed compactly according to the bubble size function.
304
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324
Fig. 2. Delaunay triangulation with the advancing front concept.
Step 3. Surface packing of bubbles: The center of the domain is given in the input data and then the line connecting the center and the vertex is packed in the same way as in Step 2. We call this line diagonal line for surface packing (DLSP). After each DLSP is packed, bubbles are inserted at the line between nodes on diIerent DLSPs (Fig. 1(c)). The nodes generated in a hollow domain such as a void are removed by the vertical ray shooting algorithm (Fig. 1(d)) [22]. Step 4. Physically based dynamic simulation: The above procedures (Steps 1–3) are for the initial bubble placement. The bubble distribution from the procedures may make gaps or overlaps between bubbles. To overcome this problem, bubble con guration balanced with the interbubble attractive=repulsive forces is found in this step (Fig. 1(e)). The force function is de ned as a cubic function of the distance between the centers of two bubbles. The forces applied at boundary bubbles are projected to the boundary line to constrain the bubbles on the boundary. The xed node speci ed at Step 1 applies the force to other bubbles, but it does not take any forces by other bubbles. The equations of motion of each bubble are integrated by the fourth-order Runge–Kutta method [7]. After integration, a bubble is deleted when two bubbles are too close for their bubble sizes and when the bubble exists outside the domain. 2.2. Delaunay triangulation with the advancing front concept Triangular meshes are generated by connecting nodes obtained from the BPM, and for this procedure Delaunay triangulation with the advancing front concept is incorporated. As the triangulation is started from boundary edges and advances into the interior domain until the whole domain is covered completely, we carry out Delaunay triangulation using the largest interior angle property of the Voronoi nearest neighbors [23] (Fig. 2). Newly generated edges are considered as new boundaries and the nearest node is searched by the largest interior theorem. This node connection algorithm can be as e5cient as the point insertion method [24], if nodes are organized properly [23]. Moreover, it is easy to extend the algorithm in the case that the problem domain is multiply connected because the generation of fronts from the inner boundaries is straightforward. 2.3. Additional algorithms for checking node position The nodes lying outside the domain in Step 3 or 4 of the BPM are removed using the vertical ray shooting algorithm. In the Delaunay triangulation stage, the node search algorithm (NSA) on boundary edge is incorporated keeping in mind that the candidate node should be located at only one side. If we assume that the direction vectors (vx ; vy ) of boundary lines are in a counterclockwise
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324
305
Fig. 3. Sketch for the node search algorithm.
direction as in Fig. 3, the admissible region of the outer (inner) boundary can be determined from ◦ the normal vector obtained by rotating the direction vector 90 counterclockwise. Since the normal vector of a boundary line is (−vy ; vx ), the line equation is as follows: − vy x + vx y + c = 0:
(1)
If a node (x1 ; y1 ) is on the left side of the direction vector, − vy x1 + vx y1 + c ¿ 0;
(2)
and if a node is on the right side of the direction vector, − vy x1 + vx y1 + c ¡ 0:
(3)
As the admissible region for the outer=inner boundary is on the left=right side, the node satisfying inequality (2)=(3) is allowable as the candidate node for the boundary lines. 2.4. Recursive=adaptive bubble population control To increase the regularity of a triangular mesh, adaptive and recursive bubble population controls are used. Adaptive control: The number of surrounding meshes of each interior node is calculated, and if the average value of the numbers is smaller than 5.9 or larger than 6.1, bubble addition or deletion is performed, because the number for regular meshes is 6. If the number of surrounding meshes of a node is smaller than 5, a new bubble is inserted at the midpoint of the edge with maximum interior angle (Fig. 4(a)). If the number is larger than 7, a node which is adjacent to a neighboring mesh with a smaller interior angle than the other neighboring mesh ( ¡ ) is deleted between two bubbles on the edge with minimum interior angle (Fig. 4(b)). Recursive control: When the adaptive bubble control is used, this algorithm returns to the dynamic simulation (Step 4 of BPM) to nd a new equilibrium position of the bubbles. This recursion is repeated until the adaptive control is not activated or the number of repetitions reaches 5. 2.5. Bandwidth minimization Node and element numberings from automatic mesh generation are distributed randomly, and the maximum bandwidth of global stiIness matrix may be very large. Therefore, a bandwidth minimization technique is needed if the band solver is used to solve the simultaneous equations resulted by FE approximation. The bandwidth minimization technique in this study is based on Puttonen’s
306
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324
(a)
(b)
Fig. 4. Adaptive bubble control: (a) bubble addition and (b) bubble deletion.
bandwidth reduction algorithm [19]. The basic idea is to compact the rows and columns of the symmetric node connectivity matrix. If the node connectivity matrix is a square matrix, the required memory in the problem with 10 000 nodes is 100 Mbytes, for example. However, if the condensed connectivity matrix (CCM) which contains only nonzero column numbers on every row is used for row=column change, the required memory is at most 100 kbytes. In Ref. [19], the use of the CCM is commented on, but a full description is not presented. Hence, in this study the algorithm using the CCM is explained and the enhancement of CPU time by this algorithm is shown in several examples. Since Puttonen’s algorithm uses row=column exchanges, row=column exchanges in the CCM can be classi ed into four cases according to the existence of nodes at the corresponding row=column. Fig. 5 shows four cases taking place in the CCM when the ith column and the jth column are exchanged in any kth row. In Cases 1 and 2, no changes appear in the CCM, and the changes in the CCM in other cases are as follows: Case 3: CCM (k; 1) = i; CCM (k; 2) = i + 1; CCM (k; 3) = i + 3 ⇒ CCM (k; 1) = i + 1; CCM (k; 2) = i + 3; CCM (k; 3) = j, Case 4: CCM (k; 1) = i + 1; CCM (k; 2) = i + 2; CCM (k; 3) = j ⇒ CCM (k; 1) = i; CCM (k; 2) = i + 1; CCM (k; 3) = i + 2. To show the e5ciency of the CCM, the CPU time is compared with that by using square connectivity matrix in three models (Fig. 6) through a PC with an Alpha-533 MHz CPU (Table 1). From the
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324 i
307
j
Case 1
Case 2
Case 3
Case 4
Fig. 5. Four cases of row=column exchange for the CCM.
Fig. 6. Models used for comparison of CPU time during bandwidth minimization.
Table 1 Comparison of CPU time
Model 1 Model 2 Model 3
Number of nodes
Number of elements
Square connectivity matrix (s)
CCM (s)
Change in half-bandwidth
173 399 593
304 716 1104
0.141 0.640 3.080
0.068 0.139 0.425
129 → 22 273 → 14 534 → 52
results, CPU time with the use of CCM is found to be much smaller than the case without the CCM, and the e5ciency is more prominent in the model with the larger number of nodes. As the node number is reordered through the bandwidth minimization procedure, the extra storage for relating the node numbers before and after bandwidth minimization is added. The Oowchart of the above automatic mesh generation algorithm is presented in Fig. 7.
308
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324 Input Minimal geometric data, Bubble size function, Constants(m,c)
Bubble Packing Method
Boundary line packing of bubbles
Surface packing of bubbles
Physically-based dynamic simulation
Delaunay triangulation with advancing front concept
Yes (Recursive bubble control)
Is adaptive bubble population control used? No Bandwidth minimization : CCM
Mesh generation completed!
Fig. 7. Flowchart of automatic mesh generation algorithm.
3. Nonlinear nite element analysis To use the remeshing algorithm based on the BPM in the large deformation problem, we perform the nonlinear nite element analysis. A kinematic strain measure [25] utilizing polar decomposition of the deformation gradient and constitutive equations based on the theory of material of type-N [26] are used to simulate geometrical and materially nonlinear behaviors. To realize contact conditions, the extended interior penalty method [27] is adopted. The resulting highly nonlinear equilibrium equation is reduced to the incremental weak form and approximated by the total Lagrangian formulation.
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324
309
3.1. Kinematics Assuming the unrotated current con guration Ct − and an intermediate con guration Cp− ; U can be decomposed additively as follows. [25]: U = Ue + Up − I;
(4)
where Ue and Up represent the elastic stretch and the plastic stretch tensors, respectively. The velocity gradient tensor and the velocity strain tensor take the following forms: ˙ t + R(U˙ e U−1 + U˙ p U−1 )Rt ; ˙ −1 = RR L = FF
(5)
D = L|sym = R(U˙ e U−1 |sym + U˙ p U−1 |sym )Rt :
(6)
From Eqs. (5) and (6), one can de ne the elastic and plastic strain rates as E˙ = U˙ e U−1 |sym ;
P˙ = U˙ p U−1 |sym :
(7)
3.2. Elasto-viscoplastic material model The theory of materials of type N consists of ve eld equations (conservation of mass, balance of linear and angular momenta, conservation of energy, law of entropy production) and the following two potentials: • Free energy functional: 1 1 1 2 2 {(tr E) + 2(tr E )} + h1 z + (h1 − h0 )exp(−mz) : = 0 2 m • Flow potential functional: = D0
∞ i=0
(−1)
i+1
i h2i−1 i!(2i − 1)J2i−1=2
1 where = n 3
The nal constitutive equations can be obtained as follows: S = [(tr E)I + 2E]; 0 h = [h1 + (h0 − h1 )exp(−mz)]; 0
h2 1 ˙ √ P = D0 S ; exp − J 2h J2 2 1 ˙ z˙ = S : P: h
n+1 n
(8)
:
(9)
(10) (11) (12) (13)
To integrate the above constitutive equations (12), (13) with respect to time, the predictor–corrector method for a hardening variable (z) and the backward Euler method for plastic strain (P) are used.
310
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324
4. Adaptive remeshing algorithm Whenever FEM is used to solve a diIerential equation, discretization error is introduced and the distorted mesh due to the large deformation process, in particular, may induce large approximation errors. In order to decrease this error, the remeshing procedure is performed over the whole domain and the re ned meshes are generated at the element with large error in the present study. For the re nement criterion, the maximum permissible error and the bubble size control factor are used. The mesh re nement can be implemented easily by reducing the bubble size in the BPM because the mesh size is determined by bubble size. Since node regeneration by the BPM is started at the boundary and the nodes may be added at the boundary lines during the re nement stage, the update of boundary information is important. The extra storage for boundary nodes commented in Section 2.5 is needed in this re nement stage of boundary lines in order to identify the boundary nodes among the nodes that resulted from the bandwidth minimization technique. After the boundary update is nished, the mesh generation procedure described in Section 2 is performed again. The Oowchart of this remeshing algorithm is shown in Fig. 8. 4.1. Recovery technique and error estimation To estimate a posterior error, the nodal point value is computed from the Gauss point value by the recovery technique and the recovered value obtained by interpolating the nodal point value is
Finite Element Analysis at (i) - th increment step
Error estimation
Refinement criterion
i=i+1
Is remeshing necessary?
Yes Automatic mesh generation
Mapping of state variables (stress, strain, displacement, et al)
Fig. 8. Flowchart of remeshing algorithm.
No i=i+1
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324
311
Node determined by recovery procedure Sampling point selected firstly Sampling point selected secondly(either of)
Fig. 9. Recovery at boundary nodes.
compared with the nite element solution. As the recovery technique, SPR [10] is adopted in this research. For the linear triangular element, the linear polynomial of strain is used as shown below: !1∗ = Pa;
(14) T
where P = [1; x; y] and a = [a1 ; a2 ; a3 ] . The unknown parameters a of Eq. (14) are calculated by a least-squares t of the set of superconvergent points. After minimizing the diIerence between the discontinuous strain (!) ˆ derived by nite element analysis and the linear expansion !1∗ , we can obtain the following matrix form: (15) a = A−1 b; n n T T where A = i=1 P (xi ; yi )P(xi ; yi ) and b = i=1 P (xi ; yi )!(x ˆ i ; yi ) and n is the number of sampling points. As the number of unknowns is 3, the number of sampling points must be 3 or more. To prevent the abrupt increase of recovered values at the boundary, we select the number of sampling points to be more than 3. If the number is less than 4 at the boundary nodes, the neighboring elements connected with any two nodes on the element patch are added (Fig. 9). After a is known, the nodal parameter (!S∗ ) is found out by substituting the coordinates of the node in Eq. (14), and the smoothed strain eld is interpolated by the nodal parameter. !∗ = N !S∗ :
(16)
Here, the interpolation function N is the same as the one used for the interpolation of displacements in FEM. The strain error using the recovered value is de ned by e! ≡ !∗ − !:ˆ The L2 norm of strain error e! and the predicted relative error can be written as m m
e! 2 ∗ e! = e! 2i = (!∗ − !) ˆ 2i d%; = ; !∗ 2 + e! 2 %i i=1 i=1
(17)
(18)
312
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324
Refining circle d1 2
d1+d2+d3 3
d3 2 d2 2
Candidate element for refinement Fig. 10. Re ning circle de ned by bubble size.
where m is the total number of elements. If a maximum permissible relative error (∗max ) is provided, the maximum permissible error in any element i can be de ned as shown [16]:
1=2 ∗ e! i; max = √max !∗ 2 + e! 2 : (19) m 4.2. Re9nement criterion In this research, since a re nement is provided at the element with error more than the maximum permissible error in an element ( e! i; max ), only the elements with high error are re ned as ∗max of Eq. (19) is larger. In other words, if ∗max is large, the re ned region is narrow and the increase in number of elements after re nement is also small. Therefore, ∗max can be chosen according to the error distribution of the example, and we will compare the re ned meshes according to the variation of ∗max in the linear elastostatic problem. After an element with large error is found, the re nement region and the new mesh size within the region should be systematically speci ed. A re ning circle is introduced in this algorithm to indicate the re nement region centered on the candidate element for re nement (Fig. 10). The center corresponds to the centroid of the candidate element for re nement and the radius is taken as the average bubble size of three vertices. The bubble size of a new node, which will be placed within this re ning circle, is calculated as follows: e! i; max q dnew = dave : (20) e! i Here, dave is the average bubble size of three vertices, and q is the factor to control the new bubble size dnew . When q ¿ 1, the reduction ratio of bubble size (dave −dnew )=dave becomes large, and when q ¡ 1, the ratio is relatively small. That is to say, the re nement region is determined by ∗max and the bubble size generated in the re ning circle can be controlled by q. Therefore, the new bubble
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324
313
size distribution by two parameters (∗max ; q) in addition to the updated boundary information are essential to this remeshing algorithm using BPM. 4.3. Data transfer When the remeshing stage is performed, the data transfer should be treated carefully. In the linear elastic problem, the transfer of a boundary node, at which boundary conditions or loading conditions are imposed, is only necessary because the loading is applied again at the undeformed con guration. The information at the boundary node is particularly important in this algorithm, because the automatic mesh generation is started at the boundary nodes and new nodes may be inserted at the boundary lines for two cases in the remeshing stage. One case is where a new node is inserted at the middle of the boundary line when the candidate element for re nement is adjacent to the boundary line. The other case is where new nodes are inserted during Step 2 of the BPM (boundary line packing of bubbles), when the distance between two nodes is greater than the sum of two bubble diameters due to the large deformation. If the boundary conditions of two end nodes on the boundary line agree, the same condition is applied at the new node. In the remeshing stage of a nonlinear problem, not only the information at the boundary node but also the state variable such as displacement, stress and strain should be considered. A Gauss point value such as stress or strain is transferred by the following three steps, and a nodal point value such as displacement is transferred by only Step 2 [15]. (1) Transformation of Gauss point value into nodal point value in old mesh by SPR, (2) Mapping between nodal point values in old and new meshes, (3) Transformation of nodal point value into Gauss point value in new mesh. In step 2, we rst use NSA in order to search an old element to which a new node belongs. When we consider three sides of an old element as outer boundaries, a new node exists in the element if the node is located in the admissible domain, that is to say, the node satis es Eq. (2) for all three lines. The new nodal point value is calculated by interpolating the old nodal point value using the area coordinates. In step 3, the Gauss point value is the average of three nodal point values because one point Gauss quadrature is utilized. In Ref. [17], it is addressed that these steps 2 and 3 procedures are superior to the direct transfer when the nodal value of an old mesh is transferred into the Gauss point of a new mesh. The set of contact candidate nodes is also considered in the upsetting problem when the information at the boundary node is transferred. If a new node added at the boundary line is between the contact candidate nodes, the node is included in the set of contact candidate nodes, and the contact force applied at the node is calculated.
5. Numerical experiments 5.1. Linear elastic analysis We consider a linear elastostatic problem, the tension of a square plate with an elliptical hole under the plane strain assumption. The geometry of the model is shown in Fig. 11. Young’s modulus and
314
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324
Fig. 11. Geometric description of square plate with elliptical hole.
Fig. 12. Meshes generated by uniform re nement in the tension problem: (a) mesh-1 (bubble size = 0:4), (b) mesh-2 (bubble size = 0:2), (c) mesh-3 (bubble size = 0:1).
Poisson’s ratio are E = 105 ; ( = 0:3, respectively, and the error in the equivalent strain (!eq = 2 ! : !ij ) is used for the error estimation. 3 ij At rst, three meshes are generated with uniform bubble sizes: 0.4, 0.2 and 0.1, respectively. The resulting meshes are regular as shown in Figs. 12(a) – (c) and the numbers of elements are 250, 1031 and 4264, respectively. The error ( e! i ) distributions of the three meshes are presented in Fig. 13, and the error decreases as the element size decreases. Since the error at point (A) in Fig. 11 is larger than any other region, it can be expected that more re ned meshes will be necessary at point (A). The converged equivalent strain is obtained by the extrapolation based on the following equation using the strains for the three uniform bubble sizes [28]. e 6 Chmin(p; ) ≈ CN −0:5 min(p; ) :
(21)
Here, h is the size of the element, N the number of degrees of freedom, C a positive constant depending on the problem, p the polynomial degree of shape function, and the strength of the
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324
315
Fig. 13. Error (e! i ) distribution of uniformly re ned meshes in the tension problem: (a) mesh-1, (b) mesh-2, (c) mesh-3.
singularities. Finally, the relative error in the L2 norm is calculated by 1=2 ! 2 − ! e! ˆ 2 = = : ! ! 2
(22)
For the adaptive re nement, Mesh-1 is used as the initial mesh and remeshing is continued until the relative error () becomes less than 5%. We take account of ve ∗max ’s (10%, 9%, 8%, 7%, 6%) for q = 1. The re ned meshes are generated around point (A) as shown in Fig. 14, and the mesh size is reduced as the re nement progresses. When ∗max is small, the number of elements generated during each remeshing step increases. The mesh quality de ned by the following equation is calculated and the averages of mesh qualities are summarized in Table 2. r (23) Q=2 ; R where r is the radius of the inscribed circle and R the radius of the circumscribed circle. The comparison of the convergence rates in uniform and adaptive re nements is shown in Fig. 15. Obviously, the rate of convergence in adaptive re nement is higher than that in uniform re nement,
316
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324
Fig. 14. Meshes generated by adaptive re nement in the tension problem (Mesh-r{∗max }%-{remeshing order} (number of elements after remeshing), q = 1): (a) Mesh-r10-1 (348), (b) Mesh-r10-2 (499), (c) Mesh-r10-3 (664), (d) Mesh-r10-4 (776), (e) Mesh-r8-1 (437), (f) Mesh-r8-2 (734), (g) Mesh-r8-3 (991), (h) Mesh-r6-1 (683), (i) Mesh-r6-2 (1299).
Table 2 Average quality of meshes generated by remeshing in the tension problem Re nement criterion
First remeshing
Second remeshing
Third remeshing
Fourth remeshing
Fifth remeshing
∗max = 10%; q = 1 ∗max = 9%; q = 1 ∗max = 8%; q = 1 ∗max = 7%; q = 1 ∗max = 6%; q = 1
0.912 0.912 0.900 0.903 0.910
0.910 0.891 0.903 0.909 0.911
0.907 0.908 0.905 0.887
0.911 0.897 0.918
0.900
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324
317
Relative error (η) (%)
30
10
Refinement type [Uniform and * , q)] Adaptive(ηmax Uniform 10%, 1 9%, 1 8%, 1 7%, 1 6%, 1
1 100
1000
5000
Number of degree of freedom
Fig. 15. Convergence rates of relative error(%) in the tension problem.
and the tendency is consistent with respect to ∗max . From these results, it can be said that this remeshing algorithm works well at the region with high error and the remeshing=re nement procedure using the BPM is eIective in controlling the re ning area and the mesh size. 5.2. Nonlinear elasto-viscoplastic analysis As a representative nonlinear problem, we consider the metal forming process, that is, the upsetting problem. The numerical model is an Al 2024 axisymmetric billet consisting of 16 mm radius and 40 mm height, and the boundary and loading conditions are shown in Fig. 16. If the solution does not converge within the given iteration number in an increment step, we halve the increment of die movement to search the equilibrium position under a smaller loading increment. The error in equivalent strain is used for the error calculation. For the uniform re nement, three meshes are generated with the uniform bubble sizes 4, 2 and 1 mm in the BPM (Fig. 17(a-1,-2,-3)) and the resulting numbers of elements are 100, 384 and 1512, respectively. From the deformed meshes of Mesh-1–3 (Fig. 17(b-1,-2,-3)), it is recognized that the degree of distortion increases more and more with the increase of height reduction. In particular, the deformation of the coarse mesh such as Mesh-1 is unacceptable in the physical sense, because the mesh around point (E) in Fig. 16 interferes with the lower die. In other words, the initial coarse mesh cannot be compatible with the complex geometry of die when large deformation occurs. To overcome this di5culty, Liu et al. proposed the algorithm of automatic local mesh subdivision [29]. Although the degree of distortion diminishes as a whole for the ner mesh such as Mesh-2 or Mesh-3, the error remains large as shown in Fig. 20 because of the local severe distortion. Therefore, the remeshing technique is particularly necessary in the large deformation problem. To use this remeshing technique in the upsetting problem, we add the predicted relative error (∗ ) in the re nement criterion. That is, we let remeshing to take place when the predicted relative error reaches
318
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324
Fig. 16. Boundary and loading conditions of the upsetting problem.
∗max . This is because the remeshing may occur in every increment step if the maximum permissible error in an element ( e! i; max ) is only compared with the error in each element, such as in the elastic analysis. Mesh-1 is used as the initial mesh and two re nement criteria (∗max = 8%; q = 1 and ∗max = 8%; q = 1:5) are utilized. The meshes before and after remeshing and error distribution at the moment are shown in Figs. 18 and 19 for each re nement criterion. The number of elements generated in one remeshing step for q = 1:5 is larger than that for q = 1. The averages of mesh qualities for the two cases are shown in Table 3. Since the mesh quality is reduced as height reduction increases, we now study the mesh smoothing procedure to enhance the mesh quality at the boundary mesh. In Fig. 20, e! ’s for the above meshes are compared as height reduction increases. The errors of Mesh-1–3 increase monotonically but the errors of Mesh-r8-1 and Mesh-r8-1.5 drop at the instant of every remeshing. The forces required during the upsetting process are plotted in Fig. 21 . It is seen that the required force after remeshing decreases because the internal work decreases after the mapping of the state variables. The slope of the force curve tends to approach from Mesh-1 to Mesh-3 as the re nement progresses. If the total number of elements is not changed during the remeshing, the force curve after remeshing will follow the force curve of Mesh-1 [17]. From the above results, it can be concluded that this remeshing algorithm based on the BPM generates well-shaped ne meshes controlled by two parameters in the large deformation problem, also.
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324
319
Fig. 17. Initial and deformed shapes of uniform meshes in the upsetting problem: (a-1) Mesh-1 (bubble size = 4), (a-2) Mesh-2 (bubble size = 2), (a-3) Mesh-3 (bubble size = 1), (b-1) 84.88% reduction, (b-2) 81.76% reduction, (b-3) 84.02% reduction.
6. Concluding remarks In this paper, we have proposed a new remeshing algorithm using the BPM in two-dimensional nite element analysis and applied this algorithm to the large deformation problem. The optimal node position was determined by the BPM and the advancing front concept in Delaunay triangulation was adopted to construct the triangular meshes. To use the information of nodes and elements generated by the BPM in FEA, the eIective bandwidth minimization technique using the CCM was developed. In the remeshing procedure, a re ning circle was introduced to specify the region where ne and regular meshes are generated. The new bubble size to be located within the re ning circle is determined by the maximum permissible relative error ∗max and the bubble size control factor q. The remeshing=re nement algorithm using the BPM is simple and useful because the global remeshing is carried out automatically if only ∗max and q are assigned. The capability of this algorithm was presented in the linear and nonlinear problems through several re nement criteria. From the elastic examples, it is seen that ne meshes are generated at the point where strain gradients are high, and a high rate of convergence in the adaptive re nement is obtained. Above all, we are sure that the present algorithm is also helpful in nding the converged and suitable solutions in a large deformation process such as the upsetting problem. The development of a three-dimensional remeshing algorithm based on the BPM is now in progress.
320
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324
Fig. 18. Finite element meshes before and after remeshing (left: before remeshing, middle: error distribution, right: after remeshing, ∗max = 8%; q = 1): (a) rst remeshing (1.95% reduction, number of elements after remeshing: 185), (b) second remeshing (19.38% reduction, number of elements after remeshing: 303), (c) third remeshing (38.75% reduction, number of elements after remeshing: 398).
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324
321
Fig. 19. Finite element meshes before and after remeshing (left: before remeshing, middle: error distribution, right: after remeshing, ∗max =8%; q=1:5): (a) rst remeshing (1.95% reduction, number of elements after remeshing: 291), (b) second remeshing (29.77% reduction, number of elements after remeshing: 750), (c) third remeshing (50.76% reduction, number of elements after remeshing: 1072).
322
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324
Table 3 Average quality of meshes generated by remeshing in the upsetting problem Re nement criterion
First remeshing
Second remeshing
Third remeshing
∗max = 8%; q = 1 ∗max = 8%; q = 1:5
0.883 0.862
0.880 0.847
0.869 0.844
Fig. 20. L2 norm in strain error (e! ) vs. height reduction.
Fig. 21. Required force vs. height reduction.
Acknowledgements The authors gratefully acknowledge the nancial support of the Ministry of Science and Technology under the National Research Laboratory program (Grant Number 00-N-NL-01-C-026).
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324
323
References [1] C.K. Lee, R.E. Hobbs, Automatic adaptive nite element mesh generation over arbitrary two-dimensional domain using advancing front technique, Comput. & Structures 71 (1999) 9–34. [2] S.H. Lo, Delaunay triangulation of non-convex planar domains, Internat. J. Numer. Methods Eng. 28 (1989) 2695–2707. [3] T.D. Blacker, M.B. Stephenson, Paving: a new approach to automated quadrilateral mesh generation Internat. J. Numer. Methods Eng. 32 (1991) 811–847. [4] P.L. Baehmann, S.L. Wittchen, M.S. Shephard, K.R. Grice, M.A. Yerry, Robust, geometrically based, automatic two-dimensional mesh generation, Internat. J. Numer. Methods Eng. 24 (1987) 1043–1078. [5] D.O. Potyondy, P.A. Wawrzynek, A.R. IngraIea, An algorithm to generate quadrilateral or triangular element surface meshes in arbitrary domains with applications to crack propagation, Internat. J. Numer. Methods Eng. 38 (1995) 2677–2701. [6] S.B. Petersen, J.M.C. Rodrigues, P.A.F. Martins, Automatic generation of quadrilateral meshes for the nite element analysis of metal forming processes, Finite Elem. Anal. Des. 35 (2000) 157–168. [7] K. Shimada, D.C. Gossard, Automatic triangular mesh generation of trimmed parametric surfaces for nite element analysis, Comput. Aided Geom. Design 15 (1998) 199–222. [8] O.C. Zienkiewicz, J.Z. Zhu, A simple error estimator and adaptive procedure for practical engineering analysis, Internat. J. Numer. Methods Eng. 24 (1987) 337–357. [9] H.-S. Oh, R.C. Batra, Application of Zienkiewicz–Zhu’s error estimate with superconvergent patch recovery to hierarchical p-re nement, Finite Elem. Anal. Des. 31 (1999) 273–280. [10] O.C. Zienkiewicz, J.Z. Zhu, The superconvergent patch recovery and a posteriori error estimates, Part 1: the recovery technique, Internat. J. Numer. Methods Eng. 33 (1992) 1331–1364. [11] O.C. Zienkiewicz, J.Z. Zhu, The superconvergent patch recovery and a posteriori error estimates, Part 2: error estimates and adaptivity, Internat. J. Numer. Methods Eng. 33 (1992) 1365–1382. [12] B. Boroomand, O.C. Zienkiewicz, Recovery by equilibrium in patches, Internat. J. Numer. Methods Eng. 40 (1997) 137–164. [13] H.T.Y. Yang, M. Heinstein, J.M. Shih, Adaptive 2D nite element simulation of metal forming processes, Internat. J. Numer. Methods Eng. 28 (1989) 1409–1428. [14] A.R. Khoei, R.W. Lewis, Adaptive nite element remeshing in a large deformation analysis of metal powder forming, Internat. J. Numer. Methods Eng. 45 (1999) 801–820. [15] J.-H. Cheng, Automatic adaptive remeshing for nite element simulation of forming process, Internat. J. Numer. Methods Eng. 26 (1988) 1–18. [16] G.H. Paulino, I.F.M. Menezes, J.B. Cavalcante Neto, L.F. Martha, A methodology for adaptive nite element analysis: towards an integrated computational environment, Comput. Mech. 23 (1999) 361–388. [17] Ch. Pavana Chand, R. Krishna Kumar, Remeshing issues in the nite element analysis of metal forming problems, J. Mater. Process. Technol. 75 (1998) 63–74. [18] K. Shimada, N. Mori, T. Kondo, T. Itoh, K. Kase, A. Makinouchi, Automated mesh generation for nite element analysis of sheet metal forming, Internat. J. Vehicle Design 21 (2=3) (1999) 278–291. [19] J. Puttonen, Simple and eIective bandwidth reduction algorithm, Internat. J. Numer. Methods Eng. 19 (1983) 1139–1152. [20] S.J. Owen, S. Saigal, Surface mesh sizing control, Internat. J. Numer. Methods Eng. 47 (2000) 497–511. [21] K. Mulmuley, Computational Geometry: an Introduction through Randomized Algorithms, Prentice-Hall, Englewood CliIs, NJ, 1994, pp. 288–289. [22] C. Du, A note on nding nearest neighbours and constructing Delaunay triangulation in the plane, Comm. Numer. Methods Eng. 14 (1998) 871–877. [23] D.F. Watson, Computing the n-dimensional Delaunay tesselation with application to Voronoi polytopes, Comput. J. 24 (1981) 167–172. [24] S.J. Kim, S.W. Chung, Numerical simulation of nite rotation based on strain kinematics by polar decomposition, International Conference on Computational Engineering and Science, Atlanta, GA, October 6 –9, 1998, pp. 691– 696. [25] S.J. Kim, J.T. Oden, Finite element analysis of a class of problems in nite elastoplasticity based on the thermodynamical theory of materials of type N , Comput. Methods Appl. Mech. Eng. 53 (1985) 277–307.
324
S.W. Chung, S.J. Kim / Finite Elements in Analysis and Design 39 (2003) 301 – 324
[26] S.J. Kim, J.H. Kim, Finite element analysis of laminated composites with contact constraint by extended interior penalty method, Internat. J. Numer. Methods Eng. 36 (1993) 3421–3439. [27] O.C. Zienkiewicz, J.Z. Zhu, N.G. Gong, EIective and practical h-p-version adaptive analysis procedure for the nite element method, Internat. J. Numer. Methods Eng. 28 (1989) 879–891. [28] D. Liu, Z.J. Luo, M.X. Gu, The algorithm of automatic local mesh subdivision and its application to nite-element analysis of a large deformation forming process, J. Mater. Process. Technol. 83 (1998) 164–169. [29] J.Z. Zhu, O.C. Zienkiewicz, E. Hinton, J. Wu, A new approach to the development of automatic quadrilateral mesh generation, Internat. J. Numer. Methods Eng. 32 (1991) 849–866.