- Email: [email protected]
A novel meshing algorithm for dynamic finite element analysis
Precision Engineering 27 (2003) 245–257
A novel meshing algorithm for dynamic finite element analysis Lihui Wang a,∗ , Toshimichi Moriwaki b a
Integrated Manufacturing Technologies Institute, National Research Council of Canada, 800 Collip Circle, London, Ont., Canada N6G 4X8 b Department of Mechanical Engineering, Kobe University, 1-1 Rokkodai-cho, Nada-ku, Kobe 657-8501, Japan Received 12 July 2002; received in revised form 31 October 2002; accepted 11 December 2002
Abstract This paper describes a new algorithm to handle problems in dynamic finite element analysis and run-time simulation, where mesh re-generation or dynamic adjustment is required. Based on a concept called coded box cell (CBC) substitution, this algorithm can be applied to both initial mesh generation and dynamic mesh adjustment along the border zones of multiple primitives that form an entire model. During the initial mesh generation, appropriate labels are assigned to the nodes and the faces of each finite element. These labels are used to facilitate decision-making in dynamic mesh adjustment. A mapping technique is adopted to transform curved surfaces to plain surfaces for the ease of automatic mesh adjustment while still using the same algorithm. The results of a case study show that a finite element mesh can be adjusted dynamically and locally around its border zone; and the algorithm can be utilized effectively to simulate the thermal behavior of a device under real operating conditions. Crown Copyright © 2003 Published by Elsevier Science Inc. All rights reserved. Keywords: Meshing; CAE; Finite element analysis; Mapping and inverse mapping; Thermal simulation
1. Introduction Today, simulation of production processes is becoming much more important in manufacturing, where complex machining processes and material handling are encountered. A particular emphasis has been given to research and development of virtual manufacturing having various computer-aided engineering (CAE) packages for analysis and simulation of machine behavior and dynamics, aiming at realizing an optimized production, right the first time, with high quality, low cost, and short lead time. It often requires an advanced system capability to analyze and simulate the dynamic behaviors of products, machines, or production lines, as if they are under real operating conditions. In the last two decades, a large number of CAE systems have been developed and put into practical use in product development, such as analyses of stress, strain, vibration, and thermal behavior of the product by applying finite element method (FEM) and boundary element method (BEM). Recently, there is a growing interest in the development of meshless methods for numerical solutions of partial differential equations. Some typical meshless methods include the reproducing kernel particle method (RKPM) [1], element ∗ Corresponding
author. Tel.: +1-519-430-7084; fax: +1-519-430-7064. E-mail address: [email protected] (L. Wang).
free Galerkin method (EFG) [2], finite point method (FPM) [3], and others [4–8]. However, considering the accuracy and stability of the analytical results, FEM remains as the most popular and practical approach available today for solving various complex engineering problems. Unlike the meshless methods, FEM or BEM inherently requires a properly meshed geometry, before any analysis can take place. Iterative analysis usually requires the mesh being partially refined or even regenerated to increase analytical accuracy. Despite a large number of meshing methods developed in the last three decades, algorithms for dynamic mesh generation and localized mesh adjustment are still missing from the literature. It is urgently required to establish a systematic approach that not only can handle the re-meshing problem effectively, but also can simulate a machine’s behavior dynamically and accurately, especially under real operating conditions. A complete understanding of a machine’s dynamics or its real life behavior is vital to both the structural design and the run-time machine control. The objective of this research is to develop an integrated design/analysis system that is capable of dynamic simulation of machines in real life. This paper introduces a new algorithm that is able to handle problems in FEM analysis and simulation, where mesh re-generation or dynamic adjustment is required. Based on a concept called coded box cell (CBC) substitution [9], this algorithm can be used to adjust those finite elements whose nodes disagree with each other during
0141-6359/03/$ – see front matter. Crown Copyright © 2003 Published by Elsevier Science Inc. All rights reserved. doi:10.1016/S0141-6359(03)00005-9
246
L. Wang, T. Moriwaki / Precision Engineering 27 (2003) 245–257
Nomenclature M F D
movable object fixed object direction of normal vector on the contacting plane of the objects F and M P direction of motion of the object M against F V direction perpendicular to both the directions P and D N(φ) number of nodes in the direction P or V in the FEM mesh of object φ δF , δM mesh sizes (lengths between adjacent nodes) of the objects F and M, respectively VF , VM sets of nodes of the objects F and M along their border zone, respectively LF , LM node labels of the objects F and M, respectively SF , SM face labels of the objects F and M, respectively IDCBC set of identification numbers of CBCs to be used for substitution R3 real object space RP3 projective space (x, y, z) coordinates in the real object space R3 (ξ, η, ζ) coordinates in the projective space RP3 f mapping function g inverse mapping function Mid mapping function identity pr instance of primitives Ipr image of pr
analysis at their contacting surfaces (border zones), in addition to the initial mesh generation. Following a brief review and classification of the available meshing methods, the processing procedures of the algorithm are described in detail. A mapping technique is adopted to transform curved surfaces to plain ones for the ease of automatic mesh adjustment by applying the same algorithm. The algorithm is verified through an example. The results of the case study show that an FEM mesh can be adjusted dynamically and locally around its border zone; and the algorithm can be utilized effectively to simulate its thermal behavior under real operating conditions.
2. A brief review Most of mesh generation tasks involve a great deal of effort to generate numerical data, such as node numbers, nodal point coordinates, parameter values, and element connectivity. This work could be very time-consuming and error-prone if done manually, especially when a large number of elements or very complex geometry are encountered. There are many bottlenecks in the meshing process: geometry cleanup, geometry simplification or de-featuring, and finite element mesh
generation. During the last few decades, much research have been carried out to automate the FEM mesh generation process [10], based on the geometric models created by using CAD systems. Initial efforts in the literature were focused on the twodimensional mesh generation algorithms based on triangulation [11–13] and quadrilateral elements [13,14]. Attention then shifted to the meshing algorithms for three-dimensional geometry [15–18] and geometry with curved surfaces [18–22] using tetrahedra, hexahedra, and other polyhedra like triangular prisms. Despite those accomplishments in the last few decades, research on new meshing algorithms remains active. Owen et al. [23,24] use an algorithm called H-Morph to generate a hexahedral-dominant finite element mesh for arbitrary volumes. As an extension of the Q-Morph algorithm [25], the H-Morph method starts with an initial tetrahedral mesh and systematically transforms and combines tetrahedra into hexahedra, without the need of decomposing an arbitrary volume. On the contrary, Lu et al. [26] use a feature-based approach for volume decomposition before applying an appropriate and well-established meshing algorithm to each of the meshable pieces. Li et al. [27] presented a unified scheme for simultaneously refining and coarsening a mesh. This adaptive meshing algorithm is suitable for a domain that has a moving boundary. Recently, notable efforts have been given to the automated hexahedral mesh generation, because hexahedral mesh generally provides higher “artificial stiffness” and better solutions than the tetrahedral one for many physical problems [28]. Sheffer and Bercovier [29] introduced a hexahedral meshing algorithm based on the embedded Voronoi graph (EVG), where the EVG is used to decompose a non-linear volume into simple sub-volumes meshable by basic meshing techniques. Lai et al. [30] proposed an enhanced sweeping mesh generation algorithm by creating all hexahedral element meshes between multiple source surfaces and multiple target surfaces, while most other sweeping techniques [31–33] create all hexahedral element meshes by projecting an existing single-surface mesh along a specified trajectory to a specified single target surface. Dhondt presented a method that can generate an unstructured 20-node brick element mesh for arbitrary structure, based on a triangulation of the structure’s surface [34,35]. Starting from a hexahedral master mesh encompassing the structure, the elements that are intersected by the triangulation are determined, cut, and re-meshed according to their cutting topologies. Whisker Weaving is another advancing front algorithm for all-hexahedral mesh generation [36,37]. It uses global information derived from grouping the mesh dual into surfaces to construct the connectivity of the mesh, and then positions the nodes afterwards. Based on the previous literature surveys [28,33,38] and authors’ investigation, major meshing algorithms available in public domain can be classified and summarized, as shown in Fig. 1, in terms of the techniques used and the topologies applied. Most meshing approaches can be applied to 3D cases.
L. Wang, T. Moriwaki / Precision Engineering 27 (2003) 245–257
247
Fig. 1. Classification of mesh generation methods.
3. Concept of dynamic mesh generation When a milling machine is boring an engine block, the relative positions among its components (table, column, and spindle, etc.) vary with time. It is generally difficult to analyze and simulate the dynamic behavior of the milling machine using conventional FEA tools while its components are moving relatively, because most of the FEA tools require the mesh to be regenerated once the corresponding geometry model is changed. It becomes impossible if the geometry model keeps changing rapidly due to the relative movements between its components. In contrast, dynamic mesh generation does not require the entire mesh to be regenerated, and it considers mesh generation in a constructive way when the geometry model of an assembly having movable components is designed. The fundamental process of the concept is shown in Fig. 2. In the first step, a set of primitives (basic geometric shapes of a component or an assembly) can be generated interactively by using a solid modeling system, followed by the initial mesh generation for each primitive. Hexahedral elements are selected in this approach. When the primitives are combined in the third step, the elements near their border zone are adjusted or modified locally in the fourth step, so that the mesh
of one primitive connects with that of another smoothly across the border. This operation needs to be performed only once if the geometry model consists of non-movable primitives only. Otherwise, as shown in the fifth step, the mesh along the border of the movable primitives has to be adjusted on the fly whenever its relative position is changed. This dynamic mesh modification is limited to those finite elements whose nodes disagree with each other within the border zone and the elements being affected. It is not difficult to adjust the mesh along the border zone if the sizes of mesh elements (across the border) of two primitives are the same. The mesh can be modified, in this case, by adjusting the nodal coordinates of those conflicting elements in the border zone. However, in most cases, the mesh elements of one primitive must be subdivided in order that the numbers of the elements of both primitives agree with each other within the border zone. A requirement to be met in subdividing an initial mesh is to localize the modification only in the area near the border zone so as to increase the efficiency of dynamic adjustment. A set of pre-cut hexahedral cells, called CBC as shown in Fig. 3, are defined and used to subdivide an element into several smaller cells. The CBCs are defined based on the methods of subdividing a hexahedral mesh element viewing from different directions.
Fig. 2. Concept of dynamic mesh generation.
248
L. Wang, T. Moriwaki / Precision Engineering 27 (2003) 245–257
gies of the elements or moving their nodes within the border zone to ensure a correct nodal connectivity. Fig. 4 illustrates the mechanism of the CBC substitution approach for both the localized mesh alignment (initial mesh generation) and the dynamic mesh adjustment through a three-element simple example. As shown in the lower-left of the figure, the six CBCs defined in Fig. 3 are used to assume the substitution. There are two steps to be followed in order to adjust an ill connection between two elements, the mid element and the right element in this case, as shown in the middle-left of Fig. 4: 1. Choose an appropriate CBC that has the same topological structure to be achieved, and 2. Substitute the CBC with the initial (thick-lined) mesh element.
Fig. 3. Topological structures of CBCs.
The modified mesh is shown in the lower-right of the figure. It is worthy to mention that a pre-defined CBC is substituted ‘virtually’ with a real element to replace its topology. The dimension and orientation of the CBC are adjusted (rotated, scaled, or sheared) relatively according to real situations. The algorithm used for CBC substitution is based on the earlier study of the first author [39], and presented in the following subsections in detail.
The first row and the first column of the figure describe how to subdivide an element, if viewing from Z and X directions. In Fig. 3, the Y-axis coincides the normal vector of the contacting surface (border zone) of any two primitives. The X- and Z-axes show the directions along which the hexahedral element is subdivided. Since CBC34 is identical to CBC43 in topological structure, only the six CBCs within the bold line are used for mesh adjustment. If a mesh element is required to be subdivided, it is simply substituted by an appropriate CBC to change its topology. This substitution method is named CBC substitution [9] by the authors, and presented in detail in Section 4.
Assume that two primitives (a fixed primitive F with larger mesh size and a movable primitive M with smaller mesh size) contact each other on a flat surface. The contacting surface is the border zone of the mesh adjustment in this case. The issue is how to cut the large mesh elements of F, near the border zone, into smaller ones and match them with the elements of M. The procedures consist of five steps.
4. CBC substitution
Step 1. Mesh regeneration for primitives with extra large elements.
The mesh elements in the border zone of two adjacent primitives are adjusted locally by either changing the topolo-
This procedure applies to primitive F only if its mesh size δF is greater than two times of the mesh size δM of primitive
4.1. Basic procedures of CBC substitution
Fig. 4. Mesh adjustment by CBC substitution.
L. Wang, T. Moriwaki / Precision Engineering 27 (2003) 245–257
M. If it is the case, the mesh elements of F are cut to satisfy the following condition: 1 ≤ δF /δM ≤ 1.5. Step 2. Extraction of elements and nodes to be modified. The nodes that are not shared by both primitives F and M in the border zone are detected, and the corresponding elements are extracted based on the geometry data of the primitives. The coordinates of the unmatched nodes must be adjusted, or the extracted elements must be substituted with a CBC, to make one-to-one correspondence between the nodes of F and M. Step 3. Determination of N-Labels. N-Labels (node labels) specify the relationships between the nodes of F and M, and provide a base for localizing mesh adjustment. Initial values of 1/2 are assigned to the nodes inside of the border zone, while keeping 0 as the N-Label for other nodes outside of border zone. These initial values are then modified based on the relative positions of the nodes of F against those of M (see Section 4.2 for further details). The value of an N-Label L(Vi ) describes the following conditions for node Vi of the primitive F. 0 Vi is outside of border zone. Not to be adjusted (initial value). 1 Vi is inside of border zone. 2 To be adjusted (initial value). n Vi is inside of border zone, and has L(Vi ) = been matched with node Vn of M. 1 n + 2 Vi is inside of border zone. A new node to be created between Vn and Vn+1 of M. (1) Some examples of N-Labels are illustrated in Fig. 5. After initial values are assigned to the nodes around the border zone, mesh adjustment will focus on the nodes with a label of 1/2. This figure shows a case where the node Vi of F has just found its corresponding node Vj of M to be matched.
249
Step 4. Determination of F-Labels. F-Labels (face labels) are assigned to the elements that are in need of fine adjustments or subdivisions. Corresponding to the numbers shown in the first row and the first column of Fig. 3, an F-Label specifies how the face of an element should be subdivided. The value of an F-Label Sj is as follows: 0 the face is outside of border zone (initial value). 1 the face should not be subdivided (initial value). Sj = 2 the face should be divided into two subareas. 3 the face should be divided into three subareas. 4 the face should be divided into four subareas. (2) Each mesh element happening to be adjusted near border zone should have two F-Labels, which specify the ways of subdividing the element in two different directions. The F-Labels are determined based on the N-Labels of the element (see Section 4.3 for further details). Step 5. Selection of suitable CBCs. A suitable type of CBC is selected for each of the mesh elements that require subdivisions, by referring to the F-Labels of the element. The element is then substituted with the selected CBC to replace its topology and to increase/decrease its subdivisions. The nodal coordinates of the newly modified elements are changed accordingly to match nodes of different primitives across the border zone. This final CBC substitution is applied only to the mesh elements near the border zone. Most of the initial mesh elements are not affected by this operation. Although it is possible to divide elements into same size for the two components to achieve the ease of nodal matching, the node-splitting problem of moving components still needs to be solved to allow boundary conditions to pass from one element to another. As a simplified case, it is covered by the same algorithm. 4.2. Determination of N-Labels Initial values of N-Labels are assigned to the nodes of primitives F and M around their border zone, based on the
Fig. 5. Examples of N-Labels during node adjustment.
250
L. Wang, T. Moriwaki / Precision Engineering 27 (2003) 245–257
following equations. LF (ViF ) = LM (VjM ) =
1 2
0 1 2
0
(Continued ) M ) ∧ LM (V M ) ≡ 1 (case 7 in general) if exist(Vj−1 j−1 2
if ViF ∈ V M ∩ V F otherwise if VjM ∈ V M ∩ V F otherwise
F , V F , V F ); create between(Vi−1 i i−1 F ) = j − 1, M ) = i − 1; LF (Vi−1 LM (Vj−1 2
(3)
The initial values are set to be 1/2, if the nodes are inside of the border zone between primitive F and M. Otherwise, they are set to be 0. The nodes with N-Labels of 0s are not considered in the following process, unless affected by other nodes. For a node with an N-Label of 1/2, it is first adjusted based on the nearest approximation to the nodes of the adjacent primitive across the border zone. In the case that no correspondence can be found, a new node is created, as shown in Fig. 5, to ensure one-to-one matched nodal connectivity. Fig. 6 summarizes the procedures of how to find or create a node in seven possible cases. For example, the node VjM is moved to and matched with ViF in case 1. Their N-Labels are changed to i and j, respectively, to reflect this adjustment. On the other hand, as shown in the case of 5(a) of Fig. 6, M at the right edge of primitive M is left unthe node Vj+1 matched; a new node V Fi is therefore created to pair with M , instead of moving node V F backward. This simplifies Vj+1 i+1 the algorithm processing effectively for unidirectional node adjustment (from left to right along the border). The algorithm of N-Label assignment is developed to cover those possible cases illustrated in Fig. 6. It is worthy to mention that this operation is localized in the border zone and performed from one end to the other along the border. Its processing procedure is given in the following C-like language, for ∀ViF (ViF ∈ V F ∧ LF (ViF ) ≡ 21 ) do M )∃V F ; search a VjM ∈ V M that is between(VjM , Vj+1 i F sw = check status(Vi ); switch (sw) ( see Fig. 6) case 1: LM (VjM ) ≡ 0 ∧ is near(VjM )∃ViF move to(VjM , ViF ); LM (VjM ) = i, LF (ViF ) = j; M )∃V F case 2: LM (VjM ) ≡ 0 ∧ is near(Vj+1 i ...... M )∃V F case 6: LM (VjM ) ≡ 21 ∧ is near(Vj+1 i M ); move to(ViF , Vj+1 M ) = i; LF (ViF ) = j + 1, LM (Vj+1 F ) ≡ 0 ∧ is near(V M )∃V F if LF (Vi−1 j i−1 F ,V M ); move to(Vi−1 j F ) = j, LM (VjM ) = i − 1; LF (Vi−1 else F , V F , V F ); create between(Vi−1 i i−1 F F M L (Vi−1 ) = j, L (VjM ) = i − 21 ; end if end switch
end if end for where, is between, check status, is near, move to, create between, and exist are defined macro functions for nodal status checking and node alignment. 4.3. Determination of F-Labels The values of F-Labels of a primitive are largely determined based on the N-Labels of another primitive across their border. The value assignments are carried out along two directions; they are, the direction designated by P (the direction of the motion between the primitives F and M), and the direction designated by V (the direction perpendicular to P and the normal vector of the contacting surface shared by F and M). Fig. 7 shows a typical example of the process of F-Label determination. The first row in the figure gives schematic illustration of the N-Labels of M and the nodal connectivity between the nodes of M and F across the border. The initial values of F-Labels are assigned to each mesh element of F as shown in the second row in the figure, using the following equations. As the mesh size of F is larger than that of M, only the mesh elements of F should be subdivided in this case. S F (i)
=
S M (j)
=
0 1
if N-Label of ith node LF (Vi F ) = 0 otherwise
0 1
if N-Label of jth node LM (Vj M ) = 0 otherwise (4)
The N-Labels of M are then used to adjust the initial F-Labels of F. As shown in the third row, some F-Labels are changed from 1 to 3 by following the procedure outlined below, for ∀VjM ∈ V M do if LM (VjM ) ≡ ξ + S F (ξ) = 3; end if end for
1 2
∧ ξ ∈ NF
where, ξ is a parameter giving the ID numbers of nodes, and NF is the node set of F. The elements with F-Labels of 3 are paired together (see Fig. 7). The F-Labels of those nodes in between the pairs are changed from 1 to 2 as shown in the fourth row of the figure, by applying the following procedure.
L. Wang, T. Moriwaki / Precision Engineering 27 (2003) 245–257
251
Fig. 6. Seven possible cases of node adjustment.
for ∀i ∈ N F that S F (i) ≥ 1 do search a i that S F (i) ≡ 3;; search a λ that S F (i + λ) ≡ 3 ∧ (i + λ) ∈ N F ; for ∀ (1 < < λ) do S F (i + γ) = 2; end for end for search a i if that S F (i) ≡ 0 ∧ S F (i − 1) ≡ 2 do; S F (i) = 3; end all
The final results are given in the last row of Fig. 7. The procedures mentioned above are carried out along the two directions designated by P and V to cover three-dimensional cases of mesh adjustments. The same procedure can also be used to find out the F-Labels of M, if needed, based on the N-Labels of F. 4.4. Selection of suitable CBCs Suitable CBCs are selected for those mesh elements whose nodes disagree with each other, by referring to the F-Labels
252
L. Wang, T. Moriwaki / Precision Engineering 27 (2003) 245–257
Fig. 7. F-Label determination through an example.
of the elements in the directions of P and V. Assuming that SP and SV are set of F-Labels of the elements in the two directions, the ID number of an appropriate CBC can be obtained by using the following equation. IDCBC = {(SP (x), SV (y))|SP (x) ∈ SP ∧ SV (y) ∈ SV } (5)
replaced elements are adjusted properly to achieve smooth nodal connectivity. For example, assuming that the F-Label SPF (x) = 3 for the xth element of F in the direction of P and SVF (y) = 2 for the yth element of F in the direction of V, the ID of an appropriate CBC to be substituted with the (x, y)th element can be found as follows. def
The individual mesh elements, which are identified for subdivisions, are then simply substituted with the CBCs to modify their topologies, and the nodal coordinates of the
IDFCBC (x, y) = (SPF (x), SVF (y)) = (3, 2)→32
(6)
It is worthy of mentioning that the elements with ID of 11 are not changed in topology. This can largely eliminate the number of elements being affected during dynamic mesh adjustment and hence increase the efficiency. Fig. 8 shows one example of dynamic mesh adjustment.
5. Extension of CBC substitution 5.1. Concept of two-space meshing
Fig. 8. Example of dynamic mesh adjustment.
The CBC substitution approach described in Section 4 is basically developed for the geometry with plane surfaces. The algorithm is proven to be simple and effective when primitives contact with each other on plane border zones. An extension of the approach is needed to handle the more complex cases of mesh adjustment for arbitrary primitives with curved surfaces, or curved primitives. Mapping technique in projective geometry [40] is adopted for integrating the curved primitives with the CBC substitution approach. Fig. 9 illustrates the concept of two-space meshing and the extended CBC substitution through a meshing example of metal-milling. As shown in Fig. 9, the geometric model of a tool-workpiece is first generated in real object space R3 , where all the actual models are represented. Simultaneously, the geometric model is transformed to its corresponding projective space RP3 , where initial mesh elements are generated to each
L. Wang, T. Moriwaki / Precision Engineering 27 (2003) 245–257
253
5.2. Mapping and inverse mapping Mapping and inverse mapping bridge the real object space and its projective space, making the two-way projective transformations between a real model and its image possible. Fig. 10 shows the relationship between mapping and inverse mapping, and how the mapping functions can be found. A circle C2 = {(x, y) ∈ R2 |x2 + y2 ≤ 1} in two-dimensional Euclidean space R2 (real object space) and a square S 2 = {(ξ, η) ∈ RP2 ||ξ| ≤ 1, |η| ≤ 1} in projective space RP2 are taken as an example. Since both circle and square are compact set, continuous projection exists between C2 and S2 . The procedures of finding a mapping function f:C2 → S 2 and its inverse mapping function g:S 2 → C2 are as follows. (a) Draw a projection line outward from any point P that is inside of both C2 and S2 . (b) Find the intersections S and I with C2 and S2 , respectively. (c) Derive the relationship between S and I so as to establish the functions f and g. It is not difficult to derive the mapping and inverse mapping functions for this simple example as: x2 + y 2 (λx, λy) λ = (x, y) = (0, 0) f(x, y) = max{|x|, |y|} (0, 0) (x, y) = (0, 0) |η|} (λξ, λη) λ = max{|ξ|, (ξ, η) = (0, 0) 2 g(ξ, η) = ξ + η2 (0, 0) (ξ, η) = (0, 0) Fig. 9. Mapping for extended CBC substitution.
of the projected primitives (projective images). Since the curved surfaces can be transformed to plain ones by using appropriate mapping functions, the needed mesh adjustment can be done easily in RP3 by the proposed CBC substitution approach. Upon completion, the projected primitives are transformed inversely to R3 , by selecting the corresponding inverse mapping functions properly. It should be emphasized that the mesh adjustment is done in RP3 while the FEM analysis or simulation is performed in R3 , by repeating the process of Step 3, shown in Fig. 9, whenever relative motion occurred between the primitives.
(7) Since f · g = g · f = 1, if regarding (ξ, η) as (x, y), the circle C2 is said topological isomorphic to the square S2 or C2 ∼ = S 2 . This fact is of vital importance to the uniformity of projective transformation between R2 and RP2 without any ill cases (e.g. f: X → Y , g: Y → Z, but X = Z). Two pairs of commonly used mapping and inverse mapping functions are defined and listed in Table 1. The mapping functions of other geometries (cube, cone, torus, or any other analytical surfaces) can be found in the same way. These functions are used to calculate the node coordinates during mappings.
Fig. 10. Mapping and inverse mapping.
254
L. Wang, T. Moriwaki / Precision Engineering 27 (2003) 245–257
Table 1 Typical mapping and inverse mapping functions Cylinder
x2 + y 2 (λx, λy, z), λ = max{|x|, |y|} f(x, y, z) = (0, 0, z) max{|ξ|, |η|} (λξ, λη, ζ), λ = ξ 2 + η2 g(ξ, η, ζ) = (0, 0, ζ)
(x, y, z) = (0, 0, z) (x, y, z) = (0, 0, z) (ξ, η, ζ) = (0, 0, ζ) (ξ, η, ζ) = (0, 0, ζ)
Sphere
x2 + y2 + z2 (λx, λy, λz) λ = (x, y, z) = (0, 0, 0) max{|x|, |y|, |z|} f(x, y, z) = (0, 0, 0) (x, y, z) = (0, 0, 0) max{|ξ|, |η|, |ζ|} (λξ, λη, λζ), λ = (ξ, η, ζ) = (0, 0, 0) ξ 2 + η2 + ζ 2 g(ξ, η, ζ) = (0, 0, 0) (ξ, η, ζ) = (0, 0, 0)
5.3. Algorithm of extended CBC substitution As shown in Fig. 9, the dynamic mesh adjustment of a model with curved primitives includes projective transformations (mapping and inverse mapping), initial mesh alignment, and dynamic mesh adjustment. The processing algorithm of the extended approach is divided into three steps accordingly and described in a C-like language in detail. Step 1. Establishing a projective space and mapping the original model to the projective space. The procedure in Step 1 is to transform all curved primitives from R3 to RP3 through well-defined mapping functions. Initial meshes are generated to the projected primitives (images) before meshes are aligned across borders in the next step. for ∀pr(pr ∈ R3 ∧ pr ∈ {Curved Primitive Set}) do Mid ⇐ select mapping function(pr); map to(Mid , pr, Ipr ); generate initial mesh(Ipr ); pushIpr → RP3 ; for ∀pr(pr ∈ R3 ∧ pr ∈ / {Curved Primitive Set}) do generate initial mesh(pr); duplicate pr → RP3 ; Step 2. Adjusting initial meshes and mapping them inversely to the real space. The initial meshes generated in the previous step are adjusted locally in the projective space. A complete mesh model for analysis is obtained by inverse mapping. Any meshing data changes should be reflected in both spaces.
for ∀Ipr (Ipr ∈ RP3 ) do adjust mesh(Ipr , CBC); refreshIpr → RP3 ; Mid ⇐ select inverse mapping function(Ipr ); inverse map to(Mid , Ipr , pr); refresh pr → R3 ; Step 3. Dynamic mesh adjustment between the two spaces. The relative positions among model components may vary with time during dynamic FEA or simulation. It is necessary to modify the FEM mesh model dynamically whenever the relative positions are changed. The dynamic mesh adjustment is described as follows. Whenever (relative motion occured==true) for ∀pr(pr ∈ R3 ) do status ⇐ check motion status(pr); if (status==moved) Ipr ⇐ select projected object(pr); adjust mesh;(Ipr , CBC) refresh Ipr →RP3 ; Mid ⇐ select inverse mapping function(Ipr ); inverse map to(Mid , pr, Ipr ); refresh pr→R3 ; Special attentions to those nodes on the border zones of the primitives should be given to the inverse mapping in Steps 2 and 3. Fig. 11 illustrates the nodal connectivity before and after the inverse mapping. Since the nodes indicated by 䊊 differ in nodal connectivity after the inverse mapping, each of them are divided into two and assigned to the primitives that previously shared the node. The other nodes given by remain unchanged. Whether to divide a node or not is
L. Wang, T. Moriwaki / Precision Engineering 27 (2003) 245–257
255
three nodes approach to collinear, the accuracy of FEA may decrease. To prevent it from happening, we use 20-node elements instead of 8-node elements for FEA calculations. The details of how to avoid this problem and what techniques to be applied will be reported separately. 6. Case study The proposed CBC substitution is developed for dynamic analysis using finite element method (FEM). The FEM is known as a general technique for constructing an approximation of solution over a collection of finite elements, by using variation concept. When associated with the CBC substitution while finite elements are changing, one problem to be solved is how to handle the nodal analytical data for those adjusted and newly generated nodes. Since the thermal behavior of machines is more influential than most of the other factors [41], a method for dynamic thermal analysis is considered in this case study. The basic procedures are: Fig. 11. Node division during inverse mapping.
simple—it is judged by the nodal connectivity of the sharing primitives in the real object space R3 . After inverse mapping as shown in Fig. 11, there are cases that three nodes of an element share one arc segment. If the
1. Calculating temperature distribution at time t, and recording all nodal analytical data including nodal temperature Tt and its incremental rate T˙ t . 2. Adjusting nodal connectivity when the positions of some elements are changed during a time interval of (t, and recalculating those nodal coordinates accordingly. 3. Deriving and interpolating new nodal temperature for those affected nodes based on Tt and T˙ t .
Fig. 12. Temperature distribution and deformation with moving heat source.
256
L. Wang, T. Moriwaki / Precision Engineering 27 (2003) 245–257
4. Calculating temperature distribution Tt +(t and T˙ t+(t at time t + (t. The above procedures are repeated until the thermal analysis is completed. However, Steps 2 and 3 are needed only if the dynamic mesh adjustment is involved in that cycle. An interpolating algorithm [42] consisting of a primary interpolation and a secondary interpolation is used in Step 3. Generally, there exists a complex phenomenon across the interface between two movable components. How to handle the problem is considered critical to the accuracy of a finite element analysis. Instead of using a specific type of finite elements, we developed a concept called virtual interface [43], across which a virtual heat conduction is taken place between the contacting surfaces to simulate the heat conduction in reality. Frictional heat between relatively moving components is also taken into consideration, as it may serve as an additional heat source. With this virtual interface introduced to accumulate the complex behavior across components, our meshing algorithm becomes generic and simple. A good convergence in this case study is achieved manually by changing step intervals between calculations. Fig. 12 shows nine typical analytical results of both the temperature distribution and the thermal deformation of a simple table-base model of machine tools, by utilizing the methods proposed. In this case, a heat source is fixed on the top surface of the table that moves one round-trip over the base in 30 min during the thermal analysis. Although the case study is given to a simple model with only flat surfaces, other complicated geometries can be dealt with in the same way in terms of analytical procedures, because the major problem of dynamic meshing is handled by the extended CBC substitution approach. Further details on this will be reported separately.
7. Conclusions A new approach called CBC substitution is proposed to deal with initial alignments of constructive mesh generation and dynamic mesh adjustments during FEM analysis. A mapping technique is adopted to transform curved geometries to plane ones, so as to handle the complicated geometries in a consistent way. The meshing process is localized. Only the disagreed elements are adjusted dynamically. Common to other meshing algorithms, this approach also has its limitation. For a complex geometry, if unique mapping and inverse mapping functions cannot be found, the geometry must be decomposed and approximated until it can be described by certain mathematic equations. Other than this limitation, the proposed methodology is found effective and efficient in dealing with dynamic meshing situations. It is expected that the new methodology can be applied to dynamic analyses and simulations for better product design or more accurate device control based on their real life behaviors captured by the simulations.
References [1] Liu WK, Jun S, Li S, Adee J, Belytschko T. Reproducing kernel particle methods for structural dynamics. Int J Numeric Methods Eng 1995;28:1655–79. [2] Lu YY, Belytschko T, Gu L. A new implementation of the element free Galerkin method. Comput Methods Appl Mech Eng 1994;113:397– 414. [3] Oñate E, Idelsohn S. A mesh-free finite point method for advectivediffusive transport and fluid flow problems. Computat Mech 1998;21:283–92. [4] Liszka T, Orikisz J. Finite difference method at arbitrary irregular grids and its application in applied mechanics. Comput Struct 1980;11:83– 95. [5] Mukherjee YX, Mukherjee S. The boundary node method for potential problems. Int J Numeric Methods Eng 1997;40:797–815. [6] Duarte CA, Oden JT. An h-p adaptive method using clouds. Comput Methods Appl Mech Eng 1996;139:237–62. [7] Zhu T, Zhang J, Atluri SN. A meshless local boundary integral equation (LBIE) method for solving nonlinear problems. Computat Mech 1998;22:174–86. [8] Aluru NR, Li G. Finite cloud method: a true meshless technique based on a fixed reproducing kernel approximation. Int J Numeric Methods Eng 2001;50:2373–410. [9] Moriwaki T, Sugimura N, Wang L. A modelling system for finite element analysis of machine products. Trans North Am Manufact Res Inst SME 1993;21:383–9. [10] International Meshing Roundtable. http://www.imt.sandia.gov/. [11] Cavendish JC. Automatic triangulation of arbitrary planar domains for the finite element method. Int J Numeric Methods Eng 1974;8:679–96. [12] Sadek EA. A scheme for the automatic generation of triangular finite elements. Int J Numeric Methods Eng 1980;15:1813–22. [13] Liu Y, Chen K. A two-dimensional mesh generator for variable order triangular and rectangular elements. Comput Struct 1988;29(6):1033– 53. [14] Sarrate J, Huerta A. Efficient unstructured quadrilateral mesh generation. Int J Numeric Methods Eng 2000;49:1327–50. [15] Van Phai N. Automatic mesh generation with tetrahedron elements. Int J Numeric Methods Eng 1982;18:273–89. [16] Cavendish JC, Field DA, Frey WH. An approach to automatic three-dimensional finite element mesh generation. Int J Numeric Methods Eng 1985;21:329–47. [17] Mochizuki T, Kojima T. A shape data processing system for 3-dimensional finite element analysis using solid model. J Jpn Soc Precision Eng 1986;52(11):1910–5. [18] Calvo NA, Idelsohn SR. All-hexahedral element meshing: generation of the dual mesh by recurrent subdivision. Comput Methods Appl Mech Eng 2000;182:371–8. [19] Lo SH. Finite element mesh generation over curved surfaces. Comput Struct 1988;29(5):731–42. [20] Cheng F, Jaromczyk JW, Lin JR, Chang SS, Lu JY. A parallel mesh generation algorithm based on the vertex label assignment scheme. Int J Numeric Methods Eng 1989;28:1429–48. [21] Karamete BK, Beall MW, Shephard MS. Triangulation of arbitrary polyhedra to support automatic mesh generators. Int J Numeric Methods Eng 2000;49:167–91. [22] Dey S, O’Bara RM, Shephard MS. Towards curvilinear meshing in 3D: the case of quadratic simplices. Comput Aided Design 2001;33:199– 209. [23] Owen SJ, Saigal S. H-Morph an indirect approach to advancing front hex meshing. Int J Numeric Methods Eng 2000;49:289–312. [24] Owen SJ. Hex-dominant mesh generation using 3D constrained triangulation. Comput Aided Design 2001;33:211–20. [25] Owen SJ, Staten ML, Canann SA, Saigal S. Q-Morph: an indirect approach to advancing front quad meshing. Int J Numeric Methods Eng 1999;44:1314–40.
L. Wang, T. Moriwaki / Precision Engineering 27 (2003) 245–257 [26] Lu Y, Gadh R, Tautges TJ. Feature based hex meshing methodology: feature recognition and volume decomposition. Comput Aided Design 2001;33:221–32. [27] Li XY, Teng SH, Üngör A. Simultaneous refinement and coarsening for adaptive meshing. Eng Comput 1999;15:280–91. [28] Halpbern M. Industrial requirements and practices in finite element meshing: a survey of trends. In: Proceedings of 6th International Meshing Roundtable, SAND97-2399, Sandia National Laboratories, 1997. [29] Sheffer A, Bercovier M. Hexahedral meshing of non-linear volumes using Voronoi faces and edges. Int J Numeric Methods Eng 2000;49:329– 51. [30] Lai M, Benzley S, White D. Automated hexahedral mesh generation by generalized multiple source to multiple target sweeping. Int J Numeric Methods Eng 2000;49:261–75. [31] Staten ML, Canann SA, Owen SJ. BMSweep: locating interior nodes during sweeping. Eng Comput 1999;15:212–8. [32] Knupp P. Applications of mesh smoothing: copy, morph, and sweep on unstructured quadrilateral meshes. Int J Numeric Methods Eng 1999;45:37–45. [33] Tautges TJ. The generation of hexahedral meshes for assembly geometry: survey and progress. Int J Numeric Methods Eng 2001;50:2617–42. [34] Dhondt G. Unstructured 20-node brick element meshing. Comput Aided Design 2001;33:233–49.
257
[35] Dhondt G. A new automatic hexahedral mesher based on cutting. Int J Numeric Methods Eng 2001;50:2109–26. [36] Tautges T, Blacker T, Mitchell S. The Whisker Weaving algorithm: a connectivity-based method for constructing all-hexahedral finite element meshes. Int J Numeric Methods Eng 1996;39:3327–49. [37] Folwell NT, Mitchell SA. Reliable Whisker Weaving via Curve Contraction. Eng Comput 1999;15:292–302. [38] Ho-Le K. Finite element mesh generation methods: a review and classification. Comput Aided Design 1988;20(1):27–38. [39] Wang L. Development of an integrated CAD/CAE system for machine tool design. Doctoral Dissertation, Kobe University, Japan, March 1993. [40] Yokota I. From topological geometry to projective geometry. Modern Mathematics Press, 1993. [41] Bryan J. International status of thermal error research. Annals of the CIRP 1990;39(2):645–56. [42] Wang L. An interpolating algorithm for dynamic thermal analysis of machine tools. In: Proceedings of the First International Symposium on Advances in Intelligent Computer Integrated Manufacturing System; 1994. p. 213–8. [43] Moriwaki T, Wang L. Study on thermal behaviour of machine tool elements across movable joints under operating states. In: Proceedings of International Conference on Precision Engineering; 1995. p. 150–4.
A novel meshing algorithm for dynamic finite element analysis Lihui Wang a,∗ , Toshimichi Moriwaki b a
Integrated Manufacturing Technologies Institute, National Research Council of Canada, 800 Collip Circle, London, Ont., Canada N6G 4X8 b Department of Mechanical Engineering, Kobe University, 1-1 Rokkodai-cho, Nada-ku, Kobe 657-8501, Japan Received 12 July 2002; received in revised form 31 October 2002; accepted 11 December 2002
Abstract This paper describes a new algorithm to handle problems in dynamic finite element analysis and run-time simulation, where mesh re-generation or dynamic adjustment is required. Based on a concept called coded box cell (CBC) substitution, this algorithm can be applied to both initial mesh generation and dynamic mesh adjustment along the border zones of multiple primitives that form an entire model. During the initial mesh generation, appropriate labels are assigned to the nodes and the faces of each finite element. These labels are used to facilitate decision-making in dynamic mesh adjustment. A mapping technique is adopted to transform curved surfaces to plain surfaces for the ease of automatic mesh adjustment while still using the same algorithm. The results of a case study show that a finite element mesh can be adjusted dynamically and locally around its border zone; and the algorithm can be utilized effectively to simulate the thermal behavior of a device under real operating conditions. Crown Copyright © 2003 Published by Elsevier Science Inc. All rights reserved. Keywords: Meshing; CAE; Finite element analysis; Mapping and inverse mapping; Thermal simulation
1. Introduction Today, simulation of production processes is becoming much more important in manufacturing, where complex machining processes and material handling are encountered. A particular emphasis has been given to research and development of virtual manufacturing having various computer-aided engineering (CAE) packages for analysis and simulation of machine behavior and dynamics, aiming at realizing an optimized production, right the first time, with high quality, low cost, and short lead time. It often requires an advanced system capability to analyze and simulate the dynamic behaviors of products, machines, or production lines, as if they are under real operating conditions. In the last two decades, a large number of CAE systems have been developed and put into practical use in product development, such as analyses of stress, strain, vibration, and thermal behavior of the product by applying finite element method (FEM) and boundary element method (BEM). Recently, there is a growing interest in the development of meshless methods for numerical solutions of partial differential equations. Some typical meshless methods include the reproducing kernel particle method (RKPM) [1], element ∗ Corresponding
author. Tel.: +1-519-430-7084; fax: +1-519-430-7064. E-mail address: [email protected] (L. Wang).
free Galerkin method (EFG) [2], finite point method (FPM) [3], and others [4–8]. However, considering the accuracy and stability of the analytical results, FEM remains as the most popular and practical approach available today for solving various complex engineering problems. Unlike the meshless methods, FEM or BEM inherently requires a properly meshed geometry, before any analysis can take place. Iterative analysis usually requires the mesh being partially refined or even regenerated to increase analytical accuracy. Despite a large number of meshing methods developed in the last three decades, algorithms for dynamic mesh generation and localized mesh adjustment are still missing from the literature. It is urgently required to establish a systematic approach that not only can handle the re-meshing problem effectively, but also can simulate a machine’s behavior dynamically and accurately, especially under real operating conditions. A complete understanding of a machine’s dynamics or its real life behavior is vital to both the structural design and the run-time machine control. The objective of this research is to develop an integrated design/analysis system that is capable of dynamic simulation of machines in real life. This paper introduces a new algorithm that is able to handle problems in FEM analysis and simulation, where mesh re-generation or dynamic adjustment is required. Based on a concept called coded box cell (CBC) substitution [9], this algorithm can be used to adjust those finite elements whose nodes disagree with each other during
0141-6359/03/$ – see front matter. Crown Copyright © 2003 Published by Elsevier Science Inc. All rights reserved. doi:10.1016/S0141-6359(03)00005-9
246
L. Wang, T. Moriwaki / Precision Engineering 27 (2003) 245–257
Nomenclature M F D
movable object fixed object direction of normal vector on the contacting plane of the objects F and M P direction of motion of the object M against F V direction perpendicular to both the directions P and D N(φ) number of nodes in the direction P or V in the FEM mesh of object φ δF , δM mesh sizes (lengths between adjacent nodes) of the objects F and M, respectively VF , VM sets of nodes of the objects F and M along their border zone, respectively LF , LM node labels of the objects F and M, respectively SF , SM face labels of the objects F and M, respectively IDCBC set of identification numbers of CBCs to be used for substitution R3 real object space RP3 projective space (x, y, z) coordinates in the real object space R3 (ξ, η, ζ) coordinates in the projective space RP3 f mapping function g inverse mapping function Mid mapping function identity pr instance of primitives Ipr image of pr
analysis at their contacting surfaces (border zones), in addition to the initial mesh generation. Following a brief review and classification of the available meshing methods, the processing procedures of the algorithm are described in detail. A mapping technique is adopted to transform curved surfaces to plain ones for the ease of automatic mesh adjustment by applying the same algorithm. The algorithm is verified through an example. The results of the case study show that an FEM mesh can be adjusted dynamically and locally around its border zone; and the algorithm can be utilized effectively to simulate its thermal behavior under real operating conditions.
2. A brief review Most of mesh generation tasks involve a great deal of effort to generate numerical data, such as node numbers, nodal point coordinates, parameter values, and element connectivity. This work could be very time-consuming and error-prone if done manually, especially when a large number of elements or very complex geometry are encountered. There are many bottlenecks in the meshing process: geometry cleanup, geometry simplification or de-featuring, and finite element mesh
generation. During the last few decades, much research have been carried out to automate the FEM mesh generation process [10], based on the geometric models created by using CAD systems. Initial efforts in the literature were focused on the twodimensional mesh generation algorithms based on triangulation [11–13] and quadrilateral elements [13,14]. Attention then shifted to the meshing algorithms for three-dimensional geometry [15–18] and geometry with curved surfaces [18–22] using tetrahedra, hexahedra, and other polyhedra like triangular prisms. Despite those accomplishments in the last few decades, research on new meshing algorithms remains active. Owen et al. [23,24] use an algorithm called H-Morph to generate a hexahedral-dominant finite element mesh for arbitrary volumes. As an extension of the Q-Morph algorithm [25], the H-Morph method starts with an initial tetrahedral mesh and systematically transforms and combines tetrahedra into hexahedra, without the need of decomposing an arbitrary volume. On the contrary, Lu et al. [26] use a feature-based approach for volume decomposition before applying an appropriate and well-established meshing algorithm to each of the meshable pieces. Li et al. [27] presented a unified scheme for simultaneously refining and coarsening a mesh. This adaptive meshing algorithm is suitable for a domain that has a moving boundary. Recently, notable efforts have been given to the automated hexahedral mesh generation, because hexahedral mesh generally provides higher “artificial stiffness” and better solutions than the tetrahedral one for many physical problems [28]. Sheffer and Bercovier [29] introduced a hexahedral meshing algorithm based on the embedded Voronoi graph (EVG), where the EVG is used to decompose a non-linear volume into simple sub-volumes meshable by basic meshing techniques. Lai et al. [30] proposed an enhanced sweeping mesh generation algorithm by creating all hexahedral element meshes between multiple source surfaces and multiple target surfaces, while most other sweeping techniques [31–33] create all hexahedral element meshes by projecting an existing single-surface mesh along a specified trajectory to a specified single target surface. Dhondt presented a method that can generate an unstructured 20-node brick element mesh for arbitrary structure, based on a triangulation of the structure’s surface [34,35]. Starting from a hexahedral master mesh encompassing the structure, the elements that are intersected by the triangulation are determined, cut, and re-meshed according to their cutting topologies. Whisker Weaving is another advancing front algorithm for all-hexahedral mesh generation [36,37]. It uses global information derived from grouping the mesh dual into surfaces to construct the connectivity of the mesh, and then positions the nodes afterwards. Based on the previous literature surveys [28,33,38] and authors’ investigation, major meshing algorithms available in public domain can be classified and summarized, as shown in Fig. 1, in terms of the techniques used and the topologies applied. Most meshing approaches can be applied to 3D cases.
L. Wang, T. Moriwaki / Precision Engineering 27 (2003) 245–257
247
Fig. 1. Classification of mesh generation methods.
3. Concept of dynamic mesh generation When a milling machine is boring an engine block, the relative positions among its components (table, column, and spindle, etc.) vary with time. It is generally difficult to analyze and simulate the dynamic behavior of the milling machine using conventional FEA tools while its components are moving relatively, because most of the FEA tools require the mesh to be regenerated once the corresponding geometry model is changed. It becomes impossible if the geometry model keeps changing rapidly due to the relative movements between its components. In contrast, dynamic mesh generation does not require the entire mesh to be regenerated, and it considers mesh generation in a constructive way when the geometry model of an assembly having movable components is designed. The fundamental process of the concept is shown in Fig. 2. In the first step, a set of primitives (basic geometric shapes of a component or an assembly) can be generated interactively by using a solid modeling system, followed by the initial mesh generation for each primitive. Hexahedral elements are selected in this approach. When the primitives are combined in the third step, the elements near their border zone are adjusted or modified locally in the fourth step, so that the mesh
of one primitive connects with that of another smoothly across the border. This operation needs to be performed only once if the geometry model consists of non-movable primitives only. Otherwise, as shown in the fifth step, the mesh along the border of the movable primitives has to be adjusted on the fly whenever its relative position is changed. This dynamic mesh modification is limited to those finite elements whose nodes disagree with each other within the border zone and the elements being affected. It is not difficult to adjust the mesh along the border zone if the sizes of mesh elements (across the border) of two primitives are the same. The mesh can be modified, in this case, by adjusting the nodal coordinates of those conflicting elements in the border zone. However, in most cases, the mesh elements of one primitive must be subdivided in order that the numbers of the elements of both primitives agree with each other within the border zone. A requirement to be met in subdividing an initial mesh is to localize the modification only in the area near the border zone so as to increase the efficiency of dynamic adjustment. A set of pre-cut hexahedral cells, called CBC as shown in Fig. 3, are defined and used to subdivide an element into several smaller cells. The CBCs are defined based on the methods of subdividing a hexahedral mesh element viewing from different directions.
Fig. 2. Concept of dynamic mesh generation.
248
L. Wang, T. Moriwaki / Precision Engineering 27 (2003) 245–257
gies of the elements or moving their nodes within the border zone to ensure a correct nodal connectivity. Fig. 4 illustrates the mechanism of the CBC substitution approach for both the localized mesh alignment (initial mesh generation) and the dynamic mesh adjustment through a three-element simple example. As shown in the lower-left of the figure, the six CBCs defined in Fig. 3 are used to assume the substitution. There are two steps to be followed in order to adjust an ill connection between two elements, the mid element and the right element in this case, as shown in the middle-left of Fig. 4: 1. Choose an appropriate CBC that has the same topological structure to be achieved, and 2. Substitute the CBC with the initial (thick-lined) mesh element.
Fig. 3. Topological structures of CBCs.
The modified mesh is shown in the lower-right of the figure. It is worthy to mention that a pre-defined CBC is substituted ‘virtually’ with a real element to replace its topology. The dimension and orientation of the CBC are adjusted (rotated, scaled, or sheared) relatively according to real situations. The algorithm used for CBC substitution is based on the earlier study of the first author [39], and presented in the following subsections in detail.
The first row and the first column of the figure describe how to subdivide an element, if viewing from Z and X directions. In Fig. 3, the Y-axis coincides the normal vector of the contacting surface (border zone) of any two primitives. The X- and Z-axes show the directions along which the hexahedral element is subdivided. Since CBC34 is identical to CBC43 in topological structure, only the six CBCs within the bold line are used for mesh adjustment. If a mesh element is required to be subdivided, it is simply substituted by an appropriate CBC to change its topology. This substitution method is named CBC substitution [9] by the authors, and presented in detail in Section 4.
Assume that two primitives (a fixed primitive F with larger mesh size and a movable primitive M with smaller mesh size) contact each other on a flat surface. The contacting surface is the border zone of the mesh adjustment in this case. The issue is how to cut the large mesh elements of F, near the border zone, into smaller ones and match them with the elements of M. The procedures consist of five steps.
4. CBC substitution
Step 1. Mesh regeneration for primitives with extra large elements.
The mesh elements in the border zone of two adjacent primitives are adjusted locally by either changing the topolo-
This procedure applies to primitive F only if its mesh size δF is greater than two times of the mesh size δM of primitive
4.1. Basic procedures of CBC substitution
Fig. 4. Mesh adjustment by CBC substitution.
L. Wang, T. Moriwaki / Precision Engineering 27 (2003) 245–257
M. If it is the case, the mesh elements of F are cut to satisfy the following condition: 1 ≤ δF /δM ≤ 1.5. Step 2. Extraction of elements and nodes to be modified. The nodes that are not shared by both primitives F and M in the border zone are detected, and the corresponding elements are extracted based on the geometry data of the primitives. The coordinates of the unmatched nodes must be adjusted, or the extracted elements must be substituted with a CBC, to make one-to-one correspondence between the nodes of F and M. Step 3. Determination of N-Labels. N-Labels (node labels) specify the relationships between the nodes of F and M, and provide a base for localizing mesh adjustment. Initial values of 1/2 are assigned to the nodes inside of the border zone, while keeping 0 as the N-Label for other nodes outside of border zone. These initial values are then modified based on the relative positions of the nodes of F against those of M (see Section 4.2 for further details). The value of an N-Label L(Vi ) describes the following conditions for node Vi of the primitive F. 0 Vi is outside of border zone. Not to be adjusted (initial value). 1 Vi is inside of border zone. 2 To be adjusted (initial value). n Vi is inside of border zone, and has L(Vi ) = been matched with node Vn of M. 1 n + 2 Vi is inside of border zone. A new node to be created between Vn and Vn+1 of M. (1) Some examples of N-Labels are illustrated in Fig. 5. After initial values are assigned to the nodes around the border zone, mesh adjustment will focus on the nodes with a label of 1/2. This figure shows a case where the node Vi of F has just found its corresponding node Vj of M to be matched.
249
Step 4. Determination of F-Labels. F-Labels (face labels) are assigned to the elements that are in need of fine adjustments or subdivisions. Corresponding to the numbers shown in the first row and the first column of Fig. 3, an F-Label specifies how the face of an element should be subdivided. The value of an F-Label Sj is as follows: 0 the face is outside of border zone (initial value). 1 the face should not be subdivided (initial value). Sj = 2 the face should be divided into two subareas. 3 the face should be divided into three subareas. 4 the face should be divided into four subareas. (2) Each mesh element happening to be adjusted near border zone should have two F-Labels, which specify the ways of subdividing the element in two different directions. The F-Labels are determined based on the N-Labels of the element (see Section 4.3 for further details). Step 5. Selection of suitable CBCs. A suitable type of CBC is selected for each of the mesh elements that require subdivisions, by referring to the F-Labels of the element. The element is then substituted with the selected CBC to replace its topology and to increase/decrease its subdivisions. The nodal coordinates of the newly modified elements are changed accordingly to match nodes of different primitives across the border zone. This final CBC substitution is applied only to the mesh elements near the border zone. Most of the initial mesh elements are not affected by this operation. Although it is possible to divide elements into same size for the two components to achieve the ease of nodal matching, the node-splitting problem of moving components still needs to be solved to allow boundary conditions to pass from one element to another. As a simplified case, it is covered by the same algorithm. 4.2. Determination of N-Labels Initial values of N-Labels are assigned to the nodes of primitives F and M around their border zone, based on the
Fig. 5. Examples of N-Labels during node adjustment.
250
L. Wang, T. Moriwaki / Precision Engineering 27 (2003) 245–257
following equations. LF (ViF ) = LM (VjM ) =
1 2
0 1 2
0
(Continued ) M ) ∧ LM (V M ) ≡ 1 (case 7 in general) if exist(Vj−1 j−1 2
if ViF ∈ V M ∩ V F otherwise if VjM ∈ V M ∩ V F otherwise
F , V F , V F ); create between(Vi−1 i i−1 F ) = j − 1, M ) = i − 1; LF (Vi−1 LM (Vj−1 2
(3)
The initial values are set to be 1/2, if the nodes are inside of the border zone between primitive F and M. Otherwise, they are set to be 0. The nodes with N-Labels of 0s are not considered in the following process, unless affected by other nodes. For a node with an N-Label of 1/2, it is first adjusted based on the nearest approximation to the nodes of the adjacent primitive across the border zone. In the case that no correspondence can be found, a new node is created, as shown in Fig. 5, to ensure one-to-one matched nodal connectivity. Fig. 6 summarizes the procedures of how to find or create a node in seven possible cases. For example, the node VjM is moved to and matched with ViF in case 1. Their N-Labels are changed to i and j, respectively, to reflect this adjustment. On the other hand, as shown in the case of 5(a) of Fig. 6, M at the right edge of primitive M is left unthe node Vj+1 matched; a new node V Fi is therefore created to pair with M , instead of moving node V F backward. This simplifies Vj+1 i+1 the algorithm processing effectively for unidirectional node adjustment (from left to right along the border). The algorithm of N-Label assignment is developed to cover those possible cases illustrated in Fig. 6. It is worthy to mention that this operation is localized in the border zone and performed from one end to the other along the border. Its processing procedure is given in the following C-like language, for ∀ViF (ViF ∈ V F ∧ LF (ViF ) ≡ 21 ) do M )∃V F ; search a VjM ∈ V M that is between(VjM , Vj+1 i F sw = check status(Vi ); switch (sw) ( see Fig. 6) case 1: LM (VjM ) ≡ 0 ∧ is near(VjM )∃ViF move to(VjM , ViF ); LM (VjM ) = i, LF (ViF ) = j; M )∃V F case 2: LM (VjM ) ≡ 0 ∧ is near(Vj+1 i ...... M )∃V F case 6: LM (VjM ) ≡ 21 ∧ is near(Vj+1 i M ); move to(ViF , Vj+1 M ) = i; LF (ViF ) = j + 1, LM (Vj+1 F ) ≡ 0 ∧ is near(V M )∃V F if LF (Vi−1 j i−1 F ,V M ); move to(Vi−1 j F ) = j, LM (VjM ) = i − 1; LF (Vi−1 else F , V F , V F ); create between(Vi−1 i i−1 F F M L (Vi−1 ) = j, L (VjM ) = i − 21 ; end if end switch
end if end for where, is between, check status, is near, move to, create between, and exist are defined macro functions for nodal status checking and node alignment. 4.3. Determination of F-Labels The values of F-Labels of a primitive are largely determined based on the N-Labels of another primitive across their border. The value assignments are carried out along two directions; they are, the direction designated by P (the direction of the motion between the primitives F and M), and the direction designated by V (the direction perpendicular to P and the normal vector of the contacting surface shared by F and M). Fig. 7 shows a typical example of the process of F-Label determination. The first row in the figure gives schematic illustration of the N-Labels of M and the nodal connectivity between the nodes of M and F across the border. The initial values of F-Labels are assigned to each mesh element of F as shown in the second row in the figure, using the following equations. As the mesh size of F is larger than that of M, only the mesh elements of F should be subdivided in this case. S F (i)
=
S M (j)
=
0 1
if N-Label of ith node LF (Vi F ) = 0 otherwise
0 1
if N-Label of jth node LM (Vj M ) = 0 otherwise (4)
The N-Labels of M are then used to adjust the initial F-Labels of F. As shown in the third row, some F-Labels are changed from 1 to 3 by following the procedure outlined below, for ∀VjM ∈ V M do if LM (VjM ) ≡ ξ + S F (ξ) = 3; end if end for
1 2
∧ ξ ∈ NF
where, ξ is a parameter giving the ID numbers of nodes, and NF is the node set of F. The elements with F-Labels of 3 are paired together (see Fig. 7). The F-Labels of those nodes in between the pairs are changed from 1 to 2 as shown in the fourth row of the figure, by applying the following procedure.
L. Wang, T. Moriwaki / Precision Engineering 27 (2003) 245–257
251
Fig. 6. Seven possible cases of node adjustment.
for ∀i ∈ N F that S F (i) ≥ 1 do search a i that S F (i) ≡ 3;; search a λ that S F (i + λ) ≡ 3 ∧ (i + λ) ∈ N F ; for ∀ (1 < < λ) do S F (i + γ) = 2; end for end for search a i if that S F (i) ≡ 0 ∧ S F (i − 1) ≡ 2 do; S F (i) = 3; end all
The final results are given in the last row of Fig. 7. The procedures mentioned above are carried out along the two directions designated by P and V to cover three-dimensional cases of mesh adjustments. The same procedure can also be used to find out the F-Labels of M, if needed, based on the N-Labels of F. 4.4. Selection of suitable CBCs Suitable CBCs are selected for those mesh elements whose nodes disagree with each other, by referring to the F-Labels
252
L. Wang, T. Moriwaki / Precision Engineering 27 (2003) 245–257
Fig. 7. F-Label determination through an example.
of the elements in the directions of P and V. Assuming that SP and SV are set of F-Labels of the elements in the two directions, the ID number of an appropriate CBC can be obtained by using the following equation. IDCBC = {(SP (x), SV (y))|SP (x) ∈ SP ∧ SV (y) ∈ SV } (5)
replaced elements are adjusted properly to achieve smooth nodal connectivity. For example, assuming that the F-Label SPF (x) = 3 for the xth element of F in the direction of P and SVF (y) = 2 for the yth element of F in the direction of V, the ID of an appropriate CBC to be substituted with the (x, y)th element can be found as follows. def
The individual mesh elements, which are identified for subdivisions, are then simply substituted with the CBCs to modify their topologies, and the nodal coordinates of the
IDFCBC (x, y) = (SPF (x), SVF (y)) = (3, 2)→32
(6)
It is worthy of mentioning that the elements with ID of 11 are not changed in topology. This can largely eliminate the number of elements being affected during dynamic mesh adjustment and hence increase the efficiency. Fig. 8 shows one example of dynamic mesh adjustment.
5. Extension of CBC substitution 5.1. Concept of two-space meshing
Fig. 8. Example of dynamic mesh adjustment.
The CBC substitution approach described in Section 4 is basically developed for the geometry with plane surfaces. The algorithm is proven to be simple and effective when primitives contact with each other on plane border zones. An extension of the approach is needed to handle the more complex cases of mesh adjustment for arbitrary primitives with curved surfaces, or curved primitives. Mapping technique in projective geometry [40] is adopted for integrating the curved primitives with the CBC substitution approach. Fig. 9 illustrates the concept of two-space meshing and the extended CBC substitution through a meshing example of metal-milling. As shown in Fig. 9, the geometric model of a tool-workpiece is first generated in real object space R3 , where all the actual models are represented. Simultaneously, the geometric model is transformed to its corresponding projective space RP3 , where initial mesh elements are generated to each
L. Wang, T. Moriwaki / Precision Engineering 27 (2003) 245–257
253
5.2. Mapping and inverse mapping Mapping and inverse mapping bridge the real object space and its projective space, making the two-way projective transformations between a real model and its image possible. Fig. 10 shows the relationship between mapping and inverse mapping, and how the mapping functions can be found. A circle C2 = {(x, y) ∈ R2 |x2 + y2 ≤ 1} in two-dimensional Euclidean space R2 (real object space) and a square S 2 = {(ξ, η) ∈ RP2 ||ξ| ≤ 1, |η| ≤ 1} in projective space RP2 are taken as an example. Since both circle and square are compact set, continuous projection exists between C2 and S2 . The procedures of finding a mapping function f:C2 → S 2 and its inverse mapping function g:S 2 → C2 are as follows. (a) Draw a projection line outward from any point P that is inside of both C2 and S2 . (b) Find the intersections S and I with C2 and S2 , respectively. (c) Derive the relationship between S and I so as to establish the functions f and g. It is not difficult to derive the mapping and inverse mapping functions for this simple example as: x2 + y 2 (λx, λy) λ = (x, y) = (0, 0) f(x, y) = max{|x|, |y|} (0, 0) (x, y) = (0, 0) |η|} (λξ, λη) λ = max{|ξ|, (ξ, η) = (0, 0) 2 g(ξ, η) = ξ + η2 (0, 0) (ξ, η) = (0, 0) Fig. 9. Mapping for extended CBC substitution.
of the projected primitives (projective images). Since the curved surfaces can be transformed to plain ones by using appropriate mapping functions, the needed mesh adjustment can be done easily in RP3 by the proposed CBC substitution approach. Upon completion, the projected primitives are transformed inversely to R3 , by selecting the corresponding inverse mapping functions properly. It should be emphasized that the mesh adjustment is done in RP3 while the FEM analysis or simulation is performed in R3 , by repeating the process of Step 3, shown in Fig. 9, whenever relative motion occurred between the primitives.
(7) Since f · g = g · f = 1, if regarding (ξ, η) as (x, y), the circle C2 is said topological isomorphic to the square S2 or C2 ∼ = S 2 . This fact is of vital importance to the uniformity of projective transformation between R2 and RP2 without any ill cases (e.g. f: X → Y , g: Y → Z, but X = Z). Two pairs of commonly used mapping and inverse mapping functions are defined and listed in Table 1. The mapping functions of other geometries (cube, cone, torus, or any other analytical surfaces) can be found in the same way. These functions are used to calculate the node coordinates during mappings.
Fig. 10. Mapping and inverse mapping.
254
L. Wang, T. Moriwaki / Precision Engineering 27 (2003) 245–257
Table 1 Typical mapping and inverse mapping functions Cylinder
x2 + y 2 (λx, λy, z), λ = max{|x|, |y|} f(x, y, z) = (0, 0, z) max{|ξ|, |η|} (λξ, λη, ζ), λ = ξ 2 + η2 g(ξ, η, ζ) = (0, 0, ζ)
(x, y, z) = (0, 0, z) (x, y, z) = (0, 0, z) (ξ, η, ζ) = (0, 0, ζ) (ξ, η, ζ) = (0, 0, ζ)
Sphere
x2 + y2 + z2 (λx, λy, λz) λ = (x, y, z) = (0, 0, 0) max{|x|, |y|, |z|} f(x, y, z) = (0, 0, 0) (x, y, z) = (0, 0, 0) max{|ξ|, |η|, |ζ|} (λξ, λη, λζ), λ = (ξ, η, ζ) = (0, 0, 0) ξ 2 + η2 + ζ 2 g(ξ, η, ζ) = (0, 0, 0) (ξ, η, ζ) = (0, 0, 0)
5.3. Algorithm of extended CBC substitution As shown in Fig. 9, the dynamic mesh adjustment of a model with curved primitives includes projective transformations (mapping and inverse mapping), initial mesh alignment, and dynamic mesh adjustment. The processing algorithm of the extended approach is divided into three steps accordingly and described in a C-like language in detail. Step 1. Establishing a projective space and mapping the original model to the projective space. The procedure in Step 1 is to transform all curved primitives from R3 to RP3 through well-defined mapping functions. Initial meshes are generated to the projected primitives (images) before meshes are aligned across borders in the next step. for ∀pr(pr ∈ R3 ∧ pr ∈ {Curved Primitive Set}) do Mid ⇐ select mapping function(pr); map to(Mid , pr, Ipr ); generate initial mesh(Ipr ); pushIpr → RP3 ; for ∀pr(pr ∈ R3 ∧ pr ∈ / {Curved Primitive Set}) do generate initial mesh(pr); duplicate pr → RP3 ; Step 2. Adjusting initial meshes and mapping them inversely to the real space. The initial meshes generated in the previous step are adjusted locally in the projective space. A complete mesh model for analysis is obtained by inverse mapping. Any meshing data changes should be reflected in both spaces.
for ∀Ipr (Ipr ∈ RP3 ) do adjust mesh(Ipr , CBC); refreshIpr → RP3 ; Mid ⇐ select inverse mapping function(Ipr ); inverse map to(Mid , Ipr , pr); refresh pr → R3 ; Step 3. Dynamic mesh adjustment between the two spaces. The relative positions among model components may vary with time during dynamic FEA or simulation. It is necessary to modify the FEM mesh model dynamically whenever the relative positions are changed. The dynamic mesh adjustment is described as follows. Whenever (relative motion occured==true) for ∀pr(pr ∈ R3 ) do status ⇐ check motion status(pr); if (status==moved) Ipr ⇐ select projected object(pr); adjust mesh;(Ipr , CBC) refresh Ipr →RP3 ; Mid ⇐ select inverse mapping function(Ipr ); inverse map to(Mid , pr, Ipr ); refresh pr→R3 ; Special attentions to those nodes on the border zones of the primitives should be given to the inverse mapping in Steps 2 and 3. Fig. 11 illustrates the nodal connectivity before and after the inverse mapping. Since the nodes indicated by 䊊 differ in nodal connectivity after the inverse mapping, each of them are divided into two and assigned to the primitives that previously shared the node. The other nodes given by remain unchanged. Whether to divide a node or not is
L. Wang, T. Moriwaki / Precision Engineering 27 (2003) 245–257
255
three nodes approach to collinear, the accuracy of FEA may decrease. To prevent it from happening, we use 20-node elements instead of 8-node elements for FEA calculations. The details of how to avoid this problem and what techniques to be applied will be reported separately. 6. Case study The proposed CBC substitution is developed for dynamic analysis using finite element method (FEM). The FEM is known as a general technique for constructing an approximation of solution over a collection of finite elements, by using variation concept. When associated with the CBC substitution while finite elements are changing, one problem to be solved is how to handle the nodal analytical data for those adjusted and newly generated nodes. Since the thermal behavior of machines is more influential than most of the other factors [41], a method for dynamic thermal analysis is considered in this case study. The basic procedures are: Fig. 11. Node division during inverse mapping.
simple—it is judged by the nodal connectivity of the sharing primitives in the real object space R3 . After inverse mapping as shown in Fig. 11, there are cases that three nodes of an element share one arc segment. If the
1. Calculating temperature distribution at time t, and recording all nodal analytical data including nodal temperature Tt and its incremental rate T˙ t . 2. Adjusting nodal connectivity when the positions of some elements are changed during a time interval of (t, and recalculating those nodal coordinates accordingly. 3. Deriving and interpolating new nodal temperature for those affected nodes based on Tt and T˙ t .
Fig. 12. Temperature distribution and deformation with moving heat source.
256
L. Wang, T. Moriwaki / Precision Engineering 27 (2003) 245–257
4. Calculating temperature distribution Tt +(t and T˙ t+(t at time t + (t. The above procedures are repeated until the thermal analysis is completed. However, Steps 2 and 3 are needed only if the dynamic mesh adjustment is involved in that cycle. An interpolating algorithm [42] consisting of a primary interpolation and a secondary interpolation is used in Step 3. Generally, there exists a complex phenomenon across the interface between two movable components. How to handle the problem is considered critical to the accuracy of a finite element analysis. Instead of using a specific type of finite elements, we developed a concept called virtual interface [43], across which a virtual heat conduction is taken place between the contacting surfaces to simulate the heat conduction in reality. Frictional heat between relatively moving components is also taken into consideration, as it may serve as an additional heat source. With this virtual interface introduced to accumulate the complex behavior across components, our meshing algorithm becomes generic and simple. A good convergence in this case study is achieved manually by changing step intervals between calculations. Fig. 12 shows nine typical analytical results of both the temperature distribution and the thermal deformation of a simple table-base model of machine tools, by utilizing the methods proposed. In this case, a heat source is fixed on the top surface of the table that moves one round-trip over the base in 30 min during the thermal analysis. Although the case study is given to a simple model with only flat surfaces, other complicated geometries can be dealt with in the same way in terms of analytical procedures, because the major problem of dynamic meshing is handled by the extended CBC substitution approach. Further details on this will be reported separately.
7. Conclusions A new approach called CBC substitution is proposed to deal with initial alignments of constructive mesh generation and dynamic mesh adjustments during FEM analysis. A mapping technique is adopted to transform curved geometries to plane ones, so as to handle the complicated geometries in a consistent way. The meshing process is localized. Only the disagreed elements are adjusted dynamically. Common to other meshing algorithms, this approach also has its limitation. For a complex geometry, if unique mapping and inverse mapping functions cannot be found, the geometry must be decomposed and approximated until it can be described by certain mathematic equations. Other than this limitation, the proposed methodology is found effective and efficient in dealing with dynamic meshing situations. It is expected that the new methodology can be applied to dynamic analyses and simulations for better product design or more accurate device control based on their real life behaviors captured by the simulations.
References [1] Liu WK, Jun S, Li S, Adee J, Belytschko T. Reproducing kernel particle methods for structural dynamics. Int J Numeric Methods Eng 1995;28:1655–79. [2] Lu YY, Belytschko T, Gu L. A new implementation of the element free Galerkin method. Comput Methods Appl Mech Eng 1994;113:397– 414. [3] Oñate E, Idelsohn S. A mesh-free finite point method for advectivediffusive transport and fluid flow problems. Computat Mech 1998;21:283–92. [4] Liszka T, Orikisz J. Finite difference method at arbitrary irregular grids and its application in applied mechanics. Comput Struct 1980;11:83– 95. [5] Mukherjee YX, Mukherjee S. The boundary node method for potential problems. Int J Numeric Methods Eng 1997;40:797–815. [6] Duarte CA, Oden JT. An h-p adaptive method using clouds. Comput Methods Appl Mech Eng 1996;139:237–62. [7] Zhu T, Zhang J, Atluri SN. A meshless local boundary integral equation (LBIE) method for solving nonlinear problems. Computat Mech 1998;22:174–86. [8] Aluru NR, Li G. Finite cloud method: a true meshless technique based on a fixed reproducing kernel approximation. Int J Numeric Methods Eng 2001;50:2373–410. [9] Moriwaki T, Sugimura N, Wang L. A modelling system for finite element analysis of machine products. Trans North Am Manufact Res Inst SME 1993;21:383–9. [10] International Meshing Roundtable. http://www.imt.sandia.gov/. [11] Cavendish JC. Automatic triangulation of arbitrary planar domains for the finite element method. Int J Numeric Methods Eng 1974;8:679–96. [12] Sadek EA. A scheme for the automatic generation of triangular finite elements. Int J Numeric Methods Eng 1980;15:1813–22. [13] Liu Y, Chen K. A two-dimensional mesh generator for variable order triangular and rectangular elements. Comput Struct 1988;29(6):1033– 53. [14] Sarrate J, Huerta A. Efficient unstructured quadrilateral mesh generation. Int J Numeric Methods Eng 2000;49:1327–50. [15] Van Phai N. Automatic mesh generation with tetrahedron elements. Int J Numeric Methods Eng 1982;18:273–89. [16] Cavendish JC, Field DA, Frey WH. An approach to automatic three-dimensional finite element mesh generation. Int J Numeric Methods Eng 1985;21:329–47. [17] Mochizuki T, Kojima T. A shape data processing system for 3-dimensional finite element analysis using solid model. J Jpn Soc Precision Eng 1986;52(11):1910–5. [18] Calvo NA, Idelsohn SR. All-hexahedral element meshing: generation of the dual mesh by recurrent subdivision. Comput Methods Appl Mech Eng 2000;182:371–8. [19] Lo SH. Finite element mesh generation over curved surfaces. Comput Struct 1988;29(5):731–42. [20] Cheng F, Jaromczyk JW, Lin JR, Chang SS, Lu JY. A parallel mesh generation algorithm based on the vertex label assignment scheme. Int J Numeric Methods Eng 1989;28:1429–48. [21] Karamete BK, Beall MW, Shephard MS. Triangulation of arbitrary polyhedra to support automatic mesh generators. Int J Numeric Methods Eng 2000;49:167–91. [22] Dey S, O’Bara RM, Shephard MS. Towards curvilinear meshing in 3D: the case of quadratic simplices. Comput Aided Design 2001;33:199– 209. [23] Owen SJ, Saigal S. H-Morph an indirect approach to advancing front hex meshing. Int J Numeric Methods Eng 2000;49:289–312. [24] Owen SJ. Hex-dominant mesh generation using 3D constrained triangulation. Comput Aided Design 2001;33:211–20. [25] Owen SJ, Staten ML, Canann SA, Saigal S. Q-Morph: an indirect approach to advancing front quad meshing. Int J Numeric Methods Eng 1999;44:1314–40.
L. Wang, T. Moriwaki / Precision Engineering 27 (2003) 245–257 [26] Lu Y, Gadh R, Tautges TJ. Feature based hex meshing methodology: feature recognition and volume decomposition. Comput Aided Design 2001;33:221–32. [27] Li XY, Teng SH, Üngör A. Simultaneous refinement and coarsening for adaptive meshing. Eng Comput 1999;15:280–91. [28] Halpbern M. Industrial requirements and practices in finite element meshing: a survey of trends. In: Proceedings of 6th International Meshing Roundtable, SAND97-2399, Sandia National Laboratories, 1997. [29] Sheffer A, Bercovier M. Hexahedral meshing of non-linear volumes using Voronoi faces and edges. Int J Numeric Methods Eng 2000;49:329– 51. [30] Lai M, Benzley S, White D. Automated hexahedral mesh generation by generalized multiple source to multiple target sweeping. Int J Numeric Methods Eng 2000;49:261–75. [31] Staten ML, Canann SA, Owen SJ. BMSweep: locating interior nodes during sweeping. Eng Comput 1999;15:212–8. [32] Knupp P. Applications of mesh smoothing: copy, morph, and sweep on unstructured quadrilateral meshes. Int J Numeric Methods Eng 1999;45:37–45. [33] Tautges TJ. The generation of hexahedral meshes for assembly geometry: survey and progress. Int J Numeric Methods Eng 2001;50:2617–42. [34] Dhondt G. Unstructured 20-node brick element meshing. Comput Aided Design 2001;33:233–49.
257
[35] Dhondt G. A new automatic hexahedral mesher based on cutting. Int J Numeric Methods Eng 2001;50:2109–26. [36] Tautges T, Blacker T, Mitchell S. The Whisker Weaving algorithm: a connectivity-based method for constructing all-hexahedral finite element meshes. Int J Numeric Methods Eng 1996;39:3327–49. [37] Folwell NT, Mitchell SA. Reliable Whisker Weaving via Curve Contraction. Eng Comput 1999;15:292–302. [38] Ho-Le K. Finite element mesh generation methods: a review and classification. Comput Aided Design 1988;20(1):27–38. [39] Wang L. Development of an integrated CAD/CAE system for machine tool design. Doctoral Dissertation, Kobe University, Japan, March 1993. [40] Yokota I. From topological geometry to projective geometry. Modern Mathematics Press, 1993. [41] Bryan J. International status of thermal error research. Annals of the CIRP 1990;39(2):645–56. [42] Wang L. An interpolating algorithm for dynamic thermal analysis of machine tools. In: Proceedings of the First International Symposium on Advances in Intelligent Computer Integrated Manufacturing System; 1994. p. 213–8. [43] Moriwaki T, Wang L. Study on thermal behaviour of machine tool elements across movable joints under operating states. In: Proceedings of International Conference on Precision Engineering; 1995. p. 150–4.