Finite Elements in Analysis and Design 12 (1992) 1-1 l
1
Elsevier FINEL 254
Optimized/adapted finite elements for structural shape optimization Srinivas Kodiyalam Solid Mechanics Laboratory, General Electric Corporate R&D Center, Schenectady, N Y 12301, USA
V.N. Parthasarathy General Electric Consulting Services, Albany, NY 12205, USA
Received January 1992 Revised March 1992 Abstract. A methodologyfor geometry-based structural shape optimization with adaptive meshing is developed. This method uses a two-level optimization problem where the first level is focussed on generation of good-quality finite elements/meshes and the second level deals with optimization of the overall structural shape based on responses computed on the adapted mesh. The first-level optimization problem uses the finite element nodal coordinates as design variables whereas the second-level problem uses design-oriented, geometry-based parameters for modifying the structural shape. An error indicator, based on effective stress variations in an element, is used to refine the finite element mesh. It is seen from initial investigation that, with mesh adaptation, the accuracy of the structural responses (displacements and stresses) are increased, resulting in a more accurate evaluation of the design objective and constraints, and therefore resulting in more conservative designs.
Introduction
Shape optimization of three-dimensional structures has been a subject of considerable recent research [1-5]. In these investigations, both finite element based and geometry-based approaches have been used for formulating and controlling the problem. The development of fully automatic mesh generators that are capable of producing a valid finite element mesh in a domain of arbitrary complexity has contributed to the growth of this technology. In addition, the use of approximation concepts during the iterative design process has made shape optimization of fairly large-scale three-dimensional structures possible. The result of a shape optimization problem at a given iteration can be adversely or positively influenced by the quality of the structural responses available. The structural responses in turn are governed by the quality of the underlying mesh. When a shape change is introduced to a set of boundary nodes, undesirable element shape distortions may result. Re-analysis on such distorted meshes may lead to poor evaluation of structural quantities. Examples of such cases and their remedy in two-dimensional problems have been reported in Ref. [6]. Furthermore, during the stages of shape design, a newly derived shape may exhibit high stress gradients (or high gradients of any response quantity) in different subregions of the domain, when compared with its previous shape. Simple re-meshing of the geometry may not yield responses of desired accuracy for such cases. Correspondence to: Dr. Srinivas Kodiyalam, General Electric Corporate R&D Center, Building K-I, Room 2A25,
P.O. Box 8, Schenectady, NY 12301, USA. 0168-874X/92/$05.00 © 1992 - Elsevier Science Publishers B.V. All rights reserved
2
S. Kodiyalam, [AN. Parthasarathy / Optimi:ed/adaptedf~nite elements
Element distortions in a mesh can independently bc reduced to a great degree by taking the approach described in Ref. [7] or its variation. Such geometric-adaptive mesh improvements use mainly the element shape information while retaining the mesh's original element density. Nevertheless, such techniques have been found to aid in alleviating the ill-conditioning of the finite element solution procedure and hence obtain a physically meaningful numerical solution for a given mesh. Solution-adaptive methods on the otherhand, increase or decrease a finite element mesh's element density in accordance with the error in the numerical solution. Solution-adaptive techniques have been successfully implemented in various studies involving stand-alone linear/non-linear structural problems. The advantage of solution adaptive mesh-adaptation is in the increased reliability of the final numerical solution to the particular problem of interest. In the context of three-dimensional shape optimization, there is little or no research into the use of mesh adaptation techniques to obtain more accurate structural responses. In this paper, results of some preliminary investigation from an integrated shape optimization system with adaptive re-meshing and analysis are presented. Mesh adaptation in this context refers primarily to the solution-adaptive method, while geometric-adaptive mesh improvements are also used, where appropriate. A two-level optimization strategy is used where the first-level problem deals with optimization for improvement of e l e m e n t / m e s h quality and the secondlevel problem deals with optimization of overall structural shape based on response computed on the adapted mesh. The implementation of a solid modeling based 3-D shape optimization system, reported in Ref. [4], forms the basis of this work.
Importance of mesh quality in analysis and its control Research in the theory of finite element analysis has placed increased importance on the size and shape of the elements generated by a mesh generator. The quality of a finite element solution, and hence the final design, has been shown to be affected by the quality of the underlying mesh. A poor mesh may lead to unstable a n d / o r inaccurate finite element approximations. This stringent requirement on the size and shape has been a contributing factor in the selection of a particular mesh generator for a given analysis application. Therefore, it is imperative that the quality of a mesh--initial or a d a p t e d - - c a n be monitored and controlled in an automatic manner [8]. Literature reveals one often used technique for mesh quality i m p r o v e m e n t - - t h e Laplacian smoothing. This technique iteratively moves every node to the centroid of the surrounding polygon. The primary disadvantages of Laplacian smoothing are: (i) the inherent insensitivity to the shape of any element incident at the node being smoothed; therefore, this procedure potentially can deteriorate mesh quality in one region, while improving in other regions of the computational domain; (ii) the potential movement of the nodes to the exterior of the domain in the presence of non-convex local geometry; and (iii) the positioning of a node inside the volume of a neighboring element in 3-D meshes. Thus, Laplacian smoothing c a n render a mesh unsuitable for analysis. Several other technique developed to avoid these problems are derivatives of Laplacian smoothing. They include: (i) area-weighted Laplacian smoothing, where a node is moved toward larger area elements in the patch; and (ii) length-weighted Laplacian smoothing, where the node is moved toward the larger contributing peripheral node. Although these modifications have been somewhat successful in ameliorating the problems, they do not guarantee success for any arbitrary mesh. Such a guarantee is essential whenever the mesh generator is under program control of an automated procedure, such as shape optimization.
S. Kodiyalam, V.N. Parthasarathy / Optimized / adapted finite elements
Laplace-smoothed surface mesh 3 inverted triangles
3
Result of optimizing Laplace-smoothed mesh on left - 0 inverted triangles
a) Enlarged view of the surface mesh patch marked by the dotted ellipse in figure (b)
Fig. 1. Exampleof first-levelmesh optimization. Clearly, for guaranteed success, the solution lies in restricting a node's movement, both in distance and direction. Constrained movement of a node requires evaluation of many geometric quantities and topological information of the nodes and elements connected to the node being smoothed. In terms of quality and cost trade-offs, this translates to the fact that one has to pay an increased computational expense for an "intelligent" mesh quality improvement algorithm. A geometric-adaptive mesh improvement technique based on constrained numerical optimization for triangular and tetrahedral meshes has been described in Ref. [7]. In the present work, this technique forms the first-level optimization problem, stated as: Find the set of design variables, D, that will: Minimize: Subject to:
F ( D ) = Finite element's shape distortion ratio Volume of each individual tetrahedral element >/6, d1<~di<~d u
i = 1 . . . . . n,
(1)
where D, the vector of design variables, corresponds to the finite element nodal coordinates. This optimization problem is implemented and solved within the automatic mesh generation process. The high computational costs for tetrahedral meshes was found to be mainly due to the elaborate geometric computations required to evaluate the standard shape distortion measure in Ref. [7]. Use of a new element shape distortion measure has reduced this computational cost by a factor of more than two [9].
Solution-adaptive m e s h i n g and analysis
Numerical solution of structural problems have always been addressed in the light of the conflicting requirements of accuracy, reliability, and efficiency. Whereas accuracy is mainly concerned with the fundamental mathematical formulation of the physical phenomenon,
reliability and efficiency of the solution depend on the techniques of numerical approximation of the said mathematical statement. Conventionally, these requirements have been satisfied for a problem by starting with an initial coarse grid analysis and subsequently performing many analyses on finer meshes. For each such analysis. the acceptance of the mesh and the resulting solution was based mainly based on the user’s previous experience. Since shape optimization is an iterative process, arriving at reliable numerical solutions has to bc performed non-interactively. Efforts in the area of finite element error estimation and its control have been summarized in Ref. [lo]. Once a reliable error estimator has been developed for a given mathematical problem and its numerical approximation, it can be used to adaptively modify the mesh in such a way as the reduce/re-distribute the error to an acceptable level. If this process can be performed in an automated manner, the results of the scheme can be used as input to the shape optimization procedure and hence avoid the mesh distortion problems discussed earlier. This realization is the motivation for the present work. Error estimation for computational mechanics problems have been performed in many ways. An error norm based on strain energy of the solution was formulated in [IO]. A method based on projected nodal stress and the actual finite element nodal stress was proposed in [ll]. Yet another adaptive method based on the strain energy density variations was introduced in [12]. An error norm based on force equilibrium satisfaction of the solution was used to adaptively solve linear/ non-linear problems in [ 1.31.An estimator based on stress gradients was used in [14]. Researchers in computational fluid dynamics use less rigorously derived indicators-usually based on gradients and/or curvature of a particular solution quantity of interest [15,16]. In the present work, a variation of the implementation of Ref. [ 141 that has been extended to 3-D problems is used. The basic premise of the current implementation can be stated as follows. If the overall error, E, defined as the difference of error indicators amongst all elements in the mesh approaches zero, then the mesh is considered to have zero solution error. The particular element error indicator used in the current work is defined as:
where ue is the effective finite elements stress at the nodes of element i. Mathematically, the above definition of error gives the maximum nodal stress variation in an element. If it is required that the maximum variation be less than a prescribed upper limit, then elements with large variation can be refined in a consistent manner to capture high local solution gradients. In the limit, the use of this indicator will yield a uniform distribution of solution error in the whole computational domain. This formulation is also consistent in a larger sense with the work in Ref. [ll] which has successfully used the point-wise measure of error defined as: EP =
u,~-
(3)
&,,
nodal stress and U, is the finite element nodal stress at node i. where u,~ is the projected Several error norms for the entire mesh can now be defined using the elemental error values. The average solution error for the whole domain is calculated as follows: ,= Ea,g =
where
111
C
E, is the error
Ei/m, in element
i, determined
using eqn. (2) and m is the total number
of
S. Kodiyalam, V.N. Parthasarathy / Optimized / adapted finite elements
5
elements in the mesh. The average variation in local error can be found for the entire domain using: i=rn I
elocaI = ~
Ei/ml,
e i > 2Emin and Emin
~
UXS(Enen),
(5)
i=1
where nen is the number of edge neighbors, e i is the error in element i, determined using eqs. (2) and m~ ~ m. In the limit, a uniform distribution of solution error is achieved when across adaptations the ratio -q = E a v g / E l o c a I approaches unity. The principle of uniform distribution of solution error for mesh adaptation has been successfully employed to obtain more accurate results [11,15].
Implementation of mesh adaptation Mesh adaptation/modification has conventionally been approached in one of the following ways: (1) r-method: the existing nodes in a mesh are clustered into regions of high solution activity, before conducting a re-analysis. This is done by first associating a weight and direction of node movement in relation to the error estimate and subsequently performing a re-positioning operation of all the nodes. (2) h-method: new, smaller-sized elements are created by injecting points into existing elements having a high error estimate and vice versa. This method is the most commonly used due to the relative ease of re-meshing. (3) p-method: the mesh is generated once and the elements and nodal locations remain invariant across adaptations, but the degree of the approximating shape function is increased for those elements whose error estimate is above the preset threshold, before conducting a re-analysis. The h-method was chosen for the adaptive scheme described in this paper, as it does not require formulation of new element types as in the p-method and it also avoids the unnecessary element shape distortion of the r-method. The implementation of an h-adaptive scheme is largely dictated by the capabilities and limitations of the mesh generator. Accordingly, several modalities can be found in the literature. The present implementation of the adaptive scheme works with a octree-based fully automatic mesh generator, OCTREE [17]. The finite element software, ANSYS, is used for analysis and is executed by forking a process from the DECstation to the C O N V E X mini supercomputer. After an initial finite element analysis has been conducted, the cut-off error value for refinement is computed from eqn. (2). Next, elements in the mesh whose error exceed the cut-off value are slated for refinement. The refinement is achieved by first creating a new finite element node at the centroid of the element and subsequently assigning a mesh control value which is 50% less than the currently set value. In essence, the octree decomposition of the domain is further locally decomposed to octants of smaller size. When all elements in the mesh have been processed for refinement, the model is re-meshed using the mesh generator. During the re-meshing stage, the previously created vertices come into play by driving the octree subdivision near its vicinity. The analysis is then carried out on the new mesh. The process of refinement and re-analysis continues until the desired convergence of the ratio of average errors (-q) between consecutive meshes is achieved. The algorithm is described schematically in Fig. 2. As designs evolve, some constraint vertices created as part of mesh adaptation for a previous design may fall outside the domain of the new design. Initial processing has been enhanced to remove such outside vertices before the start of the analysis for new design.
6
S. Kodiyalam, 1~N. Parthasarathy / Optimized / adapted finite elements
Shape optimization problem The work reported in Ref. [4] adopts a hybrid approach to the problem of 3-D shape optimization. It embraces the notion of high-level problem definition and shape control via constructive solid geometry, while at the same time exploiting the benefits of boundary representation for specifying the physics of the problem and for meshing the part. This hybrid approach is strongly coupled to an automatic mesh generator and uses to advantage the explicit association of the finite element data with the model geometry for performing shape sensitivity analysis. In order to avoid discontinuities in the sensitivity coefficients, no remeshing is performed during the sensitivity calculations instead mesh relaxation procedures are used to update the geometry to reflect design variable perturbations. Also, certain hybrid approximation methods are used to reduce the number of detailed finite element analyses. For completeness, the shape optimization problem statement and the design model in outlined below. The problem defined in eqn. (6) forms the second level of optimization: Find the set of design variables that: Minimize:
F(X)
Subject to:
g/( X ) <~0
(objective function), (inequality constraint), (side constraints),
(6)
where X is the vector of design variables. In the shape design problem, the design variables are geometry-based parameters that control the sizes of the geometric primitives that comprise an object and details of two-dimensional profiles that are used to generate swept and lofted 3-D objects. The objective function is the volume of the structure, and structural response constraints are imposed on displacements and stresses at various locations of the object. A geometry-based specification of the response constraints is adopted, whereby displacements can be constrained on a particular edge, stresses on a face of the model, etc. Structural response constraints can also be imposed on the whole region. Geometric requirements, such as preventing surface intersections to maintain a topologically valid model, can be treated in the design problem. Side constraints are also imposed on the design problem.
Geometry-Based I InputProblemS~ecifications ~"
OCTREE t ANSYS
COMPUTE ERROR INDICATOR
Exit Mesh [ Adaotation ~
No
ADD MESH REFINEMENT SPECIFICATION
Fig. 2. Flowdiagramof mesh adaptation.
S. Kodiyalarn, V.N. Parthasarathy / Optimized / adapted finite elements
The integrated shape optimization system architecture An enhanced version of the shape optimization system originally developed in Ref. [4] is shown in Fig. 3. Detailed descriptions of the different aspects of this system can be found in Ref. [4]. In this system, SDRC'S GEOMOD is used as the solid modeler for creating and modifying the geometry of the 3-D object. Since a procedural call cannot be made to execute GEOMOD, the newly developed interprocess communications capability in ~-DEAS is used to establish a link between GEOMODand the 3-D shape optimization modules. The geometry and the topology data of the 3-D object is extracted from the PEARL relational database using a boundary representation translator and written into an archive file format representation suitable for the geometric modeling utility, TAGUS [18]. Thus, TA6US provides an interaction mechanism with the geometric representation derived from the modeler. The adaptive scheme, shown in Fig. 2, is performed and the adapted mesh is used for computing the design sensitivity coefficients of the objective and constraint functions with respect to the design variables. For sensitivity analysis, the mesh topology is fixed and no new mesh is generated. The link between the DSA routine and GEOMOD is through UNIX pipes, using the SDRC/I-DEAS director-observer capability, in conjunction with their program files. The primitive's sizes, and details of 2-D profiles used to create the 3-D model, are updated in
CONVEX r ............
DECstation
................................................. SDRC, Ul~ [;.aPE - - -O P T I M ~
GEOMOD
I
CONTROLLER ]
I 3 . e a ~ - . q ~ I . . . . . . ~. . . . . . . . . . . . . . . .
~roblem Description
TAGUS
LOADS/ B. C. 's .................
GEOMETRY CONSTRAINTS
o~--°-°--o~°°--°°J
os sj
E E. ANALYSIS
I & [:f MESH ADA'I~ATi0N:ii:::]
.
DESIGN SENSITMTY ANALYSIS (DSA) I UNIX
/
OPTIMIZERI (ADS)
]
APPROXIMATE ANALYSIS
RO.SS, Not executed during Design Sensitivity Analysis~DSA)
Fig. 3. Shape optimization system architecture.
8
S. Kodiyalam, V.N. Parthasarathy / Optimized/ adapted finite elements
the program file to reflect the perturbations in the design variables, and GLOMOD is then re-executed, The response computed by the analysis program and the gradient information from DSA are used to generate a hybrid approximation to the optimization problem defined in (1). This is done in the approximate analysis routine. The numerical optimizer, ADS [19], is used to solve the approximate optimization problem. During the optimization stage, no finite element analysis is performed, instead, mathematical approximations to the objective and constraint functions are used [4]. After convergence of the approximate optimization problem, the whole design cycle is repeated until a satisfactory convergence of the original optimization problem is obtained.
Design example: Aircraft turbine disk The aircraft turbine disk, shown in Fig. 4, was originally solved for optimal shape without mesh adaptation in Ref. [4]. Invoking the variational geometry capability with >DEAS, the 2-D cross-section of the disc was described and revolved by a prescribed angle to create the 3-D model. Appropriate symmetry boundary conditions were enforced on the front and the back faces of the disk. The disk was subjected to centrifugal loading and a pressure load on the rim to simulate the forces induced by the blades. The material properties and loading are as follows: = 27.1 × 106 psi, Young's modulus Poisson's ratio = (I.3, Density -- 0.284 lb/in -~, Rotational speed = 11,500 rpm = 1204 r a d / s , Blade load = 3869.4 lb/deg, Allowable Von Mises stress = 114,000 psi, Allowable tangential stress (rim) = 93,000 psi. The design problem was to minimize the volume of the disk with constraints imposed on the Von Mises effective stresses over the whole region (rim, web and bore) and the tangential stress at the rim. The design model had a total of 5 design variables that were details of the profile used to generate the swept model. More specifically, the design variables were geometric locations, on the bore and the web of the disk, that were inturn used to define the 2-D cross-section. Table 1 provides the initial and final values and bounds on these design variables. The initial design had an volume of 3.546 in 3 with all the constraints satisfied. The finite element mesh of the initial design geometry, shown in Fig. 4, consisted of 144 parabolic tetrahedron elements. Mesh adaptation was carried out, until a convergence in the error criterion was obtained. The average to local ratio (rD iteration history is shown in Fig, 5. The adapted mesh, corresponding to the initial design, shown in Fig. 4, had a total of 350
Table 1 Design variable values for aircraft disk problem Design variable
Initial value (in)
Final value (in)
[x)wer b o u n d fin)
U p p e r bound (in)
1 2 3 4 5
1.10 0.75 0.60 0.60 0.60
0.496 0.247 0.226 O. 172 0.200
O. I O. 1 1). 1 O. 1 O. 1
10.0 10.0 10.0 10.0 10.0
S. Kodiyalam, V.N. Parthasarathy / Optimized / adapted finite elements
9
Blade ing Bore-
~
W
e
b
~-
~
~
~im~Fig. 4. Aircraft turbine disk--Pre- and post-adapted mesh for initial design.
elements. The responses obtained from the analysis on the adapted mesh was then used for design sensitivity and optimization. A lighter design, with a volume of 1.61 cu.in., obtained at the end of 4 design iterations, is shown in Figure 6. The final design obtained here is a heavier design compared to the one reported in Reference [4]. This increase in objective is attributed to the fact that with adaptive meshing a more precise estimate of the structural responses (displacements and stresses) can be obtained. The objective function and maximum constraint iteration histories are shown in Figures 7 and 8, respectively and Table 1 provides the initial and final values of the design variables with their respective bounds.
Conclusions
A geometry-based structural shape optimization methodology with adaptive meshing is developed. The emphasis of this paper is on the need for: (i) generation of good quality meshes; and (ii) adaptive meshing within the automated shape design process. A two-level optimization problem is used, where at the first level nodal coordinates are used as design variables to generate good-quality finite elements/meshes and design-oriented, geometrybased parameters are used at the second level to modify the structural shape. An error
1.5 ¸
1.4
Q
1.3
1.2
1.1 2
3
4
Iteration number Fig. 5. Error ratio (rt) history from adaptive meshing (for initial geometry).
ZIJ~X Fig. 6. Optimized shape of aircraft turbine disk.
S. Kodiyalam, KN. Parthasarathy / Optimized / adapted,finite elements
10
15r I
35
,z
1]
i
3 >
25
os] z
I
O
E
Ol-
2 1.5
8
-
i
4
~~
05
0.5 1
2
3
4
Iteration n u m b e r Fig. 7. Volume iteration history for aircraft disk.
5
15
0
1
2
3
Iteration n u m b e r Fig. 8. Constraint iteration history for aircraft disk.
indicator, based on effective stress variations in an element, is used to refine the finite element mesh. It is seen from initial investigation that, with mesh adaptation, the accuracy of the structural responses (displacements and stresses) are increased resulting in a more accurate evaluation of the constraints. Increased accuracy results in more conservative designs as compared to non-adapted designs. Future research will address the convergence issues
associated with the variable nature of mesh refinement results in the design optimization loop.
References [1] R.J. YANG and M.E. BOTKIN, "A modular approach for three-dimensional shape optimization of structures", AIAA J. 25 (3), pp. 492-497, 1987. [2] S. KODIYALAM and G.N. VANDERPLAATS, "Shape optimization of three-dimensional continuum structures via force approximation techniques", AL4A J. 27 (9), pp. 1256-1263, 1989. [3] S. KODIYALAM, G.N. VANDERPLAATS and H. MIURA, "Structural shape optimization with MSC/NASTRAN", Comput. Struct. 40 (4), pp. 821-829, 199l. [4] S. KODIYALAM, V. KUMAR and P.M. FINNIGAN, "A constructive solid geometry approach to three-dimensional shape optimization", A/AA J. 30 (5), pp. 1408-1415, 1992. [5] M . E . BOTKIN, "Shape optimization using fully-automatic 3-D mesh generation", Proc. A I A A / A S M E / A S C L / A H S / A S C 32nd Structures, Structural Dynamics and Materials Conf, AIAA, Washington, DC, pp. 815-823, 1991. [6] h . KIKu(ln, K.Y. CHUNG, T. TORIGAKI and J.E. TAYLOR, "Adaptive finite element methods for shape optimization of linearly elastic structures", Comput. Methods Appl. Mech. Eng. 57, pp. 67-89, 1986. [7] V.N. PARTHASARA'rHY and S. KODIYALAM, "A constrained optimization approach to finite element mesh smoothing", Finite Elements in Analysis and Design 9 (4), pp. 309-320, September 1991. [8] T. MITT¥, T. BAKER and A. JAMESON, "Generation and adaptive alteration of Unstructured tree dimensional meshes", Proc. 3rd Int. Conf. & Numerical Grid Generation for Computational Fluid Mechanics and Related Fields, edited by A.S. ARC'UJ.A, J. HAUSER, P.R. EElSEMANN and J.F. THO~4PSON, North-Holland, Amsterdam. 1991. [9] V.N. PARTHASARATIIY, C.M. GRAICHEN and A.F. HATHAWAY, "Fast evaluation and improvement of 3-13 tetrahedral grid quality", NASA Conf. Proc. on Software Systems for Surface Modeling and Grid Generation, edited by R.E. SMITH, NASA CP-2166, NASA Langley Research Center, Hampton, VA, April 1992. [10] 1. BABUSHKAAND W.C. RHEINBOLDT, "Error estimates for adaptive finite element computations", SlAM .L Numer. Anal. 15 (4), 1978. [11] O.C. ZIENK1EWICZ and J.Z. ZHU, "A simple error estimator and adaptive procedure for practical engineering analysis", Int. J. Numer. Methods Eng. 24, pp. 337-357, 1987.
S. Kodiyalam, V.N. Parthasarathy / Optimized / adapted finite elements
11
[12] M.E. BOTKIN, "All adaptive finite element technique for plate structures", Proc. A I A A / A S M E / A S C E / A H S / ASC 24th Structures, Structural Dynamics and Materials Conf., AIAA, Washington, DC, pp. 641-648, 1983. [13] T.D. BLACKER, J. JUNG and W.R. WITKOWSKI, "An adaptive finite element technique using element equilibrium and paving", Presented at the ASME Winter Annu. Meet., ASME Paper 90-WA/CIE-2, ASME, New York, 1990. [14] J-H., CHENG,P.M. FINNIGAN,A.F. HATHAWAY,A. KELA and W.J. SCHROEDER,"Quadtree/octree meshing with adaptive analysis", Proc. 2nd Int. Conf. in Numerical Grid Generation for Computational Fluid Mechanics, edited by S. SENGUPTAet al., Pineridge Press, Swansea, UK, 1988. [15] V.N. PARTHASARATHY and S. SENGUPTA,"Solution contour based adaptive re-gridding", Proc. 3rd Int. Conf. in Numerical Grid Generation for Computational Fluid Mechanics and Related Fields, edited by A.S. ARCILLA, J. HAUSER, P.R. EISEMANN and J.F. THOMPSON, North-Holland, Amsterdam, 1991. [16] K. MORGAN, J. PERAIRE, J. PE1RO and O. HASSAN, "The computation of 3-dimensional flows using unstructured grids", Comput. Methods Appl. Mech. Eng. 87, pp. 335-352, 1991. [17] C. GRAICHEN, A.F. HATHAWAY,P.M. FINNIGAN, A. KEAL and W.J. SCHROEDER, "A 3-D fully automatic geometry-based meshing system", Presented at ASME WinterAnnu. Meet., ASME Paper 89-WA/CIE-4, ASME, New York, 1989. [18] P.M. FINNIGAN, A. KELA and J.E. DAVIS, "Geometry as a basis for finite element automation", Eng. Comput. 5, pp. 147-160,1989. [19] G.N. VANDERPLAATS,ADS User's Manual, Version 3.0, VMA Engineering, Goleta, CA, March 1988.