The Influence of Geometric Approximation on the Accuracy of High Order Methods Xiaojuan Luo Mark S. Shephard Jean-Francois Remacle Scientific Computation Research Center Rensselaer Polytechnic Institute Troy, NY 12180
[email protected]
Abstract The importance of using properly curved mesh entities when using high order discretizations to solve partial differential equations over curved domains is examined. The results show that sizable errors are introduced into the prediction of parameters of engineering interest when the geometric approximation is too low with respect to the polynomial order of the discretization.
1. Introduction High-order (p-version) finite element methods are capable of exponential rates of convergence [7]. This increased approximation power allows the use of coarse meshes in which relatively few elements cover large portions of the computational domain. However, preserving the convergence rate for problems over curved domains requires that a curvilinear mesh geometry representation be used [3]. The purpose of this paper is provide a simple demonstration of the importance of employing appropriate geometric approximations to maintain the accuracy of solution parameters in p-version finite element methods. In today’s engineering design processes the geometric domains are typically defined within the solid modeler of the CAD system. The technologies to support the automatic generation of unstructured, straight-edged, planar faced simplex meshes for domains defined in solid modelers are available. A number of these capabilities directly access the boundary representation and the functionality of the solid modeler to support the geometric operations needed by the mesh generator (e.g. [5]). These
functional interfaces can also be used to support the generation of curved meshes which, as demonstrated by example in this paper, are important for use with p-version finite elements. Section 2 briefly discusses the issues associated with the representation and creation of curved mesh entities. To demonstrate the importance of using a proper geometric representation, Section 3 considers the influence of different geometric approximations in p-version finite element ∞ analyses on the accuracy of an L norm of engineering significance for a simple curved domain problem.
2. Representation of Curved Mesh Entities Current procedures for creating a curved mesh begin from a straight sided mesh and curve the mesh edges and faces to better approximate the geometry. Figure 1a shows a straight sided mesh for a hollow cylinder while Figure 1b shows the same mesh after curving. In this example the mesh edges and faces are quadratic in shape, including curved mesh entities on the interior of the mesh as needed to ensure the validity of the resulting elements. In addition to the possibility of curving interior mesh entities, mesh modification operations are often applied to ensure the validity of the resulting curved mesh can be maintained [2].
a) straight-sided mesh
b) curved mesh
Figure 1. Simple example of mesh curving. A key question in the creation of the curved mesh entities is the selection of geometric form of the mesh entities. For those mesh entities that are on the boundary of the model one approach for defining their geometric
form is to assign them the same geometry as that portion of the model. The other approach that can be used for all curved mesh entities is to assign them an appropriate geometric form, which for those mesh entities on curved model boundaries represents an appropriate geometric approximation for that portion of the domain. A general approach to assign the mesh entities on curved boundaries the geometry of that model boundary is as follows: • Assign an appropriate portion of the model boundary parametric space to map to the parametric space of the element. • Employ the kernel of the solid modeler to provide the geometric information needed in the finite element process which is typically limited to pointwise interrogation of geometric derivatives as needed to evaluate the Jacobian mapping at numerical integration points. Although the viability of this has been demonstrated using a geometrybased analysis code linked with commercial solid modelers [3], it was found to be a computationally expensive process, particularly in non-linear problems that would require many re-evaluations of the element integrals. An alternative approach is to use appropriate geometric approximations for the mesh entities that are determined once and stored with the mesh entities. The most obvious option for such geometric approximation is appropriate order Lagrangian interpolants fit through selected points. An effective procedure for creating quadratic Lagrangian’s has been developed and implemented in a procedure that curves initially straight-sided, planar faces meshes [2]. Consideration of the extensions of this approach to higher order Lagrangians indicated that the procedure quickly became computationally demanding and geometrically complex to the point that qualifying how to control the geometric approximations in a predicable manner became infeasible. An alternative appropach under development is to employ more general geometric approximations based on Bezier curves and surfaces which possess a number of properties [4] that are advantageous to the construction of a mesh curving procedure. One advantage of Bezier geometric approximations is the shape is defined by a control polygon and there are alternative criteria that can be used in the construction of the control polygon. For example, in addition to supporting the definition of simple
interpolating mesh entities, it is possible to provide some level of control of the derivatives between mesh entities.
3. Influence of Mesh Entity Geometry on Results Although basic theoretical results do confirm [3] the need to increase the order of geometric approximation as the order of the element polynomial order is increased, the contention is often made that going to a quadratic shape for the mesh geometry is satisfactory for engineering calculations. Since the available theory is unable to specifically quantify the influence of the geometric approximation, a simple numerical study was performed. The basic problem considered is that of an infinite plate with a circular hole subjected to uniform tension in one direction at infinity [9]. To allow the numerical study of this problem a 20x20 square with a circular hole of radius one was considered with Dirichlet (displacement) boundary conditions corresponding to the exact solution applied on the outer boundaries defined as: σ∞ a r a a3 u x = ---------- --- ( κ + 1 ) cos θ + 2 --- ( ( 1 + κ ) cos θ + cos 3θ ) – 2 ----3- cos 3θ 8G a r r σ∞ a r a a3 u y = ---------- --- ( κ – 3 ) sin θ + 2 --- ( ( 1 – κ ) sin θ + sin 3 θ ) – 2 ----3- sin 3 θ 8G a r r E G = -------------------- , κ = 3 – 4ν 2(1 + ν) where E is Young’s modulus and ν Poisons ratio. In the test performed, the above parameters are selected as: 5
σ ∞ = 1MPa , E = 1.0 ×10 MPa , ν = 0.3 One solution parameter examined is the energy norm defined as: u
e(Ω)
1⁄2 1 = --- ( [ D ] { u } ) T [ K ] [ D ] { u } dΩ 2
∫
Ω
where
∂ ∂x [D] =
0 ∂ ∂y ∂ ∂y
0 ∂ ∂x
E -------------21–v Ev , [ K ] = -------------21–v 0
Ev -------------21–v E -------------21–v
0 0 .
E -------------------2(1 + v)
0
The exact value of strain energy is 3.6414711093E-3 (with 13 digits precisions) over the finite domain considered. The relative error in energy norm is defined as, e
e(Ω)
u exact e ( Ω ) – u fem e ( Ω ) = ----------------------------------------------------------exact u e(Ω) ∞
The other solution parameter investigated is an L norm of engineering interest which is the peak stress. For this problem the peak stress occurs at the top of the circular hole and has an exact value of 3.00000MPa. Figure 2a shows the basic straight-sided mesh used in this example while Figure 2b shows the mesh with the edges curved. To ensure the correct implementation of the basic expressions, the problem was solved employing non-uniform h-refinement (see Figure 2c and Figure 2d). The h-refinement produces a relative energy norm error of 1.674% and a peak stress of 3.00472MPa which is in error by 0.16%. To study the influence of the geometric approximation of circular hole on the solution results, this problem was solved using p-version polynomial finite elements of order p = 1 to p = 10 on meshes where the circular hole was approximated as four linear, quadratic, cubic or quartic curved mesh edges. In the quadratic, cubic and quartic curve mesh edge cases consideration was given to meshes where the curves interpolated the correct number of equally spaced points along the arc, and where the curves 1 were constructed to enforce slope continuity ( C ) between the elements. Figure 3 shows the two different forms of geometric approximation for the quadratic and cubic geometric mesh edge approximations. Figure 4 shows the log of relative error in the energy norm versus polynomial order for the different orders of geometric approximation. Consistent with basic theory [3], as the polynomial order of the finite element is
a) straight sided mesh
b) curved mesh
d) close-up of refined mesh c) h-refinement Figure 2. Plate with circular hole. increased past that of the geometric approximation, the error ceases to decrease at the appropriate rate and eventually stops decreasing. The value of this error value decreases with increase in the order of geometric approximation. Table 1 shows the values of stress with respect to the polynomial order for different geometric approximations. The relative error in maximum stress shown in Figure 5 is defined as e
∞
σ fem ∞ – σ exact ∞ = ---------------------------------------------σ exact ∞
1 0.95 0.9 Model Slope continuity Interpolating
Y
0.85 0.8
0.75 0
0.2
0.4
0.6
0.8
a) quadratic geometry approximation 1
0.95
Y
0.9
0.85
Model Slope continuity Interpolating
0.8
0.75
0.7 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
X
b) cubic geometry approximation Figure 3. Geometric approximation to circular hole. In the case of the linear geometry elements the stress is greatly under predicted with a linear displacement element and quickly rises well past the exact value as the displacement polynomial order is increased such that at p = 10 the stress prediction is nearly 3 times the correct value. As shown in Figure 5, the relative error in maximum stress is ± ( 10 ∼ 20 ) % for quadratic geometry approximation, ± ( 0.5 ∼ 3 ) % for cubic geometry approximation and ± ( 0.05 ∼ 0.3 ) % for quartic geometry approximation. In the case of the quadratic element, the form of geometric approximation used strongly influences the results with the interpolating geometry converging toward a value of stress that is at least 15%
22.9938 1
10
2.2789 0
10
−1
q =1 q =2 q =3 q =4
0.4158
10
0.0851
q − Geometry approximation order −2
10
0
2
4 6 Polynomial order (P)
8
10
RELATIVE ERROR IN ENERGY NORM (%)
RELATIVE ERROR IN ENERGY NORM(%)
2
10
2
10
22.9938 1
10
2.0498 q =1 q =2 q =3 q =4
0
10
−1
0.4058
10
0.04297
q − Geometry approximation order −2
10
0
2
4 6 Polynomial order (P)
8
10
a) interpolating geometry b) slope continuous geometry Figure 4. Relative Error in energy norm higher than the exact value and the case of the slope continuous geometry converging to a value that is 10% lower than the exact value. As the geometric approximation is improved to cubic and quartic the errors in the calculated stresses decrease to more acceptable levels. 0
order linear 1 1.30586 2 2.21293 3 3.15026 4 4.01508 5 4.85085 6 5.67879 7 6.49054 8 7.28651 9 8.07482 10 8.85378
1
0
1
C C C C Quadratic Quadratic Cubic Cubic 1.30586 1.30586 1.3058 1.3058 2.46639 2.50827 3.0510 2.9402 2.78196 2.63492 3.0891 3.0718 3.08403 2.81369 3.0982 3.1406 3.17727 2.95646 3.2021 3.2201 3.25421 2.83807 3.2480 3.2997 3.29829 2.66061 3.1607 3.2503 3.37931 2.65638 3.0472 3.1728 3.46451 2.62511 2.9861 3.1154 3.51207 2.63147 2.9841 3.0863 Table 1. Peak stress predictions
0
C Quartic 1.3058 3.0471 3.0953 3.0666 3.1831 3.2274 3.1591 3.0772 3.0253 3.0067
1
C Quartic 1.3058 3.0359 3.0920 3.0696 3.1836 3.2302 3.1648 3.0866 3.0417 3.0025
4. Closing Remarks This paper has considered the influence of the element geometry approximation on p-version finite element results for a simple linear elliptic problem subject to a mild stress concentration factor of three. As
RELATIVE ERROR IN MAXIMUM STRESS(%)
60
q=1 q=2 q=3 q=4
Out of range 40
20
17.06 0.22
0
−0.53 −20
q − geometry approximation order
−40
−60 1
2
3
4 5 6 7 Polynomial order (P)
8
9
10
a) Interpolating approximation
RELATIVE ERROR IN MAXIMUM STRESS
20 17.1
15
Interpolating Slope continuity
10 5 0 −5 −10 −12.2 −15 −20 2
3
4
5 6 7 Polynomial order (P)
8
9
10
b) Quadratic approximation Figure 5. Relative error in maximum stress expected, the results clearly show that the use of linear geometric approx∞ imations lead to large errors with the errors in L norms of engineering interest being so large that the predictions are meaningless. It was also shown that although going to a quadratic approximation does increase the accuracy of the predictions, key parameters of engineering interest can be over predicted by an unacceptable high level. Additional study is needed to qualify the level of geometric approximations needed for problems
with more severe behaviors as well as the influence on the solution of more complex parabolic and hyperbolic equations over curved domains analyzed with high order discretizations.
5. Acknowledgment Aspects of this work were supported by the Office of Naval Research through grant N00014-99-0725, the U.S. Department of Energy Scientific Discovery Through Advanced Computing Program agreement DEFC02-01ER25460, and Simmetrix Inc.
6. References [1] Beall, M.W. and Shephard, M.S., “A General Topology-Based Mesh Data Structure”, Int. J. Num. Meth. Engng., 40(9):1573-1596, 1997. [2] Dey, S., O’Bara, R.M. and Shephard, M.S., “Curvilinear mesh generation in 3D”, Computer-Aided Design, 33:199-209, 2001 [3] Dey, S., Shephard, M.S. and Flaherty, J.E., “Geometry-based issues associated with p-version finite element computations”, Comp. Meth. Appl. Mech. Engng., 150:39-55, 1997. [4] Gerald Farin, Curved and Surfaces for Computer Aided Geometric Design, 1992 Academic press. [5] Shephard, M.S., “Meshing environment for geometry-based analysis”, Int. J. Num. Meth. Engng., 47(1-3):169-190, 2000. [6] Shephard, M.S., Dey, S. and Flaherty, J.E., “A straight forward structure to construct shape functions for variable p-order meshes”, Comp. Meth. Appl. Mech. Engng., 147:209-233, 1997. [7] Szabo, B.A. and Babuska, I., Finite Element Analysis, J. Wiley, 1991. [8] Weiler, K.J., “The radial-edge structure: A topological representation for non-manifold geometric boundary representations”, M.J. Wozny, H.W. McLaughlin, J.L. Encarnacao, editors. Geometric modeling for CAD applications, North Holland, pp. 3-36, 1988. [9] Muskhelishvili, N.I., Some Basic Problems of the Mathematical Theory of Elasticity, 1963 Groningen, P. Noordhoff.
