INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING Int. J. Numer. Meth. Engng 2000; 00:1–6 Prepared using nmeauth.cls [Version: 2002/09/18 v2.02]

The hp-d adaptive Finite Cell Method for geometrically nonlinear problems of solid mechanics D. Schillinger1,∗,† , A. D¨ uster2 and E. Rank1 1 Lehrstuhl

f¨ ur Computation in Engineering, Technische Universit¨ at M¨ unchen, Arcisstr. 21, 80333 M¨ unchen, Germany 2 Numerische Strukturanalyse mit Anwendungen in der Schiffstechnik, Technische Universit¨ at Hamburg-Harburg, Schwarzenbergstrasse 95c, 21073 Hamburg, Germany

SUMMARY The Finite Cell Method (FCM) combines the fictitious domain approach with the p-version of the finite element method and adaptive integration. For problems of linear elasticity, it offers high convergence rates and simple mesh generation irrespective of the geometric complexity involved. The present article presents the integration of the Finite Cell Method into the framework of nonlinear finite element technology. However, the penalty parameter of the fictitious domain is restricted to few orders of magnitude in order to maintain local uniqueness of the deformation map. As a consequence of the weak penalization, nonlinear strain measures provoke excessive stress oscillations in the cells cut by geometric boundaries, leading to a low algebraic rate of convergence. Therefore, the FCM approach is complemented by a local overlay of linear hierarchical basis functions in the sense of the hp-d method, which synergetically use the h-adaptivity of the integration scheme. Numerical experiments show that the hp-d overlay effectively reduces oscillations and permits stronger penalization of the fictitious domain by stabilizing the deformation map. The hp-d adaptive Finite Cell Method is thus able to restore high convergence rates for the geometrically nonlinear case, while preserving the easy meshing property of the original FCM. Accuracy and performance of the present scheme are demonstrated by several benchmark problems in 1, 2 and 3 dimensions, and the nonlinear simulation of a complex foam c 2000 John Wiley & Sons, Ltd. sample. Copyright key words: Finite Cell Method; Fictitious domain methods; High-order finite elements; Hierarchical hp-d adaptivity; Geometric nonlinearity; Complex geometries

1. INTRODUCTION Structural analysis with standard finite elements requires the discretization of the domain of interest into a finite element mesh, whose boundaries conform to the physical boundaries of the structure [1, 2, 3]. While this constraint can be easily achieved for many applications in solid mechanics, it constitutes a severe bottleneck for structures of highly complex geometry.

∗ Correspondence to: Dominik Schillinger, Lehrstuhl f¨ ur Computation in Engineering, Department of Civil Engineering and Surveying, Technische Universit¨ at M¨ unchen, Arcisstr. 21, 80333 M¨ unchen, Germany † E-mail: [email protected]

c 2000 John Wiley & Sons, Ltd. Copyright

2

¨ D. SCHILLINGER, A. DUSTER AND E. RANK

The translation of corresponding geometrical models, usually available in the form of CAD or spatial voxel data, into conforming finite element discretizations is computationally very expensive, hard to be fully automated and often leads to error-prone meshes, which have to be improved manually by the user [4]. The meshing challenge has recently led to the rise of isogeometric methods, which directly use the spline basis of a CAD model for FE analysis (see for example [4, 5, 6] and the references therein). A more general approach to avoid time-consuming mesh generation for complex domains is provided by the fictitious domain concept [7, 8, 9]. Its main idea consists of the extension of the physical domain of interest beyond its potentially complex boundaries into a larger embedding domain of simple geometry, which can be meshed easily by a structured grid (see Figure 1). To preserve consistency with the original problem, the influence of the fictitious domain extension is extinguished by penalizing its material parameters. The geometry information is thus incorporated implicitly during evaluation of the stiffness matrix via the difference in material parameters at integration point level. The fictitious domain approach has been widely used in conjunction with penalty methods [9, 10, 11], the mortar approach [12], Lagrange multipliers [8, 13, 14, 15], Nitsche’s method [16], extended finite elements [17, 18, 19, 20], discontinuous Galerkin [21] and spectral methods [22, 23] to address problems of structural analysis, acoustics, heat conduction, fluid flow, fluid-structure interaction and optimization. It has also been combined with sensitivity analysis and level set methods for problems of topology and shape optimization [24]. The Finite Cell Method (FCM) [25, 26] combines the fictitious domain approach with the p-version of the finite element method [27, 28] and adaptive integration. For smooth problems of linear elasticity, FCM has been shown to maintain exponential rates of convergence in the energy norm known from the p-version, based on a smooth extension of the solution fields into the fictitious domain [25]. FCM thus allows for accurate structural analysis irrespective of the geometric complexity involved [29], and can be well combined with voxel based and implicit geometrical models typical for applications from biomechanics and material science [26, 30]. The advantages of FCM in terms of fast convergence at almost no meshing cost have also been successfully applied for fluid-structure interaction [31], topology optimization [32] and thin-walled structures [33]. The present contribution deals with the extension of the Finite Cell Method to geometrically nonlinear problems of solid mechanics by integrating the FCM approach into the established framework of nonlinear finite element technology [34, 35, 36]. In the resulting nonlinear FCM formulation, however, the penalization of the fictitious domain turns out to considerably affect the stability of the deformation gradient. In order to maintain the constraint of a one-to-one mapping without material penetration, the stiffness in the fictitious domain can be reduced only by few orders of magnitude in comparison to the physical stiffness, which weakens the elimination of the strain energy contribution of the fictitious domain. Consequently, solution fields develop excessive oscillations in stresses in the cells cut by geometric boundaries, which slows down convergence to a low algebraic rate. Therefore, an unfitted h-refinement around geometric boundaries is introduced, based on a linear hierarchical basis, which is added to the high-order basis in the sense of an hp-d overlay [37, 38, 39, 40]. The overlay can synergetically operate on the hierarchical sub-cell structure, built up by the adaptive integration scheme of the original FCM formulation [26], so that the full simplicity in terms of mesh generation is preserved. Numerical experiments show that the addition of hierarchical basis functions stabilizes the deformation gradient, thus permitting penalty parameters of several orders of c 2000 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2000; 00:1–6

3

hp-d FCM FOR GEOMETRICALLY NONLINEAR PROBLEMS

∂Ω

Ωfict

ΓN

α ≪ 1.0

Ω=Ωphys +Ωfict

Ωphys

α = 1.0

ΓD Figure 1: The physical domain Ωphys is extended by the fictitious domain Ωfict into the embedding domain Ω to allow easy meshing of complex geometries. The influence of Ωfict is penalized by α.

magnitude smaller than in the unrefined standard FCM. Moreover, the h-adaptivity of the overlay effectively localizes oscillations to the close vicinity of geometric boundaries. If the polynomial degree p of the high-order basis and the depth k of the adaptive hierarchical overlay are increased simultaneously, high rates of convergence can be achieved also for geometrically nonlinear problems. In summary, the main novelty of the article is the introduction of a hierarchical refinement technique in the sense of the hp-d approach, which in combination with the Finite Cell Method maintains simple mesh generation and preserves high rates of convergence for geometrically nonlinear fictitious domain problems by stabilization of the deformation gradient. The core idea of the hp-d FCM bears similarities to the recently introduced hierarchical B-spline FEM [41, 42]. The present paper is organized as follows: Section 2 gives a short review of the Finite Cell Method for linear problems. Section 3 briefly summarizes the formulation of geometric nonlinearity in principal directions and explains its influence on the standard FCM scheme. Section 4 introduces the concept of FCM with a hp-d overlay around geometric boundaries. Section 5 generalizes the hp-d FCM to 2 and 3 dimensions, outlines the imposition of unfitted boundary conditions, and introduces implicit geometrical models as a basis for fast hp-d grid generation. Finally, Section 6 presents a series of numerical benchmarks and the application example of a nonlinear compression test for a foam sample, while Section 7 briefly outlines the potential of the hp-d FCM to handle material interfaces and singularities.

2. THE FINITE CELL METHOD FOR LINEAR PROBLEMS In the following, the principles of the standard Finite Cell Method for problems of linear elasticity and its characteristic solution behavior are briefly reviewed. Details and further numerical examples can be found in [25, 26, 33]. 2.1. Variational formulation According to Figure 1, the embedding domain Ω consists of the physical domain of interest Ωphys and the fictitious domain extension Ωfict . Analogous to standard finite elements, the Finite Cell Method for linear elastic problems is derived from the principle of virtual work δW (u, δu) =

Z

σ : (∇sym δu) dV − Ω

c 2000 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls

Z

δu · b dV − Ωphys

Z

δu · t dA = 0

(1)

ΓN

Int. J. Numer. Meth. Engng 2000; 00:1–6

¨ D. SCHILLINGER, A. DUSTER AND E. RANK

4

where σ, b, u, δu and ∇sym denote the stress tensor, body forces, displacement vector, test functions and the symmetric part of the gradient, respectively [1, 2, 3]. Force boundary conditions are specified over the boundary of the embedding domain ∂Ω, where tractions are zero by definition, and over the Neumann boundary ΓN of the physical domain by traction vector t (see Figure 1). The standard elasticity tensor C [1, 2, 3] relating stresses and strains σ = αC : ε

(2)

is complemented by a scalar factor α, which leaves the material parameters unchanged in the physical domain, but penalizes the contribution of the fictitious domain as follows ( 1.0 ∀x ∈ Ωphys α (x) = (3) 10−q ∀x ∈ Ωfict In Ωfict , α must be chosen as small as possible, but large enough to prevent extreme illconditioning of the stiffness matrix [25, 26]. Typical values of α range between 10−4 and 10−8 . If α is very small or even zero, shape functions with little support in Ωphys lead to stiffness matrix entries, which are considerably smaller than entries resulting from shape functions located completely inside Ωphys . In conjunction with p-refinement, this leads to a considerable increase of the corresponding condition number, so that a solution with iterative solvers requires a prohibitive number of iterations. Therefore, direct solvers will be applied in the following. Using a structured grid of high-order finite elements (see Figure 1), which will be called finite cells in the following, kinematical quantities can be discretized as u=

n X

Na u a

(4)

Na δua

(5)

a=1

δu =

n X

a=1

The sum of Na denotes a finite set of n shape functions of the p-version of the finite element method, derived from integrated Legendre polynomials [27, 28], and ua and δua the corresponding vector of unknown coefficients [2, 3, 35]. The geometry of the structured finite cell mesh in terms of the position vector X can be represented exactly by the nodal part of the high-oder basis alone n vert X Nilin Xi (6) X= i=1

Nilin

where and Xi are linear basis functions and positions of the nvert cell vertices, respectively. The Finite Cell Method can thus be regarded as a sub-parametric finite element scheme [1]. Following the standard Bubnov-Galerkin approach [2, 3], inserting Equations (4) and (5) into the weak form Equation (1) produces a discrete finite cell representation Kua = f

(7)

with stiffness matrix K and load vector f . Displacement boundary conditions defined over the Dirichlet boundary ΓD (see Figure 1) can be integrated into the discrete system Equation (7) by penalty [9, 25, 26] or variational methods [8, 43, 44]. Due to the similarity to standard FEM, the implementation of FCM can exploit existing finite element techniques to the full. c 2000 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2000; 00:1–6

5

hp-d FCM FOR GEOMETRICALLY NONLINEAR PROBLEMS

k=0

k=1

k=2

k=3

k=4

k=5

Finite cell mesh with geometric boundary

Figure 2: 2D sub-cell structure for adaptive integration of cells cut by the geometric boundary. The quadtree is built by iterating over the sub-cells of current level k, until the max. depth is reached.

2.2. Adaptive integration of discontinuous cells The accuracy of numerical integration by Gauss quadrature [2, 3], which assumes smoothness of the integrands, is considerably influenced by discontinuities within cells introduced by the penalization parameter α of Equation (3). Therefore, the Finite Cell Method uses composed Gauss quadrature to improve integration accuracy in cells cut by the geometric boundaries, based on a hierarchical decomposition of the original cell into integration sub-cells [26]. In 2 dimensions, the sub-cell structure can be built up in the sense of a quadtree [45] (see Figure 2). Starting from the original finite cell of level k=0, each sub-cell of level k=i is first checked whether it is cut by a geometric boundary. If true, it is replaced by 4 equally spaced cells of level k=i+1, each of which is equipped with (p+1)×(p+1) Gauss points. Partitioning is repeated for all cells of current level k, until a predefined maximum depth k=m is reached, leading to an adaptive aggregation of integration sub-cells around geometric boundaries. The main idea of the adaptive integration scheme is to preserve the optimal accuracy of Gauss quadrature for all larger sub-cells of depth k
Int. J. Numer. Meth. Engng 2000; 00:1–6

replacemen ¨ D. SCHILLINGER, A. DUSTER AND E. RANK

6

fsin

∆u Parameters:

X

Ωphys L

Ωfict

Ωphys

4/3 L

2/3 L

Young’s modulus E=1.0 Poisson’s ratio ν = 0.0 Penalization parameter α = 10−q Area A=1.0 Length of each part L=1.0 Displacement load ∆u = 1.0 Sine load fsin = 1/20 sin (4πX)

Discretization with 2 finite cells:

Figure 3: Uni-axial rod example. Geometric boundaries are located at X = 1.0 and X = 2 13 .

unaffected. As shown in Figure 4, the corresponding analytical solution with finite α=10−8 , which can be derived from basic elastostatics, exhibits kinks and jumps in displacements and stresses, respectively, at geometric boundaries X=L and X=7/3. The FCM solution, however, extends smoothly into the fictitious domain despite the discontinuities of the analytical solution, illustrated for a discretization with 2 equally-spaced cells of p=12 in Figure 4. The importance of the penalization parameter and the smooth extension of the FCM solution into the fictitious domain can be explained with the help of the total strain energy U=

Z

Ψ dV = Ω

Z

σ : ε dV

(8)

Ω

where Ψ represents the strain energy function, defined over the complete domain Ω. In linear elasticity and nonlinear elasticity based on the logarithmic strain measure, Ψ can be obtained by taking the scalar product of Cauchy stresses σ and strains ε. The best approximation property to the total strain energy states that the solution of a Galerkin finite element scheme represents a least-squares best fit to the exact solution in terms of Equation (8) [2, 47]. Due to the penalization with parameter α, deviations from the exact solution in Ωfict have a considerably smaller impact on the strain energy Equation (8) than deviations in Ωphys . Therefore, a minimization of the strain energy error by the high-order basis of the FCM scheme results in an accurate approximation in Ωphys , where potential deviations lead to considerable error contributions. In Ωfict , a largely inaccurate approximation is allowed, since potential deviations lead to negligible error contributions due to penalization. In particular, this implies a smooth extension of the solution into the fictitious domain, so that its gradients in Ωphys remain unaffected up to the geometric boundary.

3. GEOMETRIC NONLINEARITY AND THE FINITE CELL METHOD The Finite Cell Method is now adapted to the geometrically nonlinear formulation based on logarithmic strains, and the resulting solution characteristics are discussed for the 1D example. Details of the nonlinear formulation can be found for instance in [35, 36, 48, 49, 50]. c 2000 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2000; 00:1–6

7

hp-d FCM FOR GEOMETRICALLY NONLINEAR PROBLEMS

4 0.0 2

-0.5

Strain ε

Axial strains ε

3

-1.0

1 0 -1

Analytical Linear FCM -2.0

εlin = λ − 1 εlog = ln λ
εGL = 21 λ2 − 1
1 εA = 2 1 − 1/λ2

-2

-1.5

0

-3

1 2 Initial rod length X

-4

3

Figure 4: Smooth extension of the FCM solution vs. discontinuous analytical solution for the linear elastic strain field (α=10−8 ).

0.0

0.5

1.0

1.5 2.0 Stretch λ

2.5

3.0

Figure 5: Engineering, logarithmic, GreenLagrange and Almansi strain measures εlin , εlog , εGL , εA , respectively, in 1 dimension.

3.1. Kinematics The deformation map ϕ describes the motion of each material particle from its initial reference configuration X to its deformed spatial configuration x. x = ϕ(X)

(9)

The displacement field u and the deformation gradient F follow as u=x−X

(10)

∂x (11) ∂X In view of its polar decomposition, F=VR can be represented by the rotation tensor R and the spatial stretch tensor V. Thus, V2 can be obtained from the left Cauchy-Green tensor F=

b = FFT = (VR)(RT V) = V2

(12)

Evaluation of the principal directions of V2 yields the orthogonal spatial Eigenvector triad {n1 ,n2 ,n3 } and the associated Eigenvalues {λ21 , λ22 , λ23 }, which are identified as the squared principal stretches. The spatial logarithmic strain tensor ε = ln V in spectral form follows as ε=

3 X

ln λa na ⊗ na

(13)

a=1

Logarithmic strains, also known as true or natural strains, represent a suitable strain measure for the entire deformation range [35, 36], and are illustrated for the 1D case in Figure 5. Due to the discontinuous penalization parameter α of Equation (3), the deformation map and corresponding displacements Equations (9) and (10) exhibit a weak discontinuity (kink) along the geometric boundaries of the physical domain Ωphys , whereas the deformation gradient and related stretches and strains Equations (11) to (13) exhibit a strong discontinuity (jump). c 2000 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2000; 00:1–6

¨ D. SCHILLINGER, A. DUSTER AND E. RANK

8

3.2. Constitutive equations in principal directions The Hencky hyperelastic model is the finite logarithmic strain based extension of the standard linear elastic material, whose strain energy function Ψ reads     E νE 2 2 2 2 Ψ=α· (14) (ln J) (ln λ1 ) + (ln λ2 ) + (ln λ3 ) + 2(1 + ν) 2(1 + ν)(1 − 2ν)

with J = det F = λ1 λ2 λ3 and material parameters Young’s modulus E and Poisson’s ratio ν. In a Finite Cell formulation, Young’s modulus E is penalized by α of Equation (3). The principal Cauchy stresses along principal axes a={1, 2, 3} follow from Equation (14) as   E α νE 1 ∂Ψ = ln λa + ln J (15) σa = J ∂ ln λa J (1 + ν) (1 + ν)(1 − 2ν) The fourth order spatial elasticity tensor in Cartesian coordinates can then be computed as cijkl =

3 3 3 X X X ∂2Ψ σa λ2b − σb λ2a 1 2σa η aaaa + (η abab + η abba ) (16) η aabb − J ∂ ln λa ∂ ln λb λ2a − λ2b a=1 α,β=1

a,b=1

implying the fourth order dyadic product η ijkl = ni ⊗ nj ⊗ nk ⊗ nl . Due to the strong discontinuity in α, all derived quantities exhibit jumps along the geometric boundaries. 3.3. Discretization, linearization and the Newton-Raphson procedure Taking into account Equations (9) to (16), the variational formulation of the FCM for geometrically nonlinear elastic problems can be derived from the principle of virtual work Z Z Z δW (ϕ, δv) = σ : ∇x δv dv − b · δv dv − t · δv da = 0 (17) ϕ(Ωphys )

ϕ(Ω)

ϕ(ΓN )

with body forces b, traction vector t and test function δv. The basic kinematical quantities u, δv and F can be discretized in the sense of Equations (4), (5) and by F=I+

n X

∇X Na ua

(18)

a=1

Inserting the displacement approximation Equation (4) and the reference configuration, interpolated subparametrically by Equation (6), into Equation (10) yields the interpolation of the deformed spatial configuation x=

n vert X

Nilin Xi +

n X

Na u a

(19)

a=1

i=1

Using these expressions in Equation (17), the discretized virtual work per coefficient a can be formulated as the difference between the internal and external equivalent force vectors f int and f ext , called residual r
(20) δW (ϕ, Na δva ) = δva · faint − faext = δva · ra Z faint = BTa σ dv (21) ϕ(Ω)

c 2000 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2000; 00:1–6

9

hp-d FCM FOR GEOMETRICALLY NONLINEAR PROBLEMS

(a) Exact solution fields.

(b) Solution fields obtained with 2 finite cells of p=12.

(c) Solution fields obtained with 2 finite cells of p=12 and a hierarchical hp-d overlay of depth 11. Figure 6: Stretch, strain and stress solutions of the geometrically nonlinear rod example of Figure 3, corresponding to a penalty parameter α=10−2 .

faext

=

Z

Na b dv + ϕ(Ωphys )

Z

Na t da

(22)

ϕ(ΓN )

where B denotes the B-matrix. The consistent linearization of the discretized weak form DδW of Equation (17) in the direction of incremental coefficients ∆ua is described in detail in [35, 36, 48, 49] and leads to a discrete system of the form (Kc + Kσ + Kt ) ∆ua = −r

(23)

where Kc , Kσ and Kt denote material, geometric and deformation dependent tangent stiffness matrices, respectively. The classical Newton-Raphson procedure iteratively solves the linearized system Equation (23) for ∆ua to update the total coefficients ua , until the norm of the residual vector r has converged below a tolerance close to zero. c 2000 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2000; 00:1–6

10

¨ D. SCHILLINGER, A. DUSTER AND E. RANK

3.4. FCM solution characteristics for a geometrically nonlinear 1D example The geometrically nonlinear formulation for the rod of Figure 3 simplifies in the 1D case to Ψ=α

E 2 (ln λ) 2

(24)

E ln λ J

(25)

E − 2σ J

(26)

σ=α

c=α

with axial stretch λ and the determinant of the deformation gradient J = λ1−2ν [35]. Taking into account the displacement load ∆u only (see Figure 3), the exact solution could even be found by a conforming discretization with 3 standard linear elements. Stretch, strain and stress are shown in Figure 6a for 10 displacement increments. A geometrically nonlinear formulation using the standard FCM of Section 2 is implemented within MATLAB [51]. Numerical experiments reveal that the lower bound of the penalty parameter for the current 1D example is α=10−2 . Corresponding solution fields obtained with the two cell discretization of Figure 3 are plotted in Figure 6b. For α<10−2 , the deformation map loses its uniqueness within the fictitious domain Ωfict . Consequently, the determinant of the deformation gradient F of Equation (11) falls below zero at some integration points, which inevitably terminates the computation. From a physical point of view, this can be interpreted as a local penetration of material, violating the principles of continuum mechanics [35, 49]. With α as large as 10−2 , the penalization of Equation (3) does not sufficiently eliminate the influence of Ωfict , so that the nature of the FCM formulation changes from a fictitious domain to an interface problem. In addition, the strain energy contribution of Ωfict is amplified by the logarithmic strain measure, which is illustrated in Figure 5 with respect to the stretch λ = l/L, where L and l denote initial and deformed lengths, respectively. Whereas engineering strains of linear elasticity are bound by definition to very small values |εlin | ≪ 1.0, logarithmic strains of nonlinear elasticity are able to grow without bounds in order to yield physically meaningful measures for very large deformation states. However, in the case of large deformations in Ωfict , nonlinear strains act as an additional counterbalance to α. As a consequence, the contribution of Ωfict to the total strain energy Equation (8) grows, so that the standard FCM scheme tries to accurately fit the solution in both the fictitious and physical domains due to the best approximation property discussed in Section 2.4. Thus, solution fields do not extend smoothly into Ωfict as shown in Figure 4, but develop large oscillations throughout the discontinuous cells (see Figure 6b). The corresponding convergence behavior deteriorates to a low algebraic rate, which is a well-known issue for high-order finite elements with inter-element discontinuities [27, 52]. Besides logarithmic strains, there exists a range of nonlinear strain measures, e.g. Green-Lagrange or Almansi strains (see [35, 36] and Figure 5 for an illustration in 1D). Numerical experiments indicate that geometrically nonlinear FCM formulations based on other strain measures in combination with corresponding constitutive equations affect the uniqueness of the deformation map in the same way. c 2000 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2000; 00:1–6

hp-d FCM FOR GEOMETRICALLY NONLINEAR PROBLEMS

Base mesh: up

Overlay mesh: uh

11

Superposition: uhp-d

Figure 7: Basic idea of the hp-d method in 1 dimension.

4. THE CONCEPT OF THE hp-d FINITE CELL METHOD IN 1 DIMENSION The present section introduces the fundamental principles of the hp-d FCM in 1D, based on a synergetic combination of the adaptivity of integration sub-cells with a hp-d overlay of hierarchical basis functions, and illustrates its improvements for the nonlinear example. 4.1. The hp-d method Figure 7 outlines the principle of the hp-d method. The starting point is a standard high-order finite element approximation up on a so-called base mesh, discretizing the whole computational domain. In order to locally improve this solution - for example, in a region where a significant error occurs - a h-version approximation uh defined on a fine overlay mesh is superimposed to the base mesh. The overall finite element approximation uhp-d = up + uh results from the superposition of base and overlay approximations. C 0 -continuity of uhp-d is guaranteed by imposing homogeneous Dirichlet boundary conditions on the overlay approximation uh . Originally proposed for adaptive domain decomposition (hence the d) [37], the hp-d method has been successfully adapted for a range of applications, such as multi-scale structural analysis [38, 39, 53] and model adaptivity [40, 54]. 4.2. Integrating hp-d into the Finite Cell Method The basic idea of the hp-d FCM is to complement the high-order finite cell discretization with an adaptive overlay of h-version basis functions around geometric boundaries in order to locally improve the approximation of discontinuities and effectively reduce oscillations. The adaptive integration scheme described in Section 2.2 already realizes h-adaptivity around geometric boundaries in the form of sub-cells, however at this point only for integration purposes. Analogous to the two-dimensional case of Figure 2, the sub-cell structure in 1D is organized in the form of a binary tree, which partitions sub-cells of level k=i cut by the geometric boundary into 2 equally spaced sub-cells of level k=i+1. The resulting sub-cell structure shows strong parallels to a hierarchical representation of a linear basis [55, 56, 57] in the sense that both are organized in levels k and carry out a bisection in each level (see Figure 8a). The key concept for a computationally efficient implementation of the hp-d FCM is to synergetically use the sub-cell structure to construct an adaptive overlay from linear hierarchical basis functions. In each level k > 0 of the binary tree, those hat functions from level k of the hierarchical basis are added, which have their nodal degrees of freedom between two newly generated sub-cells of the adaptive integration structure. This procedure is illustrated in Figure 8b for the first cell from the FCM discretization of Figure 3. In each level k, the sub-cell containing the geometric boundary is bisected and the corresponding hat function is added to the overlay. The basis c 2000 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2000; 00:1–6

12

¨ D. SCHILLINGER, A. DUSTER AND E. RANK

k=0 (High-order basis)

1 0 -1 k=1

k=1 1 0

k=2

k=2 1 0

k=3

k=3 1 0

k=4

k=4 1 0 0.0

0.5

1.0

1.5

Initial rod length X

(a) Hierarchical representation of the linear basis according to levels k.

(b) Adaptive hp-d overlay for the first cell of Figure 3. The discontinuity is located at X=1.0.

Figure 8: Principle of the hp-d Finite Cell Method.

functions of the overlay achieve a rapid concentration of nodes around the discontinuity. The hp-d FCM in this form can be easily implemented into existing FCM codes by the addition of the hierarchical overlay. Since it operates directly on the existing sub-cell structure of the adaptive integration scheme, it requires no additional mesh generation and completely retains the easy meshing property of the original FCM. Its hierarchical adaptivity considerably reduces the number of additional degrees of freedom in comparison to global h-refinement of cells or an overlay with equally spaced h-version shape functions. In Figure 6c, it is illustrated for the current example that the hp-d FCM leads to an effective reduction of oscillations, which are restricted to the close vicinity of geometric boundaries. 4.3. Re-enforced penalization by hp-d stabilization The general validity of the FCM approach largely depends on the effective penalization in the sense of Equation (3), for which parameter α needs to be chosen as small as possible. Whereas in the linear elastic case, the lower bound of α is determined by the level of illconditioning of the stiffness matrix allowed by the equation solver at hand, it is dictated in the geometrically nonlinear case by the local uniqueness of the deformation map in Ωfict . As shown in Section 3.4, this requirement limits the difference in stiffness between Ωphys and Ωfict to few orders of magnitude, which weakens the penalization of Ωfict . Numerical c 2000 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2000; 00:1–6

hp-d FCM FOR GEOMETRICALLY NONLINEAR PROBLEMS

13

experiments for the 1D example problem show that the hp-d overlay leads to a considerable stabilization in the numerical solution of the deformation map, which allows for a decrease of the penalty parameter α by several orders of magnitude in comparison to the unrefined standard FCM. While α is bound to 10−2 in an analysis with the high-order base mesh only (see Section 3.4), the superposition of overlays of m=2 and m=5 levels permit a locally unique analysis with penalty parameters α=10−4 and 10−6 , respectively. Hence, in addition to better approximating discontinuities and reducing oscillations, the hp-d adaptive FCM allows for a stronger penalization of Ωfict , thus reassuring the validity of the fictitious domain concept. 4.4. Convergence behavior of the hp-d FCM in 1D Accuracy and convergence of the FCM versions will be assessed by s |Uex − UFCM | er = × 100% Uex

(27)

denoting the relative error in terms of the total strain energy U of Equation (8) [2, 47, 58]. For fictitious domain problems, Uex denotes the exact strain energy of the original problem defined only over the physical domain Ωphys , corresponding to the first picture of Figure 1. Accordingly, UFCM denotes the strain energy contribution from Ωphys , obtained numerically with an FCM scheme. The corresponding rate of convergence q, with which er decreases from the ith to the (i+1) th refinement step compared to the number of degrees of freedom ndof is
i log10 ei+1 r /er   (28) q=− i log10 ni+1 dof /ndof

For continuous problems in 1D, optimal rates of convergence are equivalent to the polynomial degree of the basis functions involved, if h-refinement is applied [47], whereas with p-refinement, exponential rates of convergence can be achieved [27]. In the following, the uni-axial rod of Figure 3 is considered with sine load fsin . The physical reference strain energy Uex consists of the energy contribution caused by fsin only, since the displacement load at the right boundary triggers a rigid body movement of the right rod. An overkill discretization with 1000 conforming cubic finite elements taking into account the left rod loaded by fsin yields Uex,l =2.374721186 · 10−5 and Uex,nl =1.171825900 · 10−5 for the linear and geometrically nonlinear case, respectively. Figure 9 shows the convergence behavior of the presented FCM schemes, if the polynomial degree p of the two cell discretization shown in Figure 3 is increased from 1 to 20. Note that for all computations, an adaptive sub-cell structure of depth k=20 and p+1 Gauss points is used to minimize the integration error. In the linear elastic case, FCM converges exponentially with a maximum rate of q=10.14. The penalization parameter α=10−8 cannot completely erase the influence of the fictitious domain, and the corresponding error takes control at a value below 0.6%, which leads to a flattening of the convergence curve. In the geometrically nonlinear case, the one-to-one requirement of the deformation map dictates α=10−2 as discussed in Section 3.4. The weakened penalization is not able to sufficiently erase the influence of the fictitious domain, so that the error level with respect to the true solution is very high, which can be also interpreted as a modeling error. In addition, the large α in conjunction with the nonlinear strain measure slows down the corresponding convergence to an algebraic rate of q=0.78 on average. If in addition to the c 2000 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2000; 00:1–6

¨ D. SCHILLINGER, A. DUSTER AND E. RANK

Relative error in strain energy measure [%]

14

10

4

10

3

10

2

10

1

10

0

10

-1

Linear FCM, α=10-8 Nonlinear FCM, α=10-2 Nonlinear hp-d FCM, α=10-6 3

5

10 20 Degrees of freedom

50

100

Figure 9: Strain energy convergence for the uni-axial rod example.

elevation of the polynomial degree p, an overlay of k=p levels of adaptive hierarchical functions are added to the high-order finite element basis in the sense of Figure 8b, the resulting hp-d FCM scheme remains stable for a penalty parameter of α=10−6 . Due to the strengthened penalization, the error level immediately decreases by more than 2 orders of magnitude. Due to the effective reduction of oscillations around geometric boundaries, the hp-d FCM achieves exponential convergence with a maximum rate of q=7.69 for the geometrically nonlinear case. The flattening of the convergence curve due to the finite value of α occurs at a relative error of below 5%, so that for the present example, the hp-d FCM is almost 3 orders of magnitude more accurate than the unrefined standard FCM. 4.5. A note on very large deformation states and incompressibility Robustness under severe mesh distortion and in the presence of an incompressibility constraint constitute important characteristics of a geometrically nonlinear formulation. These aspects are not examined in the scope of the present paper, but some comments shall be given here. Mesh distortion as a result of very large deformation states can lead to negative values in the Jacobian of the geometrical mapping [50], which is a violation of the principles of continuum mechanics (see also Section 3.4). High-order schemes related to the p-version of the FEM have been shown to permit increased levels of mesh distortion [59, 60]. Due to its structured undistorted mesh in the reference configuration, distortion of the FCM mesh results solely from the deformation. The geometrically nonlinear Finite Cell Method can be shown to inherit the robustness of its high-order basis despite the presence of a fictitious domain, and is thus able to represent very large deformation states without further changes in the formulation (see also [61, 62]). Incompressible materials such as rubber lead to an ill-conditioned problem, whose numerical solution requires improved finite element formulations due to volumetric locking [63, 64]. In contrast to standard h-version finite elements, the p-version has been shown to offer increased robustness against locking [65, 66] and to even be locking-free, when applied to nearly incompressible Neo-Hookean materials under finite deformations [67]. Similar to the case of severe mesh distortion, the Finite Cell Method can be expected to inherit the benefits of its high-order basis in terms of locking-free behavior under incompressibility constraints. c 2000 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2000; 00:1–6

15

hp-d FCM FOR GEOMETRICALLY NONLINEAR PROBLEMS

1

0.5

0.75 0

0.5

0.25

0.25 0 -1

1

-0.25

1

0.5 0

-0.5

ξ

0

0.5

-0.5 1

-1

η

1

-0.5 0.5 -1

0

-0.5

ξ

Nodal mode

0

0.5

-0.5 1

Edge mode

η

0 -1

0.5 0

-0.5

ξ

0

0.5

-0.5 1

η

-1

-1

Internal mode

Figure 10: Examples of nodal, edge and internal modes of the two-dimensional high-order basis.

5. THE hp-d FINITE CELL METHOD IN 2 AND 3 DIMENSIONS In the following, the straightforward generalization of the hierarchical hp-d overlay for higher dimensions, the incorporation of unfitted boundary conditions and suitable implicit geometrical models for fast hp-d grid generation are addressed. 5.1. Generalizing the concept of hierarchical overlays in 2 and 3 dimensions FCM for 2 and 3 dimensional problems of solid mechanics [25, 26] uses the so-called trunk space of the p-version of the finite elements [27, 28] over a structured finite cell mesh. For the 2D case, the FCM concept, the adaptive integration scheme and examples of high-order shape functions are illustrated in Figures 1, 2 and 10, respectively. The hp-d idea introduced in the previous section creates an overlay of basis functions, selected from the hierarchical representation of a linear basis in such a way that they reproduce the adaptivity of integration sub-cells to attain a local aggregation of nodal values around geometric boundaries. The refinement procedure is achieved as follows (see also Figure 11). We start from a number of finite cells cut by a geometric boundary, which are equipped with high-order basis functions up to polynomial degree p and several levels k of adaptive integration sub-cells. Looking at the hierarchical representation of a two-dimensional linear basis in Figure 11, we note that its arrangement into levels k corresponds to the organization of quadtree based integration sub-cells (see also Figure 2). Next, we successively iterate over each level k of hierarchical sub-cells and add a hat function of the corresponding hierarchical level k to the finite cell basis, if all 4 sub-cells around a nodal position are contained in the adaptive integration structure. Hence, there are 5 potential candidates for each 2D sub-cell at the locations shown in Figure 12a, where the outer nodal positions on the 4 edges are only selected, if the corresponding neighboring subcells are also existing. In the 3D hexahedral cell of Figure 12b, each sub-cell comprises 19 potential nodes, where again the outer positions are only selectable, if all sub-cells touching the corresponding edges or faces are partitioned up to the same level k. We note that the finest level k=m of basis functions of the overlay requires a depth of the hierarchical sub-cell structure equal or larger than m. The hp-d overlay in this form fully re-uses the adaptivity provided by the sub-cell based integration structure and does not require the introduction of additional Gauss points or sub-cells. It is therefore easy to implement, maintains fast mesh generation irrespective of the geometric complexity, and achieves a high degree of local adaptivity around geometric boundaries at a limited amount of additional degrees of freedom. c 2000 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2000; 00:1–6

¨ D. SCHILLINGER, A. DUSTER AND E. RANK

16

Hierarchical representation of linear basis Nodal modes of the high-order basis Nodes of 1st hierarchical level Nodes of 2nd hierarchical level Geometric boundary k=1

k=0

Finite cell mesh

k=2

hp-d overlay on adaptive sub-cells

Figure 11: Principle of the hp-d Finite Cell Method.

5.2. Unfitted boundary conditions In case the physical domain Ωphys conforms to the boundary of the structured finite cell mesh, boundary conditions can be implemented by standard FE techniques [2, 3, 27, 28]. Geometric boundaries cutting through cells require special approaches described in the following. 5.2.1. Dirichlet boundary conditions The incorporation of homogeneous Dirichlet boundary conditions can be easily achieved by penalization. Instead of choosing a very small value, parameter α of Equation (3) is increased in Ωfict by some orders of magnitude, so that Dirichlet conditions applied to the boundaries of the finite cell mesh are directly carried on to the boundary of Ωphys [9, 25, 26]. A more accurate way to incorporate general Dirichlet boundary conditions is provided by variational methods, which lead to an enforcement in a weak sense. In the scope of the present article, the interested reader is referred to [8, 30, 44, 68, 69, 70]. 5.2.2. Neumann boundary conditions The fictitious domain concept inherently satisfies Neumann boundary conditions of zero traction, since stresses cannot be transferred beyond Ωphys due to the stiffness penalization in Ωfict [25, 26]. Non-zero Neumann boundary conditions require an integration over the boundary ΓN (see Figure 1), which therefore needs to be parameterized as a curve or surface [26]. Furthermore, in the framework of geometric nonlinearity, the deformation dependence of the loading needs to be taken into account [35, 49]. In the scope of the present paper, the implementation of deformation dependent loads in the framework of the hp-d FCM is discussed for a normal pressure pˆ along an arbitrary twodimensional polygon [73]. This procedure can be easily adapted to 3D problems by introducing a triangular approximation of boundary surfaces. From Equation (17), the virtual work due c 2000 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2000; 00:1–6

17

hp-d FCM FOR GEOMETRICALLY NONLINEAR PROBLEMS

ζ

ξ η

ξ

η

(b) 3D octree based sub-cells.

(a) 2D quadtree based sub-cells.

Figure 12: Sub-cell partitioning in 2 and 3 dimensions and corresponding potential nodal candidates for the hierarchical hp-d overlay (open circles).

to pˆ acting on ΓN follows with unit normal vector n and traction vector ˆt = pˆn as Z δW (ϕ, δv) = pˆn · δv da

(29)

ϕ(ΓN )

The algorithmic treatment of Equation (29) is based on a re-parameterization of n in terms ˆ e , which can be calculated of convective coordinates. In 2D, this leads to a normal vector n pointwise from the gradient of the deformed boundary curve with respect to the coordinate ϑ of its parameter space (see [71, 72, 73] for an in-depth discussion). In FCM, geometric boundaries are defined independently of the structured cell discretization (see Figure 1). ΓN can thus be parameterized by an arbitrary curve, approximated by the polygon  s  s  nX vert Xi X lin (30) Ni (ϑ) · = s: Ys Yis i=1

Xis

Yis

and denote the physical coordinates of the nvert vertices of the polygon in the reference configuration, which lie on ΓN . Nodal basis functions Nilin (ϑ) defined over parameter space ϑ ∈ [−1, 1] linearly interpolate ΓN between vertices. Assuming that wherever polygon s crosses a cell edge, a vertex is placed, each polygon segment can be uniquely related to a finite cell, and the parameterization of s in terms of local cell coordinates {ξ, η} follows as  s  s  nX vert ξ ξ (31) Nilin (ϑ) · is s: s = η ηi i=1

ξis

ηis

Xis

Yis

according to Equation (6). The unit tangent vector and where and correspond to t in local cell coordinates for each polygon segment can be computed as   1 ∆ξ t= p (32) ∆ξ 2 + ∆η 2 ∆η

s s with ∆ξ=ξi+1 − ξis and ∆η=ηi+1 − ηis . At any point {ξ s , η s }T between polygon vertices, ˆ e can now be computed in analogy to the deformation dependent load the normal vector n

c 2000 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2000; 00:1–6

18

¨ D. SCHILLINGER, A. DUSTER AND E. RANK

(a) Voxels with bvox = 1.

(b) STL surface mesh.

Figure 13: Spatial voxel model of an aluminium foam sample.

implementation of p-FEM [74, 75]. The necessary derivatives of basis function N with respect to ϑ can be obtained from the directional derivative along t and a chain rule contribution being the ratio between the length of the segment in local cell coordinates and the length of the parameter space ϑ p ∂N ∂N ∆ξ 2 + ∆η 2 ∂N ∆ξ ∂N ∆η ∂ξ ∆ξ + ∂η ∆η N,ϑ = p = + (33) 2 2 2 ∂ξ 2 ∂η 2 ∆ξ + ∆η The force vector entries of each basis function Na can then be obtained by summing up the standard expression [49, 73, 74] over all nseg segments  s nseg Z 1 X −y,ϑ ext dϑ (34) fa = pˆ Na xs,ϑ −1 i=1

The deformation dependent tangent stiffness matrix of the linearized system Equation (23) follows from the linearization and discretization of Equation (29) (see [35, 49] for details). Analogous to the deformation dependent load implementation of p-FEM [74, 75], the contribution of each pair of basis functions Na and Nb to the load stiffness matrix follows again by summing up the standard expression [49, 73, 74] over all nseg segments   nseg Z 1 X 0 −1 dϑ (35) kab = pˆ Na Nb,ϑ 1 0 −1 i=1

A numerical example involving homogeneous Dirichlet and deformation dependent Neumann boundary conditions along circular boundaries will be shown in Section 6. 5.3. Geometrical models for efficient hp-d mesh generation

The fundamental advantage of FCM is the very simple and fast mesh generation irrespective of the geometric complexity involved. It is based on the disconnection of the finite cell mesh c 2000 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2000; 00:1–6

hp-d FCM FOR GEOMETRICALLY NONLINEAR PROBLEMS

(a) Complete hp-d mesh.

(d) Sub-cells of level k=2.

(b) Base mesh k=0.

19

(c) Sub-cells of level k=1.

(e) Sub-cells of level k=3.

(f ) Sub-cells of level k=4.

Figure 14: hp-d mesh for the foam example, consisting of a base mesh of high-order cells and several levels of sub-cells with a hierarchical overlay of linear basis functions.

from the geometry, which is represented implicitly by the change of material parameters at integration point level. The hp-d mesh consisting of structured finite cells and adaptive sub-cells with linear hierarchical basis functions can be established by the following simple algorithm 1. Traverse all sub-cells of the currently finest level k (start with the finite cells at k = 0) and query each Gauss point if it is in Ωphys or Ωfict . 2. If Gauss points of the same sub-cell are located in different domains (hence, a geometric boundary must be present), split the sub-cell into sub-cells of the next level k = k + 1. 3. Provide all new sub-cells with (p+1) n Gauss points, where n denotes the number of space dimensions. 4. Repeat this process, until a sufficient sub-cell depth k = m is reached. The fundamental component of this procedure is the location query, which is able to determine for an arbitrary point in space, if it is located in the physical or fictitious domain. Implicit geometrical models provide a natural way for point location queries, and two examples are reviewed here as a basis for fast hp-d mesh generation. c 2000 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2000; 00:1–6

20

¨ D. SCHILLINGER, A. DUSTER AND E. RANK

5.3.1. Implicit geometry representation An easy way of defining simple geometric boundaries is provided by inequalities in the form of a hypersurface, e.g. X2 Y2 Z2 + 2 + 2 ≤1 (36) 2 a b c Depending on parameters a,b,c, Equation (36) represents holes of different geometries, for example a circle of radius R (a=b=R; Z=0), a sphere of radius R (a=b=c=R) or an ellipsoid (a 6= b 6= c). If the inequality is satisfied at an integration point of coordinates {X, Y, Z}, it is located within Ωfict and the penalization of Equation (3) applies. 5.3.2. Voxel models Voxels are volume elements, aligned in a structured spatial grid in Cartesian directions, each of which is characterized by a bit code, which determines whether there is material (bvox =1) or a void (bvox =0) at the voxel location. Thus, the geometric boundaries of Ωphys are implicitly represented by the change of bvox from one voxel to the next. Voxel models are usually derived from computed tomography (CT) scans, for example in biomedical applications [76] or material science [77], but can also be generated from a B-rep geometry of a CAD model [78]. They provide an ideal geometrical basis for the generation of hp-d adaptive meshes, since they allow an easy and efficient mapping of an arbitrary global Gauss point position onto the corresponding voxel. The computational efficiency of the voxel approach for hp-d mesh generation is illustrated by the example of an aluminium foam sample. The foam is embedded in a rectangular domain of size 20 × 20 × 20 mm, covering aluminium and voids. The internal geometry is encoded by voxels with a resolution of 210 in each Cartesian direction, obtained from a CT scan∗ . Figure 13a shows all voxels of index bvox =1 in a coarsened resolution of 28 , while the STL surface of Figure 13b gives a more detailed idea of the foam geometry. The fully automated generation of the corresponding adaptive grid structure of Figure 14a, which consists of a base mesh of 5 × 5 × 5 finite cells and 4 levels of adaptive sub-cells with hierarchical basis functions, can be accomplished in only 56 seconds† . The main costs result from the loading of voxel information encoded by 210 bits, the generation of 111.917 adaptive sub-cells and about 210 location queries for the determination of α at all integration points. Figures 14b to 14f show the base mesh as well as the adaptive sub-cells of all levels. It can be easily observed that the adaptive aggregation of sub-cells around geometric boundaries increases quickly with k. 6. NUMERICAL EXAMPLES The hp-d adaptive FCM is implemented in the framework of nonlinear elasticity in principal directions as reviewed in Section 3, and its accuracy and computational efficiency are examined for a number of example problems. Our implementation is largely based on Sandia’s library framework Trilinos [79] and uses the direct solver Pardiso [80]. Its results are compared to overkill solutions derived with standard linear quadrilateral, quadratic triangular and quadratic

∗ Courtesy of IZFP Fraunhofer Institute for Non-Destructive Testing, Saarbr¨ ucken, Germany; http://www.izfp.fraunhofer.de † On a Intel(R) Core(TM)2 T5550 @ 1.83 GHz with 3 Gb RAM

c 2000 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2000; 00:1–6

21

hp-d FCM FOR GEOMETRICALLY NONLINEAR PROBLEMS

∆u ∆u

Ωphys

a

Ωfict B

r1 r2

b A

Parameters: a=2.0 b=0.45 t=0.1 ∆u=0.2 r1 =0.6 r2 =0.2 Ωphys : E=1.0 Ωfict : α E ν=0.2

Ωphys a

r3 r1

Ωfict r2

a

Parameters: a=2.0 ∆u=0.2 r1 =r2 =0.75 r3 =0.25 Ωphys : E=1.0 Ωfict : α E ν=0.2

a

a (a) 2D plate with ellipsoidal hole.

(b) 3D cube with ellipsoidal hole.

Figure 15: Ellipsoidal perforation problem: Parameters are width/height a; plate thickness t; radii ri ; Young’s modulus E; domains Ωphys and Ωf ict ; Poisson’s ratio ν; prescribed displacements ∆u

tetrahedral elements on conforming meshes, provided by the open-source nonlinear finite element code FlagShyp [81]. Conforming isoparametric mesh generation is accomplished by the meshers Visual Domesh [82] and Netgen [83], visualization is done with ParaView [84]. 6.1. Plate with ellipsoidal hole The first example problem consists of a 2D plate in plane stress, which is perforated by an ellipsoidal hole with a ratio between major and minor radius of 3:1. Material and geometric parameters as well as boundary conditions and the location of cutline A-B are given in Figure 15a. For the computation of the corresponding reference solution, the symmetry in geometry and boundary conditions is made use of to reduce the plate to 1/4 of the original system. The considered FlagShyp discretization with a mesh conforming to the geometric boundaries consists of 14.848 standard quadratic triangular elements and 59.679 degrees of freedom, taking into account the physical domain Ωphys only. Multiplying the resulting strain energy by 4 yields a reference for the complete system of Uex = 1.13618 · 10−3 . Convergence studies with different mesh sizes indicate that the given Uex is correct up to the 4th decimal, so that relative errors in terms of Equation (27) being larger than 1% can be reliably determined. For the FCM computations, the origin of the coordinate system is placed in the center of the ellipsoidal hole, so that its geometry can be described implicitly by the inequality (X/r1 )2 +(Y /r2 )2 ≤ 1.0. The complete domain is discretized by a structured FCM mesh of 4×4 high-order finite cells, which is complemented by an h-adaptive integration grid as shown in Figures 16a and 16b, respectively. To minimize the influence of a potential integration error, an integration quadtree with an overly large depth m=9 is applied throughout all computations. For the solution of the geometrically nonlinear system, the displacement load is cut into 5 equally sized increments, each of which requires 4 to 5 Newton-Raphson iterations to converge to a norm of the residual below 10−15 . Numerical experiments indicate that standard FCM has a lower bound of α=10−2 , beyond which the local uniqueness of the deformation map is not guaranteed. Due to the weak penalization, the character of geometric boundaries changes c 2000 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2000; 00:1–6

22

¨ D. SCHILLINGER, A. DUSTER AND E. RANK

(a) Base mesh of 4 × 4 finite cells.

(c) Von Mises stress on the deformed configuration, obtained with the base mesh of p=14 only.

(b) h-adaptive integration grid / hierarchical overlay mesh, depth m=6.

(d) Von Mises stress on the deformed configuration, obtained with the base mesh of p=14 and m=6 overlay levels.

Figure 16: Comparison of standard FCM and hp-d adaptive FCM for the perforated plate example.

to material interfaces, triggering strong oscillations in the stress solution (see Figure 16c). A more detailed view is given in Figure 17, which compares von Mises stresses along cutline A-B to the reference determined by the FlagShyp discretization. The convergence in strain energy for p-refinement on the given base mesh from p=1 to 20 is plotted in Figure 18 and exhibits a low algebraic rate of around q=0.23. This situation can be considerably improved by the introduction of hierarchical hp-d adaptivity. As explained in Sections 4 and 5, the high-order basis functions defined on the base mesh of Figure 16a are now complemented by an overlay of linear hierarchical hat functions defined on the adaptive integration grid, which is shown for maximum depth m=6 in Figure 16b. Varying the penalty parameter α in steps of one order of magnitude, numerical experiments show that the resulting hp-d adaptive FCM approach permits a locally stable nonlinear analysis with α=10−4 for all polynomial degrees between p=1 and 15. Due to the strengthened penalization and the better approximation of discontinuities around geometric boundaries, corresponding stress results are free of oscillations (see Figure 16d). The von Mises c 2000 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2000; 00:1–6

23

hp-d FCM FOR GEOMETRICALLY NONLINEAR PROBLEMS

102

0.16

Von Mises stress

0.14

Relative error in energy measure [%]

Reference Standard FCM, p=14, α= 10-2 hp-d FCM, p=14, m=6, α=10-4

0.12 0.10 0.08 0.06 0.04

0.0

0.2 0.4 0.6 0.8 Distance from bottom along cut line A-B

1.0

Figure 17: Von Mises stress plotted along the undeformed cut line A-B (see Figure 15a).

101

100

Standard FCM, α=10-2 hp-d FCM, m=p-1, α=10-2 hp-d FCM, m=p-1, α=10-4 hp-d FCM, min. m, α=10-4 102

103 Degrees of freedom

104

Figure 18: Convergence behavior for the plate with ellipsoidal hole.

stresses along cutline A-B plotted in Figure 17 illustrate the improved quality of the stress solution, which comes very close to the FlagShyp reference. Two different hp-d refinement strategies are examined in the following. First, the maximum depth m of hierarchical overlay levels is simultaneously increased with the polynomial degree p, which is similar to classical hp methods (see for example [85, 86]). Convergence in strain energy for a variation of p=1 to 10 and m=0 to 9 is plotted in Figure 18 for penalization parameters α=10−2 and 10−4 , respectively. The difference between the two curves illustrates the considerable impact of the strengthened penalization on the accuracy with respect to the reference solution. It can be further observed that for lower polynomial degrees with a limited number of overlay levels, large gains in accuracy can be achieved in comparison to the standard FCM results (max. q=1.34). However, the rate of convergence slows down with ascending number of levels m to an algebraic rate of q=0.28, since the amount of additional degrees of freedom from the overlay increases faster than the corresponding increase in accuracy. This effect can be reduced, if a minimum of overlay levels m is chosen, such that sufficient stabilization of the fictitious domain Ωfict is given to allow a penalization with α=10−4 . The corresponding convergence curve in Figure 18, which is obtained by increasing p=1 to 14 and adding m=2, 5 and 6 overlay levels at p=4,7 and 8, respectively, achieves a relative error of below 2% and a final rate of convergence of q=1.38 with only 4626 degrees of freedom. Under the assumption that standard FCM with α=10−2 could keep its rate of convergence and would not flatten earlier, the prolongation of the convergence curve in Figure 17 reveals that it would need about 3 orders of magnitude more degrees of freedom (around 4 million) to reach a comparably small level of error. 6.2. Cube with ellipsoidal hole The generalization of the 2D plate of the previous sub-section to 3 dimensions leads to a cube with an ellipsoidal hole, which exhibits an axis ratio of 3:3:1. A system sketch together with geometry, material and boundary conditions is given in Figure 15b. For the computation of the corresponding reference solution, the symmetry in geometry and boundary conditions is used to reduce the cube to 1/8 of the original system. The considered FlagShyp c 2000 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2000; 00:1–6

24

¨ D. SCHILLINGER, A. DUSTER AND E. RANK

(a) Base mesh of 3 × 3 × 3 finite cells (cut in the middle).

(c) Von Mises stress on the deformed configuration, obtained with finite cells of p=6 (plot of 1/8 of the symmetric system).

(b) Corresponding integration grid / hierarchical overlay mesh, depth m=5 (cut in the middle).

(d) Von Mises stress on the deformed configuration, obtained with the base mesh of p=6 and m=5 overlay levels.

Figure 19: Comparison standard FCM and hp-d adaptive FCM for the perforated cube example.

discretization with a mesh conforming to the geometric boundaries consists of 15.300 standard quadratic tetrahedral elements with 69.862 degrees of freedom, taking into account Ωphys only. Multiplying the resulting strain energy by factor 8 yields a reference for the complete problem of Uex =3.0664 · 10−2 . Convergence studies indicate that Uex is correct up to the 3rd decimal, so that relative errors in terms of Equation (27) being larger than 3% can be reliably determined. For the FCM computations, the geometry is described implicitly by the inequality (X/r1 )2 + (Y /r2 )2 + (Z/r3 )2 ≤ 1.0. The complete domain is discretized by a structured FCM grid of 3 × 3 × 3 finite cells, complemented by an h-adaptive integration grid. They are shown in Figures 19a and 19b, respectively, where half of the symmetric discretization is cut away to uncover the refinement around the ellipsoid in the center. Throughout all computations, an octree of integration sub-cells of total depth m=5 is applied. For the solution of the geometrically nonlinear system, the displacement load is cut into 5 equally sized increments, c 2000 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2000; 00:1–6

25

hp-d FCM FOR GEOMETRICALLY NONLINEAR PROBLEMS

Relative error in energy measure [%]

80

Ω2 40

ΓD

Ωphys ΓN Ω1

20

pˆ 10 6 2 10

rout

Standard FCM, α=10-2 hp-d FCM, m=p-1, α=10-2 hp-d FCM, m=p-1, α=10-4 10 3 10 4 Degrees of freedom

rin

Parameters: ΓN : pˆ=0.25 ΓD : u=[0.0 0.0]T rin =0.25 rout =1.0 Ωphys : E=1.0 Ω1 : E1 =10−4 Ω2 : E2 =104 ν=0.0 t=1.0

10 5

Figure 20: Convergence behavior for the 3D cube with ellipsoidal hole.

Figure 21: Ring under internal pressure pˆ and fixed outer displacements.

Support of overlay elements of finest level Polygon segments Polygon vertices

Figure 22: Adaptive integration sub-cells / Hierarchical overlay mesh for the ring example with depth m=7. The zoom illustrates the approximation of the Neumann boundary by polygon segments.

each of which requires 3 to 4 Newton-Raphson iterations to converge to a norm of the residual below 10−10 . For the standard FCM formulation operating on the base mesh alone, the lower bound of the penalty parameter is α=10−2 , leading to oscillations in stresses as depicted in Figure 19c. The corresponding convergence in strain energy is plotted in Figure 20 for ascending polynomial degree p=1 to 6, exhibiting an average rate of q=0.22. In analogy to the previous sub-section, the introduction of hp-d adaptivity leads to a considerable improvement. The lower bound of α can be reduced to 10−4 already with a small number m of overlay levels, thus strengthening the penalization in Ωfict . In addition, the better approximation of discontinuities around geometric boundaries contributes to a more accurate localization of stress concentrations around geometric boundaries and to a considerable reduction of stress oscillations (see Figure 19d). The corresponding convergence in strain energy is shown Figure 20, obtained by increasing p from 1 to 6 on the base mesh and simultaneously increasing the number of hierarchical overlay levels m=p − 1. It indicates that a considerable error reduction can be achieved with the first three to four overlay levels, c 2000 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2000; 00:1–6

26

¨ D. SCHILLINGER, A. DUSTER AND E. RANK

Figure 23: Von Mises stress plotted on the deformed configuration, obtained with standard FCM on a 6 × 6 base mesh of p=8 (left) and the hp-d FCM, using p=8 and m=7 overlay levels (right).

whereas further refinement implies an overly large increase in additional degrees of freedom, not paid off by the corresponding decrease in error. For the present example, this is illustrated by comparison of the rates of convergence of the first two refinement steps, i.e. q=0.62 and q=0.57, respectively, with the final rate q=0.14. This characteristic solution behavior makes the hp-d FCM a suitable tool for costly computations of large geometries in 3D, where more than a few overlay levels would quickly exceed the hardware resources of common workstations. 6.3. Ring under internal pressure The 2D example of a circular ring in plane strain under internal pressure shown in Figure 21 illustrates the capabilities of the methods discussed in Section 5.2 for the incorporation of non-conforming boundary conditions in the framework of the hp-d FCM. The geometry of 2 . the domain enclosed by ΓD and ΓN is described implicitly by rin ≤ X 2 + Y 2 ≤ rout The FCM discretization consists of a structured base mesh of 6 × 6 cells with polynomial degree p=8 and an overlay mesh of depth m=7 as shown in Figure 22, also used for adaptive integration. Dirichlet boundary conditions along ΓD are realized by fixed displacements at the outer boundaries of the finite cell mesh, which are carried on to ΓD due to the increased stiffness E2 in Ω2 , whereas the pressure boundary condition is taken into account by integration over ΓN . Figure 22 also illustrates the approximation of the ΓN by polygon segments in the sense of Equations (30) and (31). Von Mises stress plots are shown in Figure 23, obtained with standard FCM using only the base mesh of p=8, and with the hp-d adaptive FCM adding basis functions corresponding to m=7 overlay levels. The contribution of the deformation dependent load to the external force vector and the stiffness matrix Equations (34) and (35) are obtained by integration over 12.068 polygon segments. The pressure load pˆ is divided into 25 equal increments, each of which requires 3 to 4 Newton-Raphson iterations to converge to a norm of the residual below 10−8 . Whereas stress results at the Neumann boundary are accurate in both approaches, the standard FCM discretization without overlay yields largely inaccurate stresses with strong oscillations around the Dirichlet boundary. The hp-d FCM considerably improves the stress solution at ΓD by limiting oscillations to the close vicinity of the geometric boundaries, leading to the correct circular stress pattern. The example illustrates that hp-d adaptivity helps to improve the representation of unfitted boundary conditions. c 2000 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2000; 00:1–6

hp-d FCM FOR GEOMETRICALLY NONLINEAR PROBLEMS

27

Figure 24: Von Mises stress plotted on the undeformed configuration.

6.4. Open cell aluminium foam Metal foams provide high stiffness at reduced weights, and are therefore frequently used for lightweight structures in automotive and aerospace applications [87]. The hp-d adaptive Finite Cell Method is applied to simulate a nonlinear compression test of an aluminium foam sample, whose geometry is given by the voxel model shown in Figure 13a. The corresponding hp-d mesh generation process has been described in detail in Section 5.3. As shown in Figure 14, the sample is discretized by a 5 × 5 × 5 base mesh of finite cells of polynomial degree p=6, which is complemented by an octree-based adaptive integration structure consisting of subcells up to level k=4. The aluminium foam is characterized by Young’s modulus E=70.000 N/mm2 and Poisson’s ratio ν=0.35. Modeling the influence of a testing machine, vertical top displacements are gradually increased to -1mm (5% compressive deformation), while horizontal displacements at the top surface and all displacements at the bottom surface are completely fixed. The geometrically nonlinear system is solved by applying the displacement load on top in 4 increments, each of which requires 3 to 4 Newton-Raphson iterations to converge below a norm of the residual of 10−8 . Due to hp-d stabilization, the penalization of the fictitious domain can be gradually strengthened from α=10−3 for analysis with the base mesh only down to α=10−5 for analysis with the full number of overlay levels m=4. The resulting von Mises stresses computed with all overlay levels m=4 and a total number of 207.488 degrees of freedom are displayed in Figure 24. Stresses are free of large-scale oscillations and exhibit accurate localization of stress concentrations at the convex sides of the foam members, which agrees well with engineering experience (see also the zoom in Figure 25). The convergence of the equivalent vertical top force obtained by integrating the normal stress component over the top surface, is illustrated in Figure 26 with respect to an increase in overlay levels and constant polynomial degree p=6. It can be observed that the addition of only a few hierarchical overlays around geometric boundaries in combination with the decrease of α considerably improves the prediction accuracy of the top force. Whereas the hp-d FCM in the c 2000 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2000; 00:1–6

¨ D. SCHILLINGER, A. DUSTER AND E. RANK

28

Equivalent vertical top load [kN]

-36

-3

m=0; α=10 -3 m=1; α=10

-34

-32

m=2; α=10

-30

-3

m=3; α=10 -28

-4 -5

m=4; α=10

-26 10.000

25.000 50.000 100.000 Degrees of freedom

250.000

Figure 26: Convergence of vertical top force.

Figure 25: Zoom into part of the structure.

last refinement step from m=3 to m=4 exhibits a flattening of the convergence curve towards a final value of around 27.1 kN, the standard FCM approach, using the sub-cells for integration purposes only, shows a substantial deviation of about 30%. The present foam example indicates the effectiveness of the hp-d approach for large-scale problems of very complex geometry, which can be accurately solved at a manageable computational effort, while conforming mesh generation problem is completely circumvented. Besides the direct simulation of complex structures, possible applications of the hp-d adaptive FCM comprise homogenization and FE2 approaches for repeated nonlinear RVE computations (see for example [88, 89]).

7. hp-d ADAPTIVITY FOR INTERFACE AND SINGULAR PROBLEMS To motivate further development in this direction, the potential of hierarchical hp-d adaptivity for problems beyond geometric nonlinearity is briefly addressed. Therefore, two popular benchmark problems are examined, namely the material interface problem of a 2D plate with a circular inclusion and the singular problem of the L-shaped domain. ∆u Ω1 Ω2 a

B

r

b

Parameters: a=2.0 b=0.45 t=0.1 ∆u=0.3 r=0.6 Ω1 : E=1.0 Ω2 : 0.5 E ν=0.2

A a Figure 27: Plate with circular inclusion. c 2000 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls

ΓN ΓN

a A

ΓN

Γ0

a

Parameters: a=1.0 t=1.0 E=1.0 ν=0.3 ΓN : σ·n6=0 Γ0 : σ·n=0

Γ0 ΓN a

a

Figure 28: The L-shaped domain. Int. J. Numer. Meth. Engng 2000; 00:1–6

29

hp-d FCM FOR GEOMETRICALLY NONLINEAR PROBLEMS

k=5 k=4 k=3 k=2

k=1

k=0

(a) Adaptive integration grid / hierarchical overlay mesh, depth m=9 on a base mesh of 4 × 4.

(c) Von Mises stress on the deformed configuration, obtained with the base mesh of p=10 only.

(b) Structure of hierarchical overlay levels, illustrated for levels k=0 to 5.

(d) Von Mises stress on the deformed configuration, obtained with the base mesh of p=10 and m=9 overlay levels.

Figure 29: Standard FCM and hp-d adaptive FCM for the circular inclusion problem.

7.1. Plate with circular inclusion In Figure 27, material and geometric parameters as well as boundary conditions and the location of cutline A-B for a plate in plane stress with a circular inclusion are given. Corresponding solution fields exhibit discontinuities along the material interface, i.e. kinks in displacements and jumps in stresses and strains. Assuming a geometrically nonlinear formulation as in the previous sections, the FCM implementation can be analogously applied to this class of problems with the only difference that the fictitious domain Ωphys is replaced by a second fully valid physical domain Ω2 with substantial stiffness. The complete domain is again discretized by a structured base mesh of 4 × 4 high-order finite cells (see Figure 16a), which is complemented by a quadtree based adaptive integration structure as illustrated in Figure 29a and 29b. The symmetry in geometry and boundary conditions is again made use of to reduce the plate to 1/4 of the original system for the computation of the corresponding reference solution with FlagShyp. The discretization considered consists of a mesh conforming to the c 2000 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2000; 00:1–6

¨ D. SCHILLINGER, A. DUSTER AND E. RANK

30

0.22 Reference Standard FCM, p=10 hp-d FCM, p=10, m=9

Relative error in energy measure [%]

Von Mises stress

0.20 0.18 0.16 0.14 0.12 0.10 0.0

0.2 0.4 0.6 0.8 Distance from bottom along cut line A-B

100

Standard FCM hp-d FCM, m=p-1 10-1 1 10

1.0

Figure 30: Von Mises stress plotted along undeformed cutline A-B (see Figure 27).

101

102

103 104 Degrees of freedom

105

Figure 31: Convergence behavior for the plate with circular inclusion.

interface with 106.674 standard linear quadrilaterals and 213.681 degrees of freedom, leading to a reference strain energy of Uex = 4.330559 · 10−3 for the complete system. Convergence studies with different mesh sizes indicate that the given Uex is correct up to the 4th decimal. The standard FCM approach leads to excessive stress oscillations in the cells cut by material interfaces, since the high-order basis alone is not able to accurately represent the discontinuous solution, which is illustrated for polynomial degree p=10 in Figure 29c. The hp-d adaptive FCM offers a better approximation of discontinuities along material interfaces. The corresponding stress plot of Figure 29d demonstrates that an addition of m=9 overlay levels to the base mesh of polynomial degree p=10 leads to an effective reduction of oscillations and an improved localization of stress concentrations along vertical inclusion boundaries. This observation is further corroborated by the von Mises stress plotted along cutline A-B in Figure 30. Whereas the standard FCM approach yields strong oscillations over the whole domain, the hp-d adaptive FCM solution is free of large-scale oscillations, reducing the oscillatory stress behavior to the close vicinity of the interface location. Figure 31 compares the strain energy convergence in terms of Equations (27) and (28) for the standard and hp-d adaptive FCM schemes. While pure p-refinement on the base mesh achieves only a low algebraic rate of convergence of around q=0.25, hp-d -refinement complementing the base mesh by m=p-1 hierarchical overlays of linear functions is able to improve the convergence rate to q=0.55. 7.2. The L-shaped domain The L-shaped domain shown in Figure 28 exhibits a corner singularity at the intersection point A of the two reentrant edges. For a detailed discussion of its mathematical characteristics and hp-refinement strategies for its efficient numerical solution, see [27, 28, 85, 86]. For the specific problem considered in the following, plane strain linear elasticity is assumed (see for example [2, 3]), for which an analytical solution of the strain energy exists in the form A1 a2λ1 (37) E if the following traction boundary conditions are prescribed along the outer boundaries ΓN Uex = 4.15454423

σx = A1 λ1 rλ1 −1 [(2 − Q1 (λ1 + 1)) cos ((λ1 − 1)θ) − (λ1 − 1) cos ((λ1 − 3)θ)] c 2000 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls

(38)

Int. J. Numer. Meth. Engng 2000; 00:1–6

31

hp-d FCM FOR GEOMETRICALLY NONLINEAR PROBLEMS

(a) Base mesh of 12 high-order elements.

(b) Overlay mesh of linear hierarchical functions, depth m=8.

Figure 32: Conforming hp-d adaptive discretization of the L-shaped domain.

Relative error in energy measure [%]

50.0

10.0

1.0 0.5 1 10

Figure 33: Von Mises stress (p=12, m=11).

h-extension, p=2 p-extension hp-d extension p-extension (geometr. graded mesh)

10 2 10 3 Degrees of freedom

Figure 34: Convergence for the L-shaped domain (graded mesh results from [27, 28]).

σy = A1 λ1 rλ1 −1 [(2 + Q1 (λ1 + 1)) cos ((λ1 − 1)θ) + (λ1 − 1) cos ((λ1 − 3)θ)] τxy = A1 λ1 r

λ1 −1

10 4

[(λ1 − 1) sin ((λ1 − 3)θ) + Q1 (λ1 + 1) sin ((λ1 − 1)θ)]

(39) (40)

where r and θ denote polar coordinates with origin at the corner singularity point A. The values of the constants involved are the smallest eigenvalue λ1 =0.544483737, the generalized stress intensity factor A1 =1.0 and Q1 =0.543075579. The reentrant edges Γ0 are stress-free and constants a and E are geometric and material parameters (see Figure 28). A detailed derivation of Equations (37) to (40) can be found in [27]. The L-shaped domain is discretized by a base mesh of 12 conforming high-order finite elements as shown in Figure 32a. Combining the principles illustrated in Figures 8a and 11, the base mesh is complemented by an overlay mesh of hierarchical hat functions, which realizes an adaptive h-refinement towards the singular point A as illustrated in Figure 32b. In each overlay level k, all sub-cells touching point A are refined by partitioning them into four new sub-cells of the next level k+1. Following Figure 12a, all those newly created sub-cell nodes are c 2000 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2000; 00:1–6

32

¨ D. SCHILLINGER, A. DUSTER AND E. RANK

equipped with linear hat functions of level k+1 that are either in the center of the partitioned cell or lie on a boundary that touches other sub-cells of level k+1. The resulting von Mises stress obtained with polynomial degree p=12 of the base mesh and m=11 adaptive overlay levels is plotted in Figure 33, illustrating the considerable localization of stresses in point A. The stress solution further away from the location of the singularity is smooth and can be accurately fitted by the high-order base mesh. To estimate the convergence properties of hierarchical hp-d adaptivity, an hp-d refinement strategy is applied, which increases the polynomial degree p on the given base mesh and the number of hierarchical overlay levels m=p-1 simultaneously. In Figure 34, the resulting convergence behavior is compared to those obtained by classical h-refinement with p=2, by classical p-refinement on the base mesh only and by p-refinement on a geometric mesh, which is adaptively graded towards the location of the singularity. The geometry of the graded mesh and the corresponding strain energy convergence are reported by Szab´o and Babuˇska in [27]. Whereas h-refinement, bisecting the characteristic element width h in each refinement step, achieves a rate of convergence of around q=0.2775, all other refinement strategies involving the elevation of the polynomial degree p exhibit a very similar final rate between q=0.5413 and q=0.5518. These results are in good accordance with corresponding theoretical estimates [27], predicting qest =0.272 and qest =0.544 for h- and p-refinement, respectively. The performance of the different p-refinement strategies is largely influenced by their pre-asymptotic rates, which can considerably reduce the error level before the final rate tunes in. While the efficiency of the geometrically graded mesh with maximum rates of q=1.36 cannot be reached, the hp-d adaptive refinement achieves a distinct pre-asymptotic exponential convergence with a maximum rate of q=0.84, which is considerably better than pure p-refinement.

8. SUMMARY AND CONCLUSIONS In the present paper, hierarchical hp-d adaptivity has been introduced and combined with the standard Finite Cell Method (FCM) for the accurate solution of geometrically nonlinear fictitious domain problems in solid mechanics. First, the loss of local uniqueness of the deformation map within the fictitious domain for small penalty parameters α has been identified as the main obstacle that prevents a straightforward extension of the standard FCM to geometric nonlinearity. Increasing α can stabilize the deformation map, but considerably weakens the penalization of the fictitious domain. Therefore, the nature of FCM changes to an interface problem, entailing the typical phenomena of high-order schemes with inter-element discontinuities, i.e. excessive stress oscillations and a decay of convergence in strain energy to a low algebraic rate. Hierarchical hp-d adaptivity synergetically uses the h-adaptivity of the adaptive integration scheme set up by the standard FCM formulation around geometric boundaries, thus introducing a highly adaptive overlay of linear basis functions at a limited number of additional degrees of freedom. The hp-d overlay stabilizes the deformation map in the fictitious domain part of cells cut by geometric boundaries, thus permitting a considerable decrease of the penalty parameter α in comparison to standard FCM. Due to the reenforced penalization of the fictitious domain in conjunction with the better approximation of discontinuities along geometric boundaries, oscillations can be considerably reduced and high rates of convergence can be achieved. In addition, the hp-d adaptive FCM preserves simple and fast mesh generation for complex geometries, constituting the key advantage of c 2000 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2000; 00:1–6

hp-d FCM FOR GEOMETRICALLY NONLINEAR PROBLEMS

33

the original FCM concept. These benefits in combination with its computational efficiency have been demonstrated for a range of benchmark problems in 1, 2 and 3 dimensions. In particular, it has been shown by the nonlinear simulation of a metal foam sample that the hp-d adaptive FCM can handle very complex geometries with high accuracy at a manageable computational cost, since effective hp-d stabilization can be already achieved by adding only a few hierarchical overlay levels. Finally, the potential of hierarchical hp-d adaptivity for the application to material interface and singular problems has been outlined, giving an incentive for future research in this direction.

ACKNOWLEDGEMENTS

The first author gratefully acknowledges financial support from the Centre of Advanced Computing (MAC) and the International Graduate School of Science and Engineering (IGSSE) within the Excellence Initiative of the German Federal Government at the Technische Universit¨ at M¨ unchen. The authors would like to thank J. Frisch and R.-P. Mundani for their valuable support during processing of the large voxel data sets.

REFERENCES 1. Bathe K-J. Finite Element Procedures. Prentice Hall: Upper Saddle River, 1996. 2. Hughes TJR. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis. Dover: New York, 2000. 3. Felippa CA. Introduction to finite element methods. University of Colorado, Boulder. http://www.colorado.edu/engineering/CAS/courses.d/IFEM.d/Home.html [10 October 2010]. 4. Cottrell JA, Hughes TJR, Bazilevs Y. Isogeometric Analysis: Toward Integration of CAD and FEA. Wiley: New York, 2008. 5. Hughes TJR, Cottrell JA, Bazilevs Y. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering 2005; 194:4135–4195. 6. Bazilevs Y, Calo VM, Cottrell JA, Evans JA, Hughes TJR, Lipton S, Scott MA, Sederberg TW. Isogeometric analysis using T-splines. Computer Methods in Applied Mechanics and Engineering 2010; 199:229–263. 7. Del Pino S, Pironneau, O. A fictitious domain based general pde solver. In: Numerical methods for scientific computing: Variational problems and applications, Kuznetsov Y, Neittanmaki P, Pironneau O (eds). CIMNE:Barcelona, 2003. 8. Glowinski R, Kuznetsov Y. Distributed Lagrange multipliers based on fictitious domain method for second order elliptic problems. Computer Methods in Applied Mechanics and Engineering 2007; 196:1498–1506. 9. Rami` ere I, Angot P, Belliard M. A fictitious domain approach with spread interface for elliptic problems with general boundary conditions. Computer Methods in Applied Mechanics and Engineering 2007; 196:766–781. 10. Bishop J. Rapid stress analysis of geometrically complex domains using implicit meshing. Computational Mechanics 2003; 30:460–478. 11. L¨ ohner R, Cebrala RJ, Camellia FE, Appanaboyinaa S, Baum JD, Mestreau EL, Soto OA. Adaptive embedded and immersed unstructured grid techniques. Computer Methods in Applied Mechanics and Engineering 2008; 197:2173–2197. 12. Baaijens FPT. A fictitious domain/mortar element method for fluid-structure interaction. International Journal for Numerical Methods in Fluids 2001; 35:743–761. 13. Farhat C, Hetmaniuk U. A fictitious domain decomposition method for the solution of partially axisymmetric acoustic scattering problems. Part I: Dirichlet boundary conditions. International Journal for Numerical Methods in Engineering 2002; 54:1309–1332. 14. Haslinger J, Kozubek T, Kunisch K, Peichl G. Shape Optimization and Fictitious Domain Approach for Solving Free Boundary Problems of Bernoulli Type. Computational Optimization and Applications 2003; 26:231–251. 15. Burman E, Hansbo P. Fictitious domain finite element methods using cut elements: I. A stabilized Lagrange multiplier method. Computer Methods in Applied Mechanics and Engineering 2010; 199:2680–2686. 16. Burman E, Hansbo P. Fictitious domain finite element methods using cut elements: II. A stabilized Nitsche method. Preprint submitted to Computer Methods in Applied Mechanics and Engineering 2010. c 2000 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2000; 00:1–6

34

¨ D. SCHILLINGER, A. DUSTER AND E. RANK

17. Sukumar N, Chopp DL, Mo¨ es N, Belytschko T. Modeling holes and inclusions by level sets in the extended finite-element method. Computer Methods in Applied Mechanics and Engineering 2001; 190:6183–6200. 18. Gerstenberger A, Wall WA. An eXtended Finite Element Method/Lagrange multiplier based approach for fluid-structure interaction. Computer Methods in Applied Mechanics and Engineering 2008; 197:1699–1714. 19. Haslinger J, Renard, Y. A New Fictitious Domain Approach Inspired by the Extended Finite Element Method. SIAM Journal on Numerical Analysis 2009; 47:1474–1499. 20. Becker R, Burman E, Hansbo P. A hierarchical NXFEM for fictitious domain simulations. Preprint submitted to International Journal for Numerical Methods in Engineering 2010. 21. Bastian P, Engwer C. An unfitted finite element method using discontinuous Galerkin. International Journal for Numerical Methods in Engineering 2009; 79:1557–1576. 22. Lui S. Spectral domain embedding for elliptic PDEs in complex domains. Journal of Computational and Applied Mathematics, 2009; 225:541–557. 23. Parussini L, Pediroda V. Fictitious Domain approach with hp-finite element approximation for incompressible fluid flow. Journal of Computational Physics, 2009; 228:3891–3910. 24. Allaire G, Jouve F, Toader AM. Structural optimization using sensitivity analysis and a level-set method. Journal of Computational Physics, 2004; 194:363–393. 25. Parvizian J, D¨ uster A, Rank E. Finite Cell Method: h- and p- extension for embedded domain methods in solid mechanics. Computational Mechanics 2007; 41:122–133. 26. D¨ uster A, Parvizian J, Yang Z, Rank E. The Finite Cell Method for Three-Dimensional Problems of Solid Mechanics. Computer Methods in Applied Mechanics and Engineering 2008; 197:3768–3782. 27. Szab´ o B, Babuˇska, I. Finite Element Analysis. Wiley: New York, 1991. 28. Szab´ o B, D¨ uster A, Rank, E. The p-version of the finite element method. In: Encyclopedia of Computational Mechanics Vol. 1, Stein E, de Borst R, Hughes TJR (eds). Wiley: Chichester, 2004; 119–139 (Chapter 5). 29. Rank E, D¨ uster A, Schillinger D, Yang Z. The Finite Cell Method: High order simulation of complex structures without meshing. In: Computational Structural Engineering, Yuan Y, Cui J, Mang H (eds). Springer: Heidelberg, 2009; 87–92. 30. Ruess M, Tal D, Trabelsi N, Yosibash Z, Rank E. The finite cell method for bone simulations: Verification and validation. Accepted for publication in Biomechanics and Modeling in Mechanobiology 2011. 31. Vos PEJ, van Loon R, Sherwin SJ. A comparison of fictitious domain methods appropriate for spectral/hp element discretisations. Computer Methods in Applied Mechanics and Engineering 2008; 197:2275–2289. 32. Parvizian J, D¨ uster A, Rank E. Toplogy Optimization Using the Finite Cell Method. Accepted for publication in Optimization and Engineering, 2011. 33. Rank E, Kollmannsberger S, Sorger C, D¨ uster A. Shell Finite Cell Method: A High Order Fictitious Domain Approach for Thin-Walled Structures. submitted to Computer Methods in Applied Mechanics and Engineering, 2011; 34. D¨ uster A, Hartmann S, Rank E. p-fem applied to finite isotropic hyperelastic bodies. Computer Methods in Applied Mechanics and Engineering 2003; 192:5147–5166. 35. Bonet J, Wood R. Nonlinear Continuum Mechanics for Finite Element Analysis. Cambridge University Press: Cambridge, 2008. 36. de Souza Neto EA, Peri´ c D, Owen DRJ. Computational Methods for Plasticity: Theory and Applications. Wiley: New York, 2008. 37. Rank E. Adaptive remeshing and h-p domain decomposition. Computer Methods in Applied Mechanics and Engineering, 1992; 101:299–313. 38. Rank E, Krause R. A multiscale finite-element-method. Computers & Structures, 1997; 64:139–144. 39. Krause R, Rank E. Multiscale computations with a combination of the h- and p-versions of the finiteelement method. Computer Methods in Applied Mechanics and Engineering, 2003; 192:3959–3983. 40. D¨ uster A, Niggl E, Rank E. Applying the hp-d version of the FEM to locally enhance dimensionally reduced models. Computer Methods in Applied Mechanics and Engineering, 2007; 196:3524–3533. 41. Schillinger D, Kollmannsberger S, Mundani R-P, Rank E. The finite cell method for geometrically nonlinear problems of solid mechanics. IOP Conference Series: Material Science and Engineering, 2010; 10:012170. 42. Schillinger D, Rank E. An unfitted hp-adaptive finite element method based on hierarchical B-splines for interface problems of complex geometry. Accepted for publication in Computer Methods in Applied Mechanics and Engineering, 2011; 43. Gerstenberger A, Wall WA. An embedded Dirichlet formulation for 3D continua. International Journal for Numerical Methods in Engineering 2010; 82:537–563. 44. Embar A, Dolbow J, Harari I. Imposing Dirichlet boundary conditions with Nitsche’s method and splinebased finite elements. International Journal for Numerical Methods in Engineering 2010; 83:877–898. 45. Samet H. The design and analysis of spatial data structures. Addison Wesley: Reading, 1989. 46. Zohdi T, Wriggers P. Aspects of the computational testing of the mechanical properties of microheterogeneous material samples. International Journal for Numerical Methods in Engineering 2001; 50:2573–2599. c 2000 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2000; 00:1–6

hp-d FCM FOR GEOMETRICALLY NONLINEAR PROBLEMS

35

47. Zienkiewicz OC, Taylor RL. The Finite Element Method Vol. 1: The Basis. Butterworth-Heinemann: Oxford, 2000. 48. Belytschko T, Liu WK, Moran B. Nonlinear Finite Elements for Continua and Structures. Wiley: New York, 2006. 49. Wriggers P. Nonlinear Finite Element Methods. Springer: Berlin, 2008. 50. Ibrahimbegovic A. Nonlinear Solid Mechanics: Theoretical Formulations and Finite Element Solution Methods, Springer: Heidelberg, 2009. 51. MATLAB. MATLAB user’s guide version 7.6. The MathWorks Inc.: Natick, Mass., 2008. 52. N¨ ubel V, D¨ uster A, Rank E. An rp-adaptive finite element method for elastoplastic problems. Computational Mechanics 2007; 39:557–574. 53. Krause R. Multiscale Computations with a Combined h- and p-version of the Finite-Element Method. PhD thesis, Fach Numerische Methoden und Informationsverarbeitung, Universit¨ at Dortmund, 1996. 54. D¨ uster A. High order finite elements for three-dimensional, thin-walled nonlinear continua. PhD thesis, Lehrstuhl f¨ ur Bauinformatik, Technische Universit¨ at M¨ unchen. Shaker: Aachen, 2001. 55. Yserantant H. On the multi-level splitting of finite element spaces. Numerische Mathematik 1986; 49:379– 412. 56. Krysl P, Grinspun E, Schr¨ oder P. Natural hierarchical refinement for finite element methods. International Journal for Numerical Methods in Engineering 2003; 56:1109–1124. 57. Bungartz HJ, Griebel M. Sparse grids. Acta Numerica 2004; 13:147–269. 58. Krause R, M¨ ucke R, Rank E. hp-version finite elements for geometrically non-linear problems. Communications in Numerical Methods in Engineering 1995; 11:887–897. uster A, Hartmann S, Rank E. p-fem applied to finite isotropic hyperelastic bodies. Computer Methods 59. D¨ in Applied Mechanics and Engineering 2003; 192:5147–5166. 60. Dong S, Yosibash Z. A parallel spectral element method for dynamic three-dimensional nonlinear elasticity problems. Computers & Structures 2009; 87:59–72. 61. Schillinger D, Ruess M, D¨ uster A, Rank E. Large deformation analysis with the p- and B-spline versions of the Finite Cell Method. Part I: An improved geometrically nonlinear FCM formulation. Preprint submitted to Computational Mechanics 2011. 62. Schillinger D, Ruess M, Zander N, Bazilevs Y, D¨ uster A, Rank E. Large deformation analysis with the p- and B-spline versions of the Finite Cell Method. Part II: Unfitted Dirichlet boundary conditions, severe mesh distortion and application to complex voxel based geometries. Preprint submitted to Computational Mechanics 2011. 63. Ibrahimbegovic A, Gruttmann F. A Consistent Finite Element Formulation of Nonlinear Membrane Shell Theory With Particular Reference to Elastic Rubberlike Material. Finite Element in Analysis and Design 1993; 12:75–86. 64. Gharzeddine F, Ibrahimbegovic A. Incompatible Mode Method for Finite Deformation QuasiIncompressible Elasticity. Computational Mechanics 2000; 24:419–425. 65. Suri M. Analytical and computational assessment of locking in the hp finite element method. Computer Methods in Applied Mechanics and Engineering 1996; 133 :347–371. 66. Chilton L, Suri M. On the selection of a locking-free hp element for elasticity problems. International Journal for Numerical Methods in Engineering 1997; 40:2045–2062. 67. Heisserer U, Hartmann S, D¨ uster A, Yosibash Z. On volumetric locking-free behaviour of p-version finite elements under finite deformations. Communications in Numerical Methods in Engineering 2008; 24:1019– 1032. 68. Hansbo A, Hansbo P. An unfitted finite element method, based on Nitsche’s method, for elliptic interface problems. Computer Methods in Applied Mechanics and Engineering 2002; 191:537–552. 69. Fernandez-Mendez S, Huerta A. Imposing essential boundary conditions in mesh-free methods. Computer Methods in Applied Mechanics and Engineering 2004; 193:1257–1275. 70. Fritz A, H¨ ueber S, Wohlmuth BI. A comparison of mortar and Nitsche techniques for linear elasticity. Calcolo 2004; 41:115–137. 71. Schweizerhof K, Ramm E. Displacement dependent pressure loads in nonlinear finite element analyses. Computers & Structures 1984; 18:1099–1114. 72. Simo JC, Taylor RL, Wriggers P. A note on finite-element implementation of pressure boundary loading. Communications in Applied Numerical Methods 1991; 7:513–525. 73. Mok DP, Wall WA, Bischoff M, Ramm E. Algorithmic aspects of deformation dependent loads in non-linear static finite element analysis. Engineering Computations 1999; 16:601–618. 74. Yosibash Z, Hartmann S, Heisserer U, D¨ uster A, Rank E, Szanto M. Axisymmetric pressure boundary loading for finite deformation analysis using p-FEM. Computer Methods in Applied Mechanics and Engineering 2007; 196:1261–1277. 75. Heisserer U. High-order finite elements for material and geometrical nonlinear finite strain problems. PhD thesis, Cair for Computation in Engineering, Technische Universit¨ at M¨ unchen. Shaker: Aachen, 2008. c 2000 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2000; 00:1–6

36

¨ D. SCHILLINGER, A. DUSTER AND E. RANK

76. Yosibash Z, Padan R, Joskowicz L, Milgrom C. A CT-based high-order finite element analysis of the human proximal femur compared to in-vitro experiments. Journal of Biomedical Engineering 2007; 129:297–309. ¨ 77. Fiedler T, Sol´ orzano E, Garcia-Moreno F, Ochsner A, Belova IV, Murch GE. Computed tomography based finite element analysis of the thermal properties of cellular aluminium. Materialwissenschaften und Werkstofftechnik 2009; 40: 139–143. 78. Wenisch P, Wenisch O. Fast octree-based voxelization of 3D boundary representation-objects. Technical report, Lehrstuhl f¨ ur Bauinformatik, Technische Universit¨ at M¨ unchen, 2004. 79. Trilinos Version 10.2. Sandia National Laboratories, http://trilinos.sandia.gov/, Los Alamos, NM, 2010. 80. PARDISO Solver Project. Direct solver developed by O. Schenk, K. G¨ artner et al., http://www.pardisoproject.org, 2010. 81. FlagShyp 08 Version 2.30. Nonlinear FE solver developed by J. Bonet and R. Wood. http://www.flagshyp.com, 2008. 82. Visual DoMesh 2008 Version 1.1. Mesh generator developed by C. Sorger, Chair for Computation in Engineering, Technische Universit¨ at M¨ unchen, 2010. 83. Netgen Version 4.9.13. Tetrahdral mesh generator developed by J. Sch¨ oberl, http://sourceforge.net/projects/netgen-mesher, 2010. 84. ParaView Version 3.8.1. Open-source scientific visualization package, http://www.paraview.org, Kitware Inc., Clifton Park, NY, 2010. ˇ ın P, Segeth K, Doleˇ 85. Sol´ zel I. Higher-order Finite Element Methods. Chapman&Hall/CRC:Boca Raton, 2004. 86. Demkowicz L. Computing with hp-Adaptive Finite Elements; Volume 1: One and Two Dimensional Elliptic and Maxwell Problems. Chapman&Hall/CRC:Boca Raton, 2006. 87. Banhart J. Manufacture, characterization and application of cellular metals and metal foams. Progress in Material Science 2001; 46:559–632. 88. Zohdi T, Wriggers P. Introduction to Computational Micromechanics. Springer:Berlin, 2005. 89. Sehlhorst HG, J¨ anicke R, D¨ uster A, Rank E, Steeb H, Diebels S. Numerical investigations of foam-like materials by nested high-order finite element methods. Computational Mechanics 2009; 45 45–59.

c 2000 John Wiley & Sons, Ltd. Copyright Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2000; 00:1–6