- Email: [email protected]
Applying the hp–d version of the FEM to locally enhance dimensionally reduced models
Comput. Methods Appl. Mech. Engrg. 196 (2007) 3524–3533 www.elsevier.com/locate/cma
Applying the hp–d version of the FEM to locally enhance dimensionally reduced models A. Du¨ster *, A. Niggl, E. Rank Lehrstuhl fu¨r Bauinformatik, Fakulta¨t fu¨r Bauingenieur- und Vermessungswesen, Technische Universita¨t Mu¨nchen, Arcisstrasse 21, D-80290 Mu¨nchen, Germany Received 1 November 2005; accepted 30 October 2006 Available online 12 March 2007 Dedicated to Professor Dr. Ivo Babusˇka, Austin, USA, on the occasion of his 80th birthday.
Abstract A modification of the hp–d method is presented, which allows to locally enhance a dimensionally reduced finite element approximation with an Ansatz based on a more accurate, i.e. dimensionally higher, model. The procedure is based on a hierarchical domain decomposition which is applied to couple discretizations of different types of mechanical models, like beams, plates and solids. Thereby, the error of the underlying mathematical model can be reduced in critical regions by introducing locally a physically more sophisticated model. Two examples are given to demonstrate the basic properties of the proposed method. In the first example the solution of a one-dimensional beam model is improved locally by a two-dimensional finite element computation. In the second example the local enhancement of a Reissner–Mindlin plate by a fully three-dimensional continuum approach is considered, yielding a significant reduction of the model error. 2007 Elsevier B.V. All rights reserved. Keywords: hp–d method; Model adaptivity
1. Introduction In recent years, not only the control of the discretization error of finite element computations but also the model error, i.e. the error of the underlying mathematical model has been focussed on. To control the model error Vogelius and Babusˇka [1–3] have proposed to use a sequence of hierarchic models. Szabo´ and Sahrmann [4] have presented hierarchic plate and shell models which are based on the p-version of the FEM [5,6]. In [4] it was demonstrated by several numerical examples that this approach is suitable for a large variety of practical problems which may concurrently contain thin and thick plates, shells, stiffeners, and regions where a three-dimensional representation is *
Corresponding author. Tel.: +49 89 289 25060; fax: +49 89 289 25051. E-mail addresses: [email protected] (A. Du¨ster), [email protected] (A. Niggl), [email protected] (E. Rank). 0045-7825/$ - see front matter 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2006.10.018
required. In [6–17] this formulation has been further developed and applied to (adaptive) linear as well as nonlinear computations arising in hyperelasticity, elastoplasticity, thermo-elastic problems, spring-back analysis, elastodynamics and fluid–structure interaction. Schwab [18] has presented a computable model error estimator for hierarchic plates. The proposed error estimator is based on the assumption that the resulting twodimensional plate models are solved exactly. In order to estimate the model error, the residual tractions on the faces of the plate are computed. It is demonstrated that the residual based estimators are guaranteed upper estimators which follow the model error accurately as the model order increases. Oden and Cho [19] propose a four-step hpq-adaptive finite element method, in which q and h; p control the hierarchical models and the finite element meshes, respectively. The approach is based on an a priori hierarchical model
A. Du¨ster et al. / Comput. Methods Appl. Mech. Engrg. 196 (2007) 3524–3533
error estimate, an a priori error estimate of hp-finite element methods, and an a posteriori error estimator. By considering plate- and shell-like structures, it was demonstrated that the hpq-adaptive scheme is able to control the model error as well as the approximation error. In [20,21] Stein and Ohnimus presented an error-controlled model and discretization adaptivity within a hierarchical sequence of mathematical models, based on beam, plate and shell-theories. The adaptivity is performed by expanding the model with an anisotropic h- and p-adaptive procedure. The expansion method for enhancing the model is shown to be more efficient and evident for engineers than the reduction method. For an overview of error-controlled discretization and model adaptivity the reader is referred to the recent work of Stein and Ru¨ter [22]. In this paper we present a modification of the hp–d method which was presented in the early 1990s to consistently couple local and global structural computations [23–25]. This approach was subsequently developed further and applied to problems of plane elasticity in [26,27]. It has similarities to the s-method of Fish [28–32], can yet be implemented very efficiently as an iteration of separate finite element computations on different scales. Here, we extend the idea of the hp–d method to couple structural models of different dimensions. Using this procedure, it is possible to perform a fast, dimensionally reduced analysis of structures and to improve the approximation in regions where significant deviations from the exact three-dimensional solution occur. After a short introduction to the basic concept of the hp–d method we explain how this method can be applied to enhance dimensionally reduced models with a locally higher dimensional analysis. Therefore, first we consider a one-dimensional beam model which is locally enhanced by a two-dimensional finite element computation. Next, the concept is extended to locally improve Reissner–Mindlin plate computations by a strict three-dimensional finite element approach based on anisotropic high order finite elements [6,7,11].
3525
2. The hp–d method In Fig. 1 the basic idea of the hp–d method, which was first proposed by Rank [23–25] is explained for a onedimensional model problem. Suppose that a finite element approximation on a domain X consists of two parts: ub is defined as a standard h or p-approximation 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 fine overlay mesh is superimposed to the base mesh. On this overlay mesh a second part of the approximation, uo, is defined and the overall approximation uhp–d is composed of the hierarchical sum uhp–d ¼ ub þ uo . The C0-continuity of uhp–d is guaranteed by imposing homogeneous Dirichlet boundary conditions on the overlay approximation. The same procedure can be applied to problems in two or three dimensions. The two-dimensional situation is sketched in Fig. 2. The domain of computation Xb with boundary Cb is meshed with triangular or quadrilateral elements. Curved boundaries can be taken care of by applying the blending function method [5,6] on the base as well as on the overlay mesh. A subset of Xb is defined as the domain Xo with boundary Co. Please note that Xo does not need to be a connected subset of Xb. To be more precise, we now set up the basic formulation of the hp–d method. Suppose that Bðu; vÞ ¼ FðvÞ 8v 2 V
ð1Þ
describes the weak form of equilibrium of a two-dimensional linear elasticity problem. The bilinear form is then defined as Z Z T Bðu; vÞ ¼ eT ðvÞCeðuÞt dx dy ¼ ½Lv C½Lut dx dy; ð2Þ X
where
Fig. 1. Basic idea of the hp–d method for one-dimensional problems.
X
3526
A. Du¨ster et al. / Comput. Methods Appl. Mech. Engrg. 196 (2007) 3524–3533
ub ¼ Nb Ub ; uo ¼ No Uo ;
Fig. 2. The hp–d method: domain decomposition and definition of boundaries.
oux ouy oux ouy ; eyy ¼ ;c ¼ þ eðuÞ ¼ exx ¼ ox oy xy oy ox
T ð3Þ
corresponds to the strain in Voigt notation, L is the standard strain–displacement operator, C is the two-dimensional elasticity matrix and t corresponds to the thickness of the structure. The load functional Z Z T FðvÞ ¼ v ft dx dy þ vTtt dC ð4Þ X
CN
considers volume loads f ¼ ½fx ; fy and tractions, where t ¼ ½tx ; ty T denotes the traction vector at the Neumann boundary CN Cb . In order to find a finite element approximation to the weak formulation (1) by applying the hp–d method, we again define the hierarchical sum CD Cb : uhpd ¼ u;
Co n ðCb \ Co Þ : uo ¼ 0;
ð5Þ
where ub is the approximation corresponding to the discretization of the domain Xb and u0 corresponds to the domain Xo. The essential system boundary conditions are imposed on the Dirichlet boundary CD Cb , whereas uo ¼ 0 on Co n ðCb \ Co Þ guarantees the overall C0-continuity of uhpd . Note that homogeneous Dirichlet boundary conditions need to be defined only on Co n ðCb \ Co Þ. If natural boundary conditions are to be applied on CN Cb then these conditions need to be applied also to the boundary Co \ CN . By inserting (5) in (1) we obtain Bðub þ uo ; vb Þ ¼ Fðvb Þ 8vb 2 V b ; Bðub þ uo ; vo Þ ¼ Fðvo Þ
ð6Þ
8vo 2 V o :
Considering linear problems, for example, linear elasticity, we can split the bilinear form Bð; Þ and bring the mixed terms to the right-hand side. The resulting equation system can then be solved by applying a block Gauss–Seidel iteration: ðiþ1Þ
Bðub
; vb Þ ¼ Fðvb Þ BðuðiÞ o ; vb Þ
Bðuðiþ1Þ ; vo Þ o
¼ Fðvo Þ
ðiþ1Þ Bðub ; vo Þ
where Nb, No are the matrices of the shape functions of the base and the overlay mesh, respectively. Care has to be taken that no linear dependencies in the shape function spaces are present in the discretization in order to avoid the resulting stiffness matrix to become singular. A detailed discussion on how to construct corresponding shape function spaces and a strategy to accelerate the block Gauss– Seidel-iteration can be found in [26,27]. Corresponding to a standard Bubnov–Galerkin scheme the test functions vb and vo are discretized using the same shape functions Nb and No. The resulting linear equation system corresponding to (6) Kbb Kbo Ub Fb ¼ ð10Þ T Kbo Koo Uo Fo reflects the hierarchical nature of this formulation. To solve this coupled equation system we apply a block Gauss– Seidel-iteration ðiþ1Þ
Kbb Ub
¼ Fb Kbo UðiÞ o ;
ð11Þ
ðiþ1Þ
T
uhp–d ¼ ub þ uo ;
ð8Þ ð9Þ
8vb 2 V b ; 8vo 2 V o :
ð7Þ
The displacement fields ub and uo are now discretized in a standard way by an h or a p-approximation
Koo Uoðiþ1Þ ¼ Fo KTbo Ub
being the discretized version of (7). The advantage of this approach can be seen by examining the coupling terms Kbo and KTbo together with the corresponding iterates UoðiÞ ðiþ1Þ and Ub , respectively. Considering Kbo UðiÞ o we find that this expression can be interpreted as a load vector Z Z ðiÞ T ðiÞ Fbo ¼ Kbo UðiÞ ¼ B CB t dXU ¼ BTb CðiÞ ð12Þ o o b o o t dX Xo
Xo
corresponding to negative pre-strains resulting from the displacement UðiÞ o . Bb ¼ LNb and Bo ¼ LNo are the B– matrices in the standard notation as the strain operator applied to the matrices of the shape functions of the base and the overlay mesh, respectively. Analogously, we find that Z Z ðiþ1Þ ðiþ1Þ ðiþ1Þ ðiþ1Þ Fob ¼ KTbo Ub ¼ BTo CBb t dXUb ¼ BTo Ceb t dX: Xo
Xo
ð13Þ
The advantage of the interpretation of (12) and (13) is that it is not necessary to compute the coupling matrices explicitly. It is therefore possible to implement the hp–d approximation using a given finite element program as a black box, simply exchanging the corresponding pre-strains from the local to the global mesh and vice versa. The only necessary additional program module is an interpolation of strains from one mesh to the other and the possibility of accurately incorporating pre-strains as distributed loads resulting in a pseudo-load-vector, see also Remark 1. Further details of the hp–d method can be found in [24,23,25,33,26,34,11,27]. The hp–d method can be extended by taking more than just one overlay mesh into account. Then the finite element approximation reads
A. Du¨ster et al. / Comput. Methods Appl. Mech. Engrg. 196 (2007) 3524–3533
uhp–d ¼ ub þ
n X
CD Cb : uhp–d ¼ u;
uo j ;
j¼1
Coj n ðCb \ Coj Þ : uoj ¼ 0;
j ¼ 1; . . . ; n:
ð14Þ
Inserting the hp–d approximation into the bilinear form, we obtain Bðub ; vb Þ ¼ Fðvb Þ
n X
Bðuoj ; vb Þ
8vb 2 V b ;
j¼1
Bðuo1 ; vo1 Þ ¼ Fðvo1 Þ Bðub ; vo1 Þ Bðuo2 ; vo2 Þ ¼ Fðvo2 Þ Bðub ; vo2 Þ
n X j¼1;j6¼1 n X
Bðuoj ; vo1 Þ
8vo1 2 V o1 ;
Bðuoj ; vo2 Þ
8vo2 2 V o2 ;
Bðuon ; von Þ
8von 2 V on :
j¼1;j6¼2
.. . Bðuon ; von Þ ¼ Fðvon Þ Bðub ; von Þ
n X
3527
where the pseudo-load-vector is computed without integrating across discontinuities in the strain field. In [26] the importance of an accurate integration based on a composed integration rule has been demonstrated. In the numerical examples presented below, only coupling between base and overlay meshes is considered, i.e. no coupling between overlay meshes is taken into account. Furthermore, we apply just one block Gauss–Seidel iteration step. Therefore, the computation of the pseudo-loadvector is strongly simplified, since no special quadrature technique, such as a composed integration method is necessary. The only additional program module needed is to evaluate the strains of the base mesh at the quadrature points of the overlay mesh and the possibility to compute a pseudo-load-vector on the overlay mesh due to the (pre-)strains obtained from the base mesh with a standard Gaussian integration.
j¼1;j6¼n
ð15Þ This approach can be used to ‘‘plaster’’ incompatible finite element meshes of n local sub-models. In Fig. 3 the situation with n = 2 overlay meshes is depicted. Note that the two overlay meshes might (partially) overlap. C0-continuity of uhp–d is enforced by applying homogeneous Dirichlet boundary conditions uoj ¼ 0 on Coj n ðCb \ Coj Þ. Remark 1. Applying a block Gauss–Seidel iteration allows to avoid the explicit computation of the coupling matrices. To this end, the iterates are interpreted as pseudo-loadvectors from pre-strains, as already mentioned above. This involves the extension of the finite element code to be able to interpolate strains from one mesh to the other and to compute load vectors due to pre-strains. In order to accurately compute these pseudo-load-vectors, it is necessary to perform a composed integration, accounting for the different non-matching meshes involved in the hp–d discretization. This is due to the fact that the displacement field on each of the different meshes is approximated with C0-continuous functions; therefore the strains are discontinuous across the element boundaries. However, a Gauss rule cannot integrate discontinuous functions accurately. Therefore, a composed integration rule should be applied,
3. Dimensional adaptivity based on the hp–d FEM The core of the algorithm, as shown in the last section, is the splitting of the approximate solution into a local and a global part. It is not necessary for these two parts to be computed using the same method. One could, for instance, use this idea for coupling a finite element to a boundary element approximation. Here, we apply the hp–d method to combine discretizations of different mechanical models. The aim is to improve a finite element computation of a dimensionally reduced model locally by enhancing it with a higher dimensional model. This allows to decrease the model error in those regions, where the dimensionally reduced model deviates significantly from the more precise higher dimensional model. Thereby a sequence of hierarchical models can be defined, starting from a very simple Euler–Bernoulli beam theory, for example, while the highest model corresponds to the strict three-dimensional set of equations of linear elasticity. Local nonlinear material behaviour can also be treated with the hp–d method, when an adequate domain decomposition is applied. In [35] this approach has been used to compute a strip footing near a slope under vertical loading with an elasto-plastic material model to describe the local nonlinear behaviour of the soil. In the following we will apply the hp–d method in order to couple mechanical models of different dimension. Two examples are given to demonstrate the overall approach of the proposed dimensional adaptive method. 3.1. Coupling 1D and 2D elasticity problems
Fig. 3. The hp–d method: domain decomposition and definition of boundaries in the case of two overlay meshes.
To explain the basic idea, in the first example the global approximation ub is computed from a one-dimensional beam model. In the case of uniaxial bending, the solution is given by the deflection w(x) and the rotation b(x). In regions where a critical model error is to be expected, a two-dimensional mesh is defined locally. To match the spatial dimensions of the local and the global solution parts
3528
A. Du¨ster et al. / Comput. Methods Appl. Mech. Engrg. 196 (2007) 3524–3533
the 1D solution ub is formally extended via the kinematic model assumptions to the two-dimensional case ux ðx; yÞ bðxÞy ub ðx; yÞ ¼ ¼ : ð16Þ uy ðxÞ wðxÞ From this extended displacement field the extended strains T oux ouy oux ouy eðub Þ ¼ exx ¼ ; eyy ¼ ; cxy ¼ þ ð17Þ ox oy oy ox are computed, acting as negative loads to the second part of Eq. (7) Z Bðuo1 ; vo1 Þ ¼ Fðvo1 Þ eT ðvo1 ÞCeð ub Þt dX: ð18Þ Xo1
The stresses T ðeÞ ¼ rxx ; ryy ; rxy ¼ Ceð ub Þ r
Fig. 4. One-dimensional model of the beam with two holes.
ð19Þ
are computed from the strains of the extended beam solution, where 2 3 1 m 0 E 6m 1 0 7 ð20Þ C¼ 4 5 ð1 m2 Þ ð1mÞ 0 0 2 denotes the elasticity matrix for the case of plane stress, with E being Young’s modulus and m Poisson’s ratio. Applying either the Timoshenko or the Euler–Bernoulli beam theory will yield eyy ¼ 0, since in both models it is assumed that the deflection uy ¼ wðxÞ is independent of y. On the other hand, cxy will be non zero only in the case of the Timoshenko beam theory, which accounts for shear deformations. Discretizing the weak form (18), defined on the domain Xo1 with a finite element Ansatz uo1 and solving the resulting equation system, we obtain the approximation uo1 on the overlay mesh. In order to compute the pseudo-loadvector in (18) the strains eð ub Þ need to be interpolated to the integration points of the two-dimensional elements composing the overlay mesh. Having obtained the approximation on the overlay mesh, the hp–d approximation is composed of uhp–d ¼ ub þ uo1 . Generally, the overlay solution uo1 must be projected onto the dimensionally reduced model in order to compute the next iteration step, see Eq. (7). However, the numerical tests presented in this paper show that this iteration is not necessary. Only one iteration step is needed to obtain an accurate extension to the dimensionally reduced global solution. Of course, we have to bear in mind that applying only one iteration step corresponds to an approximate solution of the resulting equation system arising from the discretization of (7). There might be problems, where more than just one iteration is needed. To obtain a fast convergence of the block Gauss–Seidel iteration one could then take advantage of the acceleration technique proposed in [27]. The basic idea of this technique is to improve the initial solution of the iterates by the hp–d method itself on a very coarse discretization. The initial solution is then
computed by an iteration with the number of degrees of freedom being very small. Therefore almost no computational effort is needed, even if this iteration needs many steps to converge, yielding improved initial iterates which can be used for a refined hp–d discretization. In the following example we consider a linear elastic beam with two holes. The beam of length l = 10 m, width b ¼ 0:3 m and height h ¼ 1:0 m is clamped at both ends and subjected to a uniform load p = 1 MN/m. Fig. 4 shows the one-dimensional model of the structure. Assuming the Euler-Bernoulli theory of beams, an analytical solution for the deflection w(x) and the rotation b(x) can be given by wðxÞ ¼
1 1 1 pl2 x2 plx3 þ px4 ; 24EI 12EI 24EI
bðxÞ ¼
ow ; ox ð21Þ
where E = 2.1 · 105 MN/m2, and I ¼ bh3 =12 is the moment of inertia. Note that this one-dimensional solution neglects the local structure, i.e. the holes in the beam completely. The procedure to improve the 1D solution in the vicinity of the left hole is sketched in Fig. 5. The first contour plot shows the undisturbed stress rxx of the 1D beam solution (21) extended to the two-dimensional case (19). To compute the extended stresses (19) Poisson’s ratio in (20) has been set to m ¼ 0:3. Based on this solution, uo1 is computed on the overlay mesh according to Eq. (18). Finally, the composed solution rxx ðuhpd Þ ¼ rxx ðub Þ þ rxx ðuo1 Þ is shown at the bottom of Fig. 5. In Figs. 6 and 7 the hp–d results are compared with a two-dimensional reference solution, plotted along two different cut-lines. The location of the cut-lines is depicted in Fig. 5. The cut-line A–A corresponds to a fixed value of x = 4 m, intersecting the left circular hole. Cut-line B–B is located at x ¼ 4:5 m. The reference solution has been obtained by discretizing the whole beam as a 2D problem on a very fine mesh with high order elements. At the leftand the right-hand side of the two-dimensional mesh, homogeneous Dirichlet boundary conditions have been applied. From Figs. 6 and 7 it is evident that the hp–d
A. Du¨ster et al. / Comput. Methods Appl. Mech. Engrg. 196 (2007) 3524–3533
3529
Fig. 5. Improvement of dimensionally reduced approximation in the vicinity of the left hole.
The same approach can be used to ‘‘plaster’’ incompatible finite element meshes of n = 2 local sub-models. Each single overlay problem (22 and 23) can be computed thereby completely independently
Fig. 6. Stress component rxx along cut-line A–A.
Bðuo1 ; vo1 Þ ¼ Fðvo1 Þ Bðub ; vo1 Þ
8vo1 2 V o1 ;
ð22Þ
Bðuo2 ; vo2 Þ ¼ Fðvo2 Þ Bðub ; vo2 Þ
8vo2 2 V o2
ð23Þ
by ignoring the coupling terms Bðuo2 ; vo1 Þ and Bðuo1 ; vo2 Þ between the two overlay discretizations, see Eq. (15). Numerical investigations have demonstrated that ignoring the coupling terms still allows for a strong improvement of the beam solution. Fig. 8 illustrates this technique for the previously presented example. The solution at the hole on the right-hand side of the beam is computed on a different mesh. The two overlay meshes need not necessarily match on their boundaries and they can overlap. Again, a very good agreement of the hp–d approximation with the reference solution can be observed, even with only one block Gauss–Seidel iteration step. Remark 2. Numerical experiments have shown that the accuracy of the two-dimensional hp–d approximation can be further increased when improving the quality of the beam solution. If one computes the shear stresses from equilibrium conditions rather than from constitutive equations, the hp–d approximation coincides even better with the reference solution. One could choose also a Timoshenko beam theory instead, which includes already the shear stresses. However, we have decided to choose the simplest possible beam model as a starting point—which is of course more challenging—in order to demonstrate the capabilities of the proposed hp–d method.
Fig. 7. Stress component rxx along cut-line B–B.
3.2. Coupling 2D and 3D elasticity problems method significantly improves the accuracy of the beam solution. The model error is strongly reduced when compared to the two-dimensional reference solution.
The second numerical example to be considered is a three-dimensional thick-walled plate which is clamped on four sides (indicated by the shaded areas) and loaded with
3530
A. Du¨ster et al. / Comput. Methods Appl. Mech. Engrg. 196 (2007) 3524–3533
Fig. 8. Improvement of the 1D beam approximation by two independent overlay solutions, stress component rxx is plotted.
a pressure tz ¼ q acting only on a small part of the upper surface, Fig. 9. In order to solve this problem with the hp–d method, we start from the weak form of equilibrium based on the Reissner–Mindlin plate theory Bðu; vÞ ¼ FðvÞ 8v 2 V ;
ð24Þ
where the bilinear form and the linear functional define the internal and external virtual work, respectively. Vector T u ¼ ½bx ðx; yÞ; by ðx; yÞ; wðx; yÞ denotes the unknowns, where bx ðx; yÞ and by ðx; yÞ are the normal rotations in x–z and y–z planes, and wðx; yÞ corresponds to the deflection of the plate. The bilinear form reads
Bðu; vÞ ¼
t3 12
Z
eTB ðvÞCB eB ðuÞdX þ t
X
Z
eTS ðvÞCS eS ðuÞdX:
X
ð25Þ Based on the kinematic assumptions of the Reissner–Mindlin plate theory, the three-dimensional Cartesian displacement components, ux ðx; y; zÞ; uy ðx; y; zÞ and uz ðx; yÞ, are computed from the rotations and the deflection as u ¼ ½ux ¼ bx ðx; yÞz; uy ¼ by ðx; yÞz; uz ¼ wðx; yÞT : The strain term in (25),
Fig. 9. Three-dimensional plate with boundary conditions.
ð26Þ
A. Du¨ster et al. / Comput. Methods Appl. Mech. Engrg. 196 (2007) 3524–3533 T
are small. The constitutive equation describing the relation between strains and stresses is given by an isotropic linear elastic law
e ¼ ½exx ; eyy ; cxy ; cyz ; cxz T obx oby obx oby ow ow ; bx þ z; z; þ ¼ z; by þ oy ox ox oy oy ox
T
T
¼ ½eTB z; eTS
ð27Þ
denotes the (small) strains obtained by computing the corresponding derivatives of the displacement field. The stress–strain relation is given by the constitutive equations CB 0 r ¼ ½rxx ; ryy ; rxy ; ryz ; rxz T ¼ e; ð28Þ 0 CS where the matrices 2 1 m E 6m 1 CB ¼ 4 ð1 m2 Þ 0 0
0
3
0 7 5;
ð1mÞ 2
0 1
ð29Þ define the assumed isotropic linear elastic behaviour of the Reissner–Mindlin plate. The warping coefficient k ¼ 56 accounts for non-uniform shear distributions. The load functional Z qv dX ð30Þ FðvÞ ¼ X
takes into account the load q acting on the surface of the plate. The idea of the hp–d method with its application to dimensional adaptivity is to enhance the dimensionally reduced model, i.e. the Reissner–Mindlin plate theory, in those regions, where a significant deviation from the three-dimensional equations of linear elasticity are expected. We therefore apply the weak form of equilibrium of linear elasticity in 3D Bðu; vÞ ¼ FðvÞ
8v 2 V ;
ð31Þ
with Bðu; vÞ ¼ Z FðvÞ ¼
Z
r ¼ ½rxx ; ryy ; rzz ; rxy ; ryz ; rxz ¼ Ce; where 2 j þ 43 l 6 j 2 l 6 3 6 6 j 23 l 6 C¼6 0 6 6 4 0 0
1 Ek CS ¼ 2ð1 þ mÞ 0
eT ðvÞCeðuÞdX;
ð32Þ
X
vTt dC
ð36Þ
j 23 l j 23 l 0 j þ 43 l j 23 l 0 j 23 l j þ 43 l 0 0 0 l 0 0
0 0
0 0
3
0
0
0
07 7 7 07 7 07 7 7 05
0 0 l 0
ð37Þ
l
is the elasticity matrix, j ¼ E=ð3ð1 2mÞÞ, the bulk modulus, l ¼ E=ð2ð1 þ mÞÞ, the shear modulus. In a first step, an analysis based on the Reissner–Mindlin plate theory is performed by applying the p-version on a base mesh consisting of 108 quadrilaterals in combination q with the tensor product space Sp;p ps ðXst Þ [6,7,11] and a polynomial degree of p = 8, Fig. 10. The mesh is refined towards the boundary and at reentrant corners in order to resolve boundary layers and singularities, respectively. Following the procedure outlined above, the two-dimensional plate solution is then extended via the Reissner– Mindlin assumptions to a three-dimensional approximation yielding the displacements (26) as ub ¼ ½ux ¼ bx ðx; yÞz; uy ¼ by ðx; yÞz; uz ¼ wðx; yÞT
ð38Þ
and the corresponding strains ðub Þ ¼ ½exx ; eyy ; ezz ; cxy ; cyz ; cxz T T obx oby obx oby ow ow ; bx þ z; z; 0; þ ¼ ; z; by þ oy ox ox oy oy ox ð39Þ with ezz ¼ 0. In the vicinity of the area where the load is applied, an overlay mesh with 17 hexahedral elements is superimposed on the base mesh. An anisotropic Ansatz space Sp;p;q ðXhst Þ [6,7,11] in conjunction with the polynomial degree template
ð33Þ
CN
ð40Þ
as a basis for a finite element discretization. The function describing the Cartesian displacement components is denoted by u ¼ ½ux ðx; y; zÞ; uy ðx; y; zÞ; uz ðx; y; zÞT :
3531
ð34Þ
It is assumed that strains e ¼ ½exx ; eyy ; ezz ; cxy ; cyz ; cxz T T oux ouy ouz oux ouy ouy ouz oux ouz ; ; ; þ ; þ ; þ ¼ ox oy oz oy ox oz oy oz ox ð35Þ
is chosen for the three-dimensional computation, where the polynomial degree for the displacement components ux ; uy in in-plane is set to p = 7 and in transversal direction to q = 4. The deflection uz is approximated with an anisotropic Ansatz, where the polynomial degree in in-plane is set to p = 8 and in transversal direction to q = 5. The polynomial degrees in (40) have been chosen such that the discretization error of the anisotropic three-dimensional approximation is small. Since a hierarchic Ansatz space is applied, the quality of the approximation can be readily
3532
A. Du¨ster et al. / Comput. Methods Appl. Mech. Engrg. 196 (2007) 3524–3533
Fig. 10. 2D base and 3D overlay mesh for a local three-dimensional discretization.
controlled by increasing the polynomial degrees. An adaptive procedure to perform this adjustment in an automatic manner has been suggested, for example, in [15,36]. Applying the anisotropic discretization, the weak form Z Bðuo ; vo Þ ¼ Fðvo Þ eT ðvo ÞCeð ub ÞdX ð41Þ Xo
including the pseudo-load-vector resulting from the Reissner–Mindlin approximation is solved. The pseudo-loadvector is computed elementwise by interpolating the strains eð ub Þ from the extended Reissner–Mindlin model to the integration points of the hexahedral elements composing the overlay mesh. Having solved (41), the local–global three-dimensional hp–d solution is then defined as the hierarchical sum uhp–d ¼ ub þ uo of the extended global plate approximation ub and the local three-dimensional approximation uo obtained on the overlay mesh. As a result of the computations, consider the von Mises stress along the cut-lines C–C and D–D plotted in Figs. 11 and 12. The reference solution is based on a three-dimensional overkill approximation with a numerical error which is negligibly small. From Figs. 11 and 12 it is evident that the Reissner–Mindlin approximation coincides very well with the reference solution for a wide range of the domain. However, in the vicinity of the area where the load is
Fig. 12. von Mises stress along cut-line D–D.
applied, the two-dimensional approximation differs significantly from the reference solution. Comparing now the hp–d approximation with the reference solution, one observes that the quality of the finite element solution is enormously improved. Deviations from the reference solution are almost invisible. It is also important to notice that only one block Gauss–Seidel iteration step was performed in this example. Compared to a three-dimensional approximation of the whole problem, the additional numerical effort of the hp–d method is very small. 4. Conclusions
Fig. 11. von Mises stress along cut-line C–C.
The outlined method of combining dimensionally reduced approximations with higher dimensional models can be an attractive method leading to strongly improved solutions using only a small additional numerical effort. Further improvements in terms of efficiency and reliability could be achieved by applying an adaptive scheme in order to localize the domain where the dimensionally reduced approximation is to be extended by a higher dimensional Ansatz automatically. An error estimator could be used to quantify the model error. Areas with a high model error would therefore indicate where to enrich the dimensionally reduced approximation by means of a higher dimensional
A. Du¨ster et al. / Comput. Methods Appl. Mech. Engrg. 196 (2007) 3524–3533
Ansatz. Future work will focus on the development of such an adaptive scheme. It should finally be noted that, although the p-version of FEM was used for the examples presented in this paper, the method is not restricted to higher order finite element approximations. Low order FEM as available in any commercial finite element program can be used likewise. The only programming effort for implementing this procedure is the interpolation of strains obtained from the dimensionally reduced, global model to the extended local approximation and the computation of load vectors due to pre-strains. References [1] M. Vogelius, I. Babusˇka, On a dimensional reduction method. I. The optimal selection of basis function, Math. Comput. 37 (1981) 31–46. [2] M. Vogelius, I. Babusˇka, On a dimensional reduction method. II. Some approximation-theoretic results, Math. Comput. 37 (1981) 47– 48. [3] M. Vogelius, I. Babusˇka, On a dimensional reduction method. III. A posteriori error estimation and an adaptive approach, Math. Comput. 37 (1981) 361–384. [4] B. Szabo´, G. Sahrmann, Hierarchic plate and shell models based on p-extension, Int. J. Numer. Methods Engrg. 26 (1988) 1855–1881. [5] B. Szabo´, I. Babusˇka, Finite Element Analysis, John Wiley & Sons, 1991. [6] B. Szabo´, A. Du¨ster, E. Rank, The p-version of the finite element method, in: E. Stein, R. de Borst, T. Hughes (Eds.), Encyclopedia of Computational Mechanics, vol. 1, John Wiley & Sons, 2004, pp. 119– 139 (Chapter 5). [7] A. Du¨ster, H. Bro¨ker, E. Rank, The p-version of the finite element method for three-dimensional curved thin walled structures, Int. J. Numer. Methods Engrg. 52 (2001) 673–703. [8] A. Du¨ster, E. Rank, A p-version finite element approach for two- and three-dimensional problems of the J2 flow theory with non-linear isotropic hardening, Int. J. Numer. Methods Engrg. 53 (2002) 49–63. [9] A. Du¨ster, A. Niggl, V. Nu¨bel, E. Rank, A numerical investigation of high-order finite elements for problems of elasto-plasticity, J. Sci. Comput. 17 (2002) 429–437. [10] A. Du¨ster, A. Niggl, E. Rank, Thermo-elastic computations of geometrically non-linear three-dimensional thin-walled continua based on high order finite elements, in: H. Mang, F. Rammerstorfer, J. Eberhardsteiner (Eds.), Proceedings of the Fifth World Congress on Computational Mechanics, Vienna, Austria, 2002. [11] A. Du¨ster, High order finite elements for three-dimensional, thinwalled nonlinear continua, Ph.D. thesis, Lehrstuhl fu¨r Bauinformatik, Fakulta¨t fu¨r Bauingenieur- und Vermessungswesen, Technische Universita¨t Mu¨nchen. 2001. [12] A. Du¨ster, S. Hartmann, E. Rank, p-fem applied to finite isotropic hyperelastic bodies, Comput. Methods Appl. Mech. Engrg. 192 (2003) 5147–5166. [13] A. Muthler, A. Du¨ster, W. Volk, M. Wagner, E. Rank, High order thin-walled solid finite elements applied to elastic spring-back computations, Comput. Methods Appl. Mech. Engrg. 195 (2006) 5377–5389. [14] E. Rank, A. Du¨ster, V. Nu¨bel, K. Preusch, O. Bruhns, High order finite elements for shells, Comput. Methods Appl. Mech. Engrg. 194 (2005) 2494–2512. [15] A. Du¨ster, D. Scholz, E. Rank, pq-Adaptive solid finite elements for three-dimensional plates and shells, Comput. Methods Appl. Mech. Engrg., submitted for publication.
3533
[16] A. Du¨ster, D. Scholz, E. Rank, Fluid–structure interaction using high-order anisotropic solid elements, in: Proceedings of the 5th International Conference on Computation of Shell and Spatial Structures, Salzburg, Austria, 2005. [17] D. Scholz, S. Kollmannsberger, A. Du¨ster, E. Rank, Thin solids for fluid–structure interaction, in: H. Bungartz, M. Scha¨fer (Eds.), Fluid–Structure Interaction, Modelling Simulation and Optimisation, Lecture Notes in Computational Science and Engineering, vol. 53, Springer, 2006, pp. 294–335. [18] C. Schwab, A-posteriori modeling error estimation for hierarchic plate models, Numer. Math. 74 (1996) 221–259. [19] J. Oden, J. Cho, Adaptive hpq-finite element methods of hierarchical models for plate- and shell-like structures, Comput. Methods Appl. Mech. Engrg. 136 (1996) 317–345. [20] E. Stein, S. Ohnimus, Coupled model- and solution-adaptivity in the finite-element method, Comput. Methods Appl. Mech. Engrg. 150 (1997) 327–350. [21] E. Stein, S. Ohnimus, Anisotropic discretization- and model-error estimation in solid mechanics by local Neumann problems, Comput. Methods Appl. Mech. Engrg. 176 (1999) 363–385. [22] E. Stein, M. Ru¨ter, Finite element methods for elasticity with errorcontrolled discretization and model adaptivity, in: E. Stein, R. de Borst, T. Hughes (Eds.), Encyclopedia of Computational Mechanics, vol. 2, John Wiley & Sons, 2004, pp. 5–58 (Chapter 2). [23] E. Rank, Adaptive remeshing and h-p domain decomposition, Comput. Methods Appl. Mech. Engrg. 101 (1992) 299–313. [24] E. Rank, A combination of hp-version finite elements and a domain decomposition method, in: Proceedings of the First European Conference on Numerical Methods in Engineering, Bru¨ssel, 1992. [25] E. Rank, A zooming-technique using a hierarchical hp-version of the finite element method, in: J. Whiteman (Ed.), The Mathematics of Finite Elements and Applications—Highlights 1993, Elsevier, Uxbridge, 1993. [26] R. Krause, Multiscale Computations with a Combined h- and p-version of the Finite-element Method, Ph.D. Thesis, Fach Numerische Methoden und Informationsverarbeitung, Universita¨t Dortmund, 1996. [27] R. Krause, E. Rank, Multiscale computations with a combination of the h- and p-versions of the finite-element method, Comput. Methods Appl. Mech. Engrg. 192 (2003) 3959–3983. [28] J. Fish, The s-version of the finite element method, Comput. Struct. 43 (3) (1992) 539–547. [29] J. Fish, A. Nath, Adaptive and hierarchical modelling of fatigue crack propagation, Int. J. Numer. Methods Engrg. 36 (1993) 2825–2836. [30] J. Fish, S. Markolefas, Adaptive s-method for linear elastostatics, Comput. Methods Appl. Mech. Engrg. 104 (1993) 363–396. [31] J. Fish, S. Markolefas, Adaptive global–local refinement strategy based on the interior error estimates of the h-method, Int. J. Numer. Methods Engrg. 37 (1994) 827–838. [32] J. Fish, M. Pandheeradi, V. Belsky, An efficient multilevel solution scheme for large scale non-linear systems, Int. J. Numer. Methods Engrg. 38 (1995) 1597–1610. [33] E. Rank, R. Krause, A-posteriori-Fehlerscha¨tzung bei Multi-SkalenProblemen, ZAMM 75 (1995), SI. [34] E. Rank, R. Krause, A multiscale finite-element-method, Comput. Struct. 64 (1997) 139–144. [35] A. Du¨ster, E. Rank, G. Steinl, W. Wunderlich, A combination of an h- and p-version of the finite element method for elastic-plastic problems, in: Proceedings of ECCM ’99, European Conference on Computational Mechanics, Mu¨nchen, Germany, 1999. [36] D. Scholz, An anisotropic p-adaptive method for linear elastostatic and elastodynamic analysis of thin-walled and massive structures, Dissertation, Lehrstuhl fu¨r Bauinformatik, Fakulta¨t fu¨r Bauingenieur- und Vermessungswesen, Technische Universita¨t Mu¨nchen, 2006.
Applying the hp–d version of the FEM to locally enhance dimensionally reduced models A. Du¨ster *, A. Niggl, E. Rank Lehrstuhl fu¨r Bauinformatik, Fakulta¨t fu¨r Bauingenieur- und Vermessungswesen, Technische Universita¨t Mu¨nchen, Arcisstrasse 21, D-80290 Mu¨nchen, Germany Received 1 November 2005; accepted 30 October 2006 Available online 12 March 2007 Dedicated to Professor Dr. Ivo Babusˇka, Austin, USA, on the occasion of his 80th birthday.
Abstract A modification of the hp–d method is presented, which allows to locally enhance a dimensionally reduced finite element approximation with an Ansatz based on a more accurate, i.e. dimensionally higher, model. The procedure is based on a hierarchical domain decomposition which is applied to couple discretizations of different types of mechanical models, like beams, plates and solids. Thereby, the error of the underlying mathematical model can be reduced in critical regions by introducing locally a physically more sophisticated model. Two examples are given to demonstrate the basic properties of the proposed method. In the first example the solution of a one-dimensional beam model is improved locally by a two-dimensional finite element computation. In the second example the local enhancement of a Reissner–Mindlin plate by a fully three-dimensional continuum approach is considered, yielding a significant reduction of the model error. 2007 Elsevier B.V. All rights reserved. Keywords: hp–d method; Model adaptivity
1. Introduction In recent years, not only the control of the discretization error of finite element computations but also the model error, i.e. the error of the underlying mathematical model has been focussed on. To control the model error Vogelius and Babusˇka [1–3] have proposed to use a sequence of hierarchic models. Szabo´ and Sahrmann [4] have presented hierarchic plate and shell models which are based on the p-version of the FEM [5,6]. In [4] it was demonstrated by several numerical examples that this approach is suitable for a large variety of practical problems which may concurrently contain thin and thick plates, shells, stiffeners, and regions where a three-dimensional representation is *
Corresponding author. Tel.: +49 89 289 25060; fax: +49 89 289 25051. E-mail addresses: [email protected] (A. Du¨ster), [email protected] (A. Niggl), [email protected] (E. Rank). 0045-7825/$ - see front matter 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2006.10.018
required. In [6–17] this formulation has been further developed and applied to (adaptive) linear as well as nonlinear computations arising in hyperelasticity, elastoplasticity, thermo-elastic problems, spring-back analysis, elastodynamics and fluid–structure interaction. Schwab [18] has presented a computable model error estimator for hierarchic plates. The proposed error estimator is based on the assumption that the resulting twodimensional plate models are solved exactly. In order to estimate the model error, the residual tractions on the faces of the plate are computed. It is demonstrated that the residual based estimators are guaranteed upper estimators which follow the model error accurately as the model order increases. Oden and Cho [19] propose a four-step hpq-adaptive finite element method, in which q and h; p control the hierarchical models and the finite element meshes, respectively. The approach is based on an a priori hierarchical model
A. Du¨ster et al. / Comput. Methods Appl. Mech. Engrg. 196 (2007) 3524–3533
error estimate, an a priori error estimate of hp-finite element methods, and an a posteriori error estimator. By considering plate- and shell-like structures, it was demonstrated that the hpq-adaptive scheme is able to control the model error as well as the approximation error. In [20,21] Stein and Ohnimus presented an error-controlled model and discretization adaptivity within a hierarchical sequence of mathematical models, based on beam, plate and shell-theories. The adaptivity is performed by expanding the model with an anisotropic h- and p-adaptive procedure. The expansion method for enhancing the model is shown to be more efficient and evident for engineers than the reduction method. For an overview of error-controlled discretization and model adaptivity the reader is referred to the recent work of Stein and Ru¨ter [22]. In this paper we present a modification of the hp–d method which was presented in the early 1990s to consistently couple local and global structural computations [23–25]. This approach was subsequently developed further and applied to problems of plane elasticity in [26,27]. It has similarities to the s-method of Fish [28–32], can yet be implemented very efficiently as an iteration of separate finite element computations on different scales. Here, we extend the idea of the hp–d method to couple structural models of different dimensions. Using this procedure, it is possible to perform a fast, dimensionally reduced analysis of structures and to improve the approximation in regions where significant deviations from the exact three-dimensional solution occur. After a short introduction to the basic concept of the hp–d method we explain how this method can be applied to enhance dimensionally reduced models with a locally higher dimensional analysis. Therefore, first we consider a one-dimensional beam model which is locally enhanced by a two-dimensional finite element computation. Next, the concept is extended to locally improve Reissner–Mindlin plate computations by a strict three-dimensional finite element approach based on anisotropic high order finite elements [6,7,11].
3525
2. The hp–d method In Fig. 1 the basic idea of the hp–d method, which was first proposed by Rank [23–25] is explained for a onedimensional model problem. Suppose that a finite element approximation on a domain X consists of two parts: ub is defined as a standard h or p-approximation 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 fine overlay mesh is superimposed to the base mesh. On this overlay mesh a second part of the approximation, uo, is defined and the overall approximation uhp–d is composed of the hierarchical sum uhp–d ¼ ub þ uo . The C0-continuity of uhp–d is guaranteed by imposing homogeneous Dirichlet boundary conditions on the overlay approximation. The same procedure can be applied to problems in two or three dimensions. The two-dimensional situation is sketched in Fig. 2. The domain of computation Xb with boundary Cb is meshed with triangular or quadrilateral elements. Curved boundaries can be taken care of by applying the blending function method [5,6] on the base as well as on the overlay mesh. A subset of Xb is defined as the domain Xo with boundary Co. Please note that Xo does not need to be a connected subset of Xb. To be more precise, we now set up the basic formulation of the hp–d method. Suppose that Bðu; vÞ ¼ FðvÞ 8v 2 V
ð1Þ
describes the weak form of equilibrium of a two-dimensional linear elasticity problem. The bilinear form is then defined as Z Z T Bðu; vÞ ¼ eT ðvÞCeðuÞt dx dy ¼ ½Lv C½Lut dx dy; ð2Þ X
where
Fig. 1. Basic idea of the hp–d method for one-dimensional problems.
X
3526
A. Du¨ster et al. / Comput. Methods Appl. Mech. Engrg. 196 (2007) 3524–3533
ub ¼ Nb Ub ; uo ¼ No Uo ;
Fig. 2. The hp–d method: domain decomposition and definition of boundaries.
oux ouy oux ouy ; eyy ¼ ;c ¼ þ eðuÞ ¼ exx ¼ ox oy xy oy ox
T ð3Þ
corresponds to the strain in Voigt notation, L is the standard strain–displacement operator, C is the two-dimensional elasticity matrix and t corresponds to the thickness of the structure. The load functional Z Z T FðvÞ ¼ v ft dx dy þ vTtt dC ð4Þ X
CN
considers volume loads f ¼ ½fx ; fy and tractions, where t ¼ ½tx ; ty T denotes the traction vector at the Neumann boundary CN Cb . In order to find a finite element approximation to the weak formulation (1) by applying the hp–d method, we again define the hierarchical sum CD Cb : uhpd ¼ u;
Co n ðCb \ Co Þ : uo ¼ 0;
ð5Þ
where ub is the approximation corresponding to the discretization of the domain Xb and u0 corresponds to the domain Xo. The essential system boundary conditions are imposed on the Dirichlet boundary CD Cb , whereas uo ¼ 0 on Co n ðCb \ Co Þ guarantees the overall C0-continuity of uhpd . Note that homogeneous Dirichlet boundary conditions need to be defined only on Co n ðCb \ Co Þ. If natural boundary conditions are to be applied on CN Cb then these conditions need to be applied also to the boundary Co \ CN . By inserting (5) in (1) we obtain Bðub þ uo ; vb Þ ¼ Fðvb Þ 8vb 2 V b ; Bðub þ uo ; vo Þ ¼ Fðvo Þ
ð6Þ
8vo 2 V o :
Considering linear problems, for example, linear elasticity, we can split the bilinear form Bð; Þ and bring the mixed terms to the right-hand side. The resulting equation system can then be solved by applying a block Gauss–Seidel iteration: ðiþ1Þ
Bðub
; vb Þ ¼ Fðvb Þ BðuðiÞ o ; vb Þ
Bðuðiþ1Þ ; vo Þ o
¼ Fðvo Þ
ðiþ1Þ Bðub ; vo Þ
where Nb, No are the matrices of the shape functions of the base and the overlay mesh, respectively. Care has to be taken that no linear dependencies in the shape function spaces are present in the discretization in order to avoid the resulting stiffness matrix to become singular. A detailed discussion on how to construct corresponding shape function spaces and a strategy to accelerate the block Gauss– Seidel-iteration can be found in [26,27]. Corresponding to a standard Bubnov–Galerkin scheme the test functions vb and vo are discretized using the same shape functions Nb and No. The resulting linear equation system corresponding to (6) Kbb Kbo Ub Fb ¼ ð10Þ T Kbo Koo Uo Fo reflects the hierarchical nature of this formulation. To solve this coupled equation system we apply a block Gauss– Seidel-iteration ðiþ1Þ
Kbb Ub
¼ Fb Kbo UðiÞ o ;
ð11Þ
ðiþ1Þ
T
uhp–d ¼ ub þ uo ;
ð8Þ ð9Þ
8vb 2 V b ; 8vo 2 V o :
ð7Þ
The displacement fields ub and uo are now discretized in a standard way by an h or a p-approximation
Koo Uoðiþ1Þ ¼ Fo KTbo Ub
being the discretized version of (7). The advantage of this approach can be seen by examining the coupling terms Kbo and KTbo together with the corresponding iterates UoðiÞ ðiþ1Þ and Ub , respectively. Considering Kbo UðiÞ o we find that this expression can be interpreted as a load vector Z Z ðiÞ T ðiÞ Fbo ¼ Kbo UðiÞ ¼ B CB t dXU ¼ BTb CðiÞ ð12Þ o o b o o t dX Xo
Xo
corresponding to negative pre-strains resulting from the displacement UðiÞ o . Bb ¼ LNb and Bo ¼ LNo are the B– matrices in the standard notation as the strain operator applied to the matrices of the shape functions of the base and the overlay mesh, respectively. Analogously, we find that Z Z ðiþ1Þ ðiþ1Þ ðiþ1Þ ðiþ1Þ Fob ¼ KTbo Ub ¼ BTo CBb t dXUb ¼ BTo Ceb t dX: Xo
Xo
ð13Þ
The advantage of the interpretation of (12) and (13) is that it is not necessary to compute the coupling matrices explicitly. It is therefore possible to implement the hp–d approximation using a given finite element program as a black box, simply exchanging the corresponding pre-strains from the local to the global mesh and vice versa. The only necessary additional program module is an interpolation of strains from one mesh to the other and the possibility of accurately incorporating pre-strains as distributed loads resulting in a pseudo-load-vector, see also Remark 1. Further details of the hp–d method can be found in [24,23,25,33,26,34,11,27]. The hp–d method can be extended by taking more than just one overlay mesh into account. Then the finite element approximation reads
A. Du¨ster et al. / Comput. Methods Appl. Mech. Engrg. 196 (2007) 3524–3533
uhp–d ¼ ub þ
n X
CD Cb : uhp–d ¼ u;
uo j ;
j¼1
Coj n ðCb \ Coj Þ : uoj ¼ 0;
j ¼ 1; . . . ; n:
ð14Þ
Inserting the hp–d approximation into the bilinear form, we obtain Bðub ; vb Þ ¼ Fðvb Þ
n X
Bðuoj ; vb Þ
8vb 2 V b ;
j¼1
Bðuo1 ; vo1 Þ ¼ Fðvo1 Þ Bðub ; vo1 Þ Bðuo2 ; vo2 Þ ¼ Fðvo2 Þ Bðub ; vo2 Þ
n X j¼1;j6¼1 n X
Bðuoj ; vo1 Þ
8vo1 2 V o1 ;
Bðuoj ; vo2 Þ
8vo2 2 V o2 ;
Bðuon ; von Þ
8von 2 V on :
j¼1;j6¼2
.. . Bðuon ; von Þ ¼ Fðvon Þ Bðub ; von Þ
n X
3527
where the pseudo-load-vector is computed without integrating across discontinuities in the strain field. In [26] the importance of an accurate integration based on a composed integration rule has been demonstrated. In the numerical examples presented below, only coupling between base and overlay meshes is considered, i.e. no coupling between overlay meshes is taken into account. Furthermore, we apply just one block Gauss–Seidel iteration step. Therefore, the computation of the pseudo-loadvector is strongly simplified, since no special quadrature technique, such as a composed integration method is necessary. The only additional program module needed is to evaluate the strains of the base mesh at the quadrature points of the overlay mesh and the possibility to compute a pseudo-load-vector on the overlay mesh due to the (pre-)strains obtained from the base mesh with a standard Gaussian integration.
j¼1;j6¼n
ð15Þ This approach can be used to ‘‘plaster’’ incompatible finite element meshes of n local sub-models. In Fig. 3 the situation with n = 2 overlay meshes is depicted. Note that the two overlay meshes might (partially) overlap. C0-continuity of uhp–d is enforced by applying homogeneous Dirichlet boundary conditions uoj ¼ 0 on Coj n ðCb \ Coj Þ. Remark 1. Applying a block Gauss–Seidel iteration allows to avoid the explicit computation of the coupling matrices. To this end, the iterates are interpreted as pseudo-loadvectors from pre-strains, as already mentioned above. This involves the extension of the finite element code to be able to interpolate strains from one mesh to the other and to compute load vectors due to pre-strains. In order to accurately compute these pseudo-load-vectors, it is necessary to perform a composed integration, accounting for the different non-matching meshes involved in the hp–d discretization. This is due to the fact that the displacement field on each of the different meshes is approximated with C0-continuous functions; therefore the strains are discontinuous across the element boundaries. However, a Gauss rule cannot integrate discontinuous functions accurately. Therefore, a composed integration rule should be applied,
3. Dimensional adaptivity based on the hp–d FEM The core of the algorithm, as shown in the last section, is the splitting of the approximate solution into a local and a global part. It is not necessary for these two parts to be computed using the same method. One could, for instance, use this idea for coupling a finite element to a boundary element approximation. Here, we apply the hp–d method to combine discretizations of different mechanical models. The aim is to improve a finite element computation of a dimensionally reduced model locally by enhancing it with a higher dimensional model. This allows to decrease the model error in those regions, where the dimensionally reduced model deviates significantly from the more precise higher dimensional model. Thereby a sequence of hierarchical models can be defined, starting from a very simple Euler–Bernoulli beam theory, for example, while the highest model corresponds to the strict three-dimensional set of equations of linear elasticity. Local nonlinear material behaviour can also be treated with the hp–d method, when an adequate domain decomposition is applied. In [35] this approach has been used to compute a strip footing near a slope under vertical loading with an elasto-plastic material model to describe the local nonlinear behaviour of the soil. In the following we will apply the hp–d method in order to couple mechanical models of different dimension. Two examples are given to demonstrate the overall approach of the proposed dimensional adaptive method. 3.1. Coupling 1D and 2D elasticity problems
Fig. 3. The hp–d method: domain decomposition and definition of boundaries in the case of two overlay meshes.
To explain the basic idea, in the first example the global approximation ub is computed from a one-dimensional beam model. In the case of uniaxial bending, the solution is given by the deflection w(x) and the rotation b(x). In regions where a critical model error is to be expected, a two-dimensional mesh is defined locally. To match the spatial dimensions of the local and the global solution parts
3528
A. Du¨ster et al. / Comput. Methods Appl. Mech. Engrg. 196 (2007) 3524–3533
the 1D solution ub is formally extended via the kinematic model assumptions to the two-dimensional case ux ðx; yÞ bðxÞy ub ðx; yÞ ¼ ¼ : ð16Þ uy ðxÞ wðxÞ From this extended displacement field the extended strains T oux ouy oux ouy eðub Þ ¼ exx ¼ ; eyy ¼ ; cxy ¼ þ ð17Þ ox oy oy ox are computed, acting as negative loads to the second part of Eq. (7) Z Bðuo1 ; vo1 Þ ¼ Fðvo1 Þ eT ðvo1 ÞCeð ub Þt dX: ð18Þ Xo1
The stresses T ðeÞ ¼ rxx ; ryy ; rxy ¼ Ceð ub Þ r
Fig. 4. One-dimensional model of the beam with two holes.
ð19Þ
are computed from the strains of the extended beam solution, where 2 3 1 m 0 E 6m 1 0 7 ð20Þ C¼ 4 5 ð1 m2 Þ ð1mÞ 0 0 2 denotes the elasticity matrix for the case of plane stress, with E being Young’s modulus and m Poisson’s ratio. Applying either the Timoshenko or the Euler–Bernoulli beam theory will yield eyy ¼ 0, since in both models it is assumed that the deflection uy ¼ wðxÞ is independent of y. On the other hand, cxy will be non zero only in the case of the Timoshenko beam theory, which accounts for shear deformations. Discretizing the weak form (18), defined on the domain Xo1 with a finite element Ansatz uo1 and solving the resulting equation system, we obtain the approximation uo1 on the overlay mesh. In order to compute the pseudo-loadvector in (18) the strains eð ub Þ need to be interpolated to the integration points of the two-dimensional elements composing the overlay mesh. Having obtained the approximation on the overlay mesh, the hp–d approximation is composed of uhp–d ¼ ub þ uo1 . Generally, the overlay solution uo1 must be projected onto the dimensionally reduced model in order to compute the next iteration step, see Eq. (7). However, the numerical tests presented in this paper show that this iteration is not necessary. Only one iteration step is needed to obtain an accurate extension to the dimensionally reduced global solution. Of course, we have to bear in mind that applying only one iteration step corresponds to an approximate solution of the resulting equation system arising from the discretization of (7). There might be problems, where more than just one iteration is needed. To obtain a fast convergence of the block Gauss–Seidel iteration one could then take advantage of the acceleration technique proposed in [27]. The basic idea of this technique is to improve the initial solution of the iterates by the hp–d method itself on a very coarse discretization. The initial solution is then
computed by an iteration with the number of degrees of freedom being very small. Therefore almost no computational effort is needed, even if this iteration needs many steps to converge, yielding improved initial iterates which can be used for a refined hp–d discretization. In the following example we consider a linear elastic beam with two holes. The beam of length l = 10 m, width b ¼ 0:3 m and height h ¼ 1:0 m is clamped at both ends and subjected to a uniform load p = 1 MN/m. Fig. 4 shows the one-dimensional model of the structure. Assuming the Euler-Bernoulli theory of beams, an analytical solution for the deflection w(x) and the rotation b(x) can be given by wðxÞ ¼
1 1 1 pl2 x2 plx3 þ px4 ; 24EI 12EI 24EI
bðxÞ ¼
ow ; ox ð21Þ
where E = 2.1 · 105 MN/m2, and I ¼ bh3 =12 is the moment of inertia. Note that this one-dimensional solution neglects the local structure, i.e. the holes in the beam completely. The procedure to improve the 1D solution in the vicinity of the left hole is sketched in Fig. 5. The first contour plot shows the undisturbed stress rxx of the 1D beam solution (21) extended to the two-dimensional case (19). To compute the extended stresses (19) Poisson’s ratio in (20) has been set to m ¼ 0:3. Based on this solution, uo1 is computed on the overlay mesh according to Eq. (18). Finally, the composed solution rxx ðuhpd Þ ¼ rxx ðub Þ þ rxx ðuo1 Þ is shown at the bottom of Fig. 5. In Figs. 6 and 7 the hp–d results are compared with a two-dimensional reference solution, plotted along two different cut-lines. The location of the cut-lines is depicted in Fig. 5. The cut-line A–A corresponds to a fixed value of x = 4 m, intersecting the left circular hole. Cut-line B–B is located at x ¼ 4:5 m. The reference solution has been obtained by discretizing the whole beam as a 2D problem on a very fine mesh with high order elements. At the leftand the right-hand side of the two-dimensional mesh, homogeneous Dirichlet boundary conditions have been applied. From Figs. 6 and 7 it is evident that the hp–d
A. Du¨ster et al. / Comput. Methods Appl. Mech. Engrg. 196 (2007) 3524–3533
3529
Fig. 5. Improvement of dimensionally reduced approximation in the vicinity of the left hole.
The same approach can be used to ‘‘plaster’’ incompatible finite element meshes of n = 2 local sub-models. Each single overlay problem (22 and 23) can be computed thereby completely independently
Fig. 6. Stress component rxx along cut-line A–A.
Bðuo1 ; vo1 Þ ¼ Fðvo1 Þ Bðub ; vo1 Þ
8vo1 2 V o1 ;
ð22Þ
Bðuo2 ; vo2 Þ ¼ Fðvo2 Þ Bðub ; vo2 Þ
8vo2 2 V o2
ð23Þ
by ignoring the coupling terms Bðuo2 ; vo1 Þ and Bðuo1 ; vo2 Þ between the two overlay discretizations, see Eq. (15). Numerical investigations have demonstrated that ignoring the coupling terms still allows for a strong improvement of the beam solution. Fig. 8 illustrates this technique for the previously presented example. The solution at the hole on the right-hand side of the beam is computed on a different mesh. The two overlay meshes need not necessarily match on their boundaries and they can overlap. Again, a very good agreement of the hp–d approximation with the reference solution can be observed, even with only one block Gauss–Seidel iteration step. Remark 2. Numerical experiments have shown that the accuracy of the two-dimensional hp–d approximation can be further increased when improving the quality of the beam solution. If one computes the shear stresses from equilibrium conditions rather than from constitutive equations, the hp–d approximation coincides even better with the reference solution. One could choose also a Timoshenko beam theory instead, which includes already the shear stresses. However, we have decided to choose the simplest possible beam model as a starting point—which is of course more challenging—in order to demonstrate the capabilities of the proposed hp–d method.
Fig. 7. Stress component rxx along cut-line B–B.
3.2. Coupling 2D and 3D elasticity problems method significantly improves the accuracy of the beam solution. The model error is strongly reduced when compared to the two-dimensional reference solution.
The second numerical example to be considered is a three-dimensional thick-walled plate which is clamped on four sides (indicated by the shaded areas) and loaded with
3530
A. Du¨ster et al. / Comput. Methods Appl. Mech. Engrg. 196 (2007) 3524–3533
Fig. 8. Improvement of the 1D beam approximation by two independent overlay solutions, stress component rxx is plotted.
a pressure tz ¼ q acting only on a small part of the upper surface, Fig. 9. In order to solve this problem with the hp–d method, we start from the weak form of equilibrium based on the Reissner–Mindlin plate theory Bðu; vÞ ¼ FðvÞ 8v 2 V ;
ð24Þ
where the bilinear form and the linear functional define the internal and external virtual work, respectively. Vector T u ¼ ½bx ðx; yÞ; by ðx; yÞ; wðx; yÞ denotes the unknowns, where bx ðx; yÞ and by ðx; yÞ are the normal rotations in x–z and y–z planes, and wðx; yÞ corresponds to the deflection of the plate. The bilinear form reads
Bðu; vÞ ¼
t3 12
Z
eTB ðvÞCB eB ðuÞdX þ t
X
Z
eTS ðvÞCS eS ðuÞdX:
X
ð25Þ Based on the kinematic assumptions of the Reissner–Mindlin plate theory, the three-dimensional Cartesian displacement components, ux ðx; y; zÞ; uy ðx; y; zÞ and uz ðx; yÞ, are computed from the rotations and the deflection as u ¼ ½ux ¼ bx ðx; yÞz; uy ¼ by ðx; yÞz; uz ¼ wðx; yÞT : The strain term in (25),
Fig. 9. Three-dimensional plate with boundary conditions.
ð26Þ
A. Du¨ster et al. / Comput. Methods Appl. Mech. Engrg. 196 (2007) 3524–3533 T
are small. The constitutive equation describing the relation between strains and stresses is given by an isotropic linear elastic law
e ¼ ½exx ; eyy ; cxy ; cyz ; cxz T obx oby obx oby ow ow ; bx þ z; z; þ ¼ z; by þ oy ox ox oy oy ox
T
T
¼ ½eTB z; eTS
ð27Þ
denotes the (small) strains obtained by computing the corresponding derivatives of the displacement field. The stress–strain relation is given by the constitutive equations CB 0 r ¼ ½rxx ; ryy ; rxy ; ryz ; rxz T ¼ e; ð28Þ 0 CS where the matrices 2 1 m E 6m 1 CB ¼ 4 ð1 m2 Þ 0 0
0
3
0 7 5;
ð1mÞ 2
0 1
ð29Þ define the assumed isotropic linear elastic behaviour of the Reissner–Mindlin plate. The warping coefficient k ¼ 56 accounts for non-uniform shear distributions. The load functional Z qv dX ð30Þ FðvÞ ¼ X
takes into account the load q acting on the surface of the plate. The idea of the hp–d method with its application to dimensional adaptivity is to enhance the dimensionally reduced model, i.e. the Reissner–Mindlin plate theory, in those regions, where a significant deviation from the three-dimensional equations of linear elasticity are expected. We therefore apply the weak form of equilibrium of linear elasticity in 3D Bðu; vÞ ¼ FðvÞ
8v 2 V ;
ð31Þ
with Bðu; vÞ ¼ Z FðvÞ ¼
Z
r ¼ ½rxx ; ryy ; rzz ; rxy ; ryz ; rxz ¼ Ce; where 2 j þ 43 l 6 j 2 l 6 3 6 6 j 23 l 6 C¼6 0 6 6 4 0 0
1 Ek CS ¼ 2ð1 þ mÞ 0
eT ðvÞCeðuÞdX;
ð32Þ
X
vTt dC
ð36Þ
j 23 l j 23 l 0 j þ 43 l j 23 l 0 j 23 l j þ 43 l 0 0 0 l 0 0
0 0
0 0
3
0
0
0
07 7 7 07 7 07 7 7 05
0 0 l 0
ð37Þ
l
is the elasticity matrix, j ¼ E=ð3ð1 2mÞÞ, the bulk modulus, l ¼ E=ð2ð1 þ mÞÞ, the shear modulus. In a first step, an analysis based on the Reissner–Mindlin plate theory is performed by applying the p-version on a base mesh consisting of 108 quadrilaterals in combination q with the tensor product space Sp;p ps ðXst Þ [6,7,11] and a polynomial degree of p = 8, Fig. 10. The mesh is refined towards the boundary and at reentrant corners in order to resolve boundary layers and singularities, respectively. Following the procedure outlined above, the two-dimensional plate solution is then extended via the Reissner– Mindlin assumptions to a three-dimensional approximation yielding the displacements (26) as ub ¼ ½ux ¼ bx ðx; yÞz; uy ¼ by ðx; yÞz; uz ¼ wðx; yÞT
ð38Þ
and the corresponding strains ðub Þ ¼ ½exx ; eyy ; ezz ; cxy ; cyz ; cxz T T obx oby obx oby ow ow ; bx þ z; z; 0; þ ¼ ; z; by þ oy ox ox oy oy ox ð39Þ with ezz ¼ 0. In the vicinity of the area where the load is applied, an overlay mesh with 17 hexahedral elements is superimposed on the base mesh. An anisotropic Ansatz space Sp;p;q ðXhst Þ [6,7,11] in conjunction with the polynomial degree template
ð33Þ
CN
ð40Þ
as a basis for a finite element discretization. The function describing the Cartesian displacement components is denoted by u ¼ ½ux ðx; y; zÞ; uy ðx; y; zÞ; uz ðx; y; zÞT :
3531
ð34Þ
It is assumed that strains e ¼ ½exx ; eyy ; ezz ; cxy ; cyz ; cxz T T oux ouy ouz oux ouy ouy ouz oux ouz ; ; ; þ ; þ ; þ ¼ ox oy oz oy ox oz oy oz ox ð35Þ
is chosen for the three-dimensional computation, where the polynomial degree for the displacement components ux ; uy in in-plane is set to p = 7 and in transversal direction to q = 4. The deflection uz is approximated with an anisotropic Ansatz, where the polynomial degree in in-plane is set to p = 8 and in transversal direction to q = 5. The polynomial degrees in (40) have been chosen such that the discretization error of the anisotropic three-dimensional approximation is small. Since a hierarchic Ansatz space is applied, the quality of the approximation can be readily
3532
A. Du¨ster et al. / Comput. Methods Appl. Mech. Engrg. 196 (2007) 3524–3533
Fig. 10. 2D base and 3D overlay mesh for a local three-dimensional discretization.
controlled by increasing the polynomial degrees. An adaptive procedure to perform this adjustment in an automatic manner has been suggested, for example, in [15,36]. Applying the anisotropic discretization, the weak form Z Bðuo ; vo Þ ¼ Fðvo Þ eT ðvo ÞCeð ub ÞdX ð41Þ Xo
including the pseudo-load-vector resulting from the Reissner–Mindlin approximation is solved. The pseudo-loadvector is computed elementwise by interpolating the strains eð ub Þ from the extended Reissner–Mindlin model to the integration points of the hexahedral elements composing the overlay mesh. Having solved (41), the local–global three-dimensional hp–d solution is then defined as the hierarchical sum uhp–d ¼ ub þ uo of the extended global plate approximation ub and the local three-dimensional approximation uo obtained on the overlay mesh. As a result of the computations, consider the von Mises stress along the cut-lines C–C and D–D plotted in Figs. 11 and 12. The reference solution is based on a three-dimensional overkill approximation with a numerical error which is negligibly small. From Figs. 11 and 12 it is evident that the Reissner–Mindlin approximation coincides very well with the reference solution for a wide range of the domain. However, in the vicinity of the area where the load is
Fig. 12. von Mises stress along cut-line D–D.
applied, the two-dimensional approximation differs significantly from the reference solution. Comparing now the hp–d approximation with the reference solution, one observes that the quality of the finite element solution is enormously improved. Deviations from the reference solution are almost invisible. It is also important to notice that only one block Gauss–Seidel iteration step was performed in this example. Compared to a three-dimensional approximation of the whole problem, the additional numerical effort of the hp–d method is very small. 4. Conclusions
Fig. 11. von Mises stress along cut-line C–C.
The outlined method of combining dimensionally reduced approximations with higher dimensional models can be an attractive method leading to strongly improved solutions using only a small additional numerical effort. Further improvements in terms of efficiency and reliability could be achieved by applying an adaptive scheme in order to localize the domain where the dimensionally reduced approximation is to be extended by a higher dimensional Ansatz automatically. An error estimator could be used to quantify the model error. Areas with a high model error would therefore indicate where to enrich the dimensionally reduced approximation by means of a higher dimensional
A. Du¨ster et al. / Comput. Methods Appl. Mech. Engrg. 196 (2007) 3524–3533
Ansatz. Future work will focus on the development of such an adaptive scheme. It should finally be noted that, although the p-version of FEM was used for the examples presented in this paper, the method is not restricted to higher order finite element approximations. Low order FEM as available in any commercial finite element program can be used likewise. The only programming effort for implementing this procedure is the interpolation of strains obtained from the dimensionally reduced, global model to the extended local approximation and the computation of load vectors due to pre-strains. References [1] M. Vogelius, I. Babusˇka, On a dimensional reduction method. I. The optimal selection of basis function, Math. Comput. 37 (1981) 31–46. [2] M. Vogelius, I. Babusˇka, On a dimensional reduction method. II. Some approximation-theoretic results, Math. Comput. 37 (1981) 47– 48. [3] M. Vogelius, I. Babusˇka, On a dimensional reduction method. III. A posteriori error estimation and an adaptive approach, Math. Comput. 37 (1981) 361–384. [4] B. Szabo´, G. Sahrmann, Hierarchic plate and shell models based on p-extension, Int. J. Numer. Methods Engrg. 26 (1988) 1855–1881. [5] B. Szabo´, I. Babusˇka, Finite Element Analysis, John Wiley & Sons, 1991. [6] B. Szabo´, A. Du¨ster, E. Rank, The p-version of the finite element method, in: E. Stein, R. de Borst, T. Hughes (Eds.), Encyclopedia of Computational Mechanics, vol. 1, John Wiley & Sons, 2004, pp. 119– 139 (Chapter 5). [7] A. Du¨ster, H. Bro¨ker, E. Rank, The p-version of the finite element method for three-dimensional curved thin walled structures, Int. J. Numer. Methods Engrg. 52 (2001) 673–703. [8] A. Du¨ster, E. Rank, A p-version finite element approach for two- and three-dimensional problems of the J2 flow theory with non-linear isotropic hardening, Int. J. Numer. Methods Engrg. 53 (2002) 49–63. [9] A. Du¨ster, A. Niggl, V. Nu¨bel, E. Rank, A numerical investigation of high-order finite elements for problems of elasto-plasticity, J. Sci. Comput. 17 (2002) 429–437. [10] A. Du¨ster, A. Niggl, E. Rank, Thermo-elastic computations of geometrically non-linear three-dimensional thin-walled continua based on high order finite elements, in: H. Mang, F. Rammerstorfer, J. Eberhardsteiner (Eds.), Proceedings of the Fifth World Congress on Computational Mechanics, Vienna, Austria, 2002. [11] A. Du¨ster, High order finite elements for three-dimensional, thinwalled nonlinear continua, Ph.D. thesis, Lehrstuhl fu¨r Bauinformatik, Fakulta¨t fu¨r Bauingenieur- und Vermessungswesen, Technische Universita¨t Mu¨nchen.
3533
[16] A. Du¨ster, D. Scholz, E. Rank, Fluid–structure interaction using high-order anisotropic solid elements, in: Proceedings of the 5th International Conference on Computation of Shell and Spatial Structures, Salzburg, Austria, 2005. [17] D. Scholz, S. Kollmannsberger, A. Du¨ster, E. Rank, Thin solids for fluid–structure interaction, in: H. Bungartz, M. Scha¨fer (Eds.), Fluid–Structure Interaction, Modelling Simulation and Optimisation, Lecture Notes in Computational Science and Engineering, vol. 53, Springer, 2006, pp. 294–335. [18] C. Schwab, A-posteriori modeling error estimation for hierarchic plate models, Numer. Math. 74 (1996) 221–259. [19] J. Oden, J. Cho, Adaptive hpq-finite element methods of hierarchical models for plate- and shell-like structures, Comput. Methods Appl. Mech. Engrg. 136 (1996) 317–345. [20] E. Stein, S. Ohnimus, Coupled model- and solution-adaptivity in the finite-element method, Comput. Methods Appl. Mech. Engrg. 150 (1997) 327–350. [21] E. Stein, S. Ohnimus, Anisotropic discretization- and model-error estimation in solid mechanics by local Neumann problems, Comput. Methods Appl. Mech. Engrg. 176 (1999) 363–385. [22] E. Stein, M. Ru¨ter, Finite element methods for elasticity with errorcontrolled discretization and model adaptivity, in: E. Stein, R. de Borst, T. Hughes (Eds.), Encyclopedia of Computational Mechanics, vol. 2, John Wiley & Sons, 2004, pp. 5–58 (Chapter 2). [23] E. Rank, Adaptive remeshing and h-p domain decomposition, Comput. Methods Appl. Mech. Engrg. 101 (1992) 299–313. [24] E. Rank, A combination of hp-version finite elements and a domain decomposition method, in: Proceedings of the First European Conference on Numerical Methods in Engineering, Bru¨ssel, 1992. [25] E. Rank, A zooming-technique using a hierarchical hp-version of the finite element method, in: J. Whiteman (Ed.), The Mathematics of Finite Elements and Applications—Highlights 1993, Elsevier, Uxbridge, 1993. [26] R. Krause, Multiscale Computations with a Combined h- and p-version of the Finite-element Method, Ph.D. Thesis, Fach Numerische Methoden und Informationsverarbeitung, Universita¨t Dortmund, 1996. [27] R. Krause, E. Rank, Multiscale computations with a combination of the h- and p-versions of the finite-element method, Comput. Methods Appl. Mech. Engrg. 192 (2003) 3959–3983. [28] J. Fish, The s-version of the finite element method, Comput. Struct. 43 (3) (1992) 539–547. [29] J. Fish, A. Nath, Adaptive and hierarchical modelling of fatigue crack propagation, Int. J. Numer. Methods Engrg. 36 (1993) 2825–2836. [30] J. Fish, S. Markolefas, Adaptive s-method for linear elastostatics, Comput. Methods Appl. Mech. Engrg. 104 (1993) 363–396. [31] J. Fish, S. Markolefas, Adaptive global–local refinement strategy based on the interior error estimates of the h-method, Int. J. Numer. Methods Engrg. 37 (1994) 827–838. [32] J. Fish, M. Pandheeradi, V. Belsky, An efficient multilevel solution scheme for large scale non-linear systems, Int. J. Numer. Methods Engrg. 38 (1995) 1597–1610. [33] E. Rank, R. Krause, A-posteriori-Fehlerscha¨tzung bei Multi-SkalenProblemen, ZAMM 75 (1995), SI. [34] E. Rank, R. Krause, A multiscale finite-element-method, Comput. Struct. 64 (1997) 139–144. [35] A. Du¨ster, E. Rank, G. Steinl, W. Wunderlich, A combination of an h- and p-version of the finite element method for elastic-plastic problems, in: Proceedings of ECCM ’99, European Conference on Computational Mechanics, Mu¨nchen, Germany, 1999. [36] D. Scholz, An anisotropic p-adaptive method for linear elastostatic and elastodynamic analysis of thin-walled and massive structures, Dissertation, Lehrstuhl fu¨r Bauinformatik, Fakulta¨t fu¨r Bauingenieur- und Vermessungswesen, Technische Universita¨t Mu¨nchen, 2006.