Computational Mechanics https://doi.org/10.1007/s00466-018-1585-6

ORIGINAL PAPER

Numerical procedure to couple shell to solid elements by using Nitsche’s method Takeki Yamamoto1

· Takahiro Yamada2

· Kazumi Matsui2

Received: 15 December 2017 / Accepted: 21 May 2018 © Springer-Verlag GmbH Germany, part of Springer Nature 2018

Abstract This paper presents a numerical procedure to couple shell to solid elements by using the Nitsche’s method. The continuity of displacements can be satisfied approximately with the penalty method, which is effective in setting the penalty parameter to a sufficiently large value. When the continuity of only displacements on the interface is applied between shell and solid elements, an unreasonable deformation may be observed near the interface. In this work, the continuity of the stress vector on the interface is considered by employing the Nitsche’s method, and hence a reasonable deformation can be obtained on the interface. The authors propose two types of shell elements coupled with solid elements in this paper. One of them is the conventional MITC4 shell element, which is one of the most popular elements in engineering applications. This approach shows the capability of discretizing the domain of the structure with the different types of elements. The other is the shell element with additional degrees of freedom to represent thickness–stretch developed by the authors. In this approach, the continuity of displacements including the deformation in the thickness direction on the interface can be considered. Several numerical examples are presented to examine the fundamental performance of the proposed procedure. The behavior of the proposed simulation model is compared with that of the whole domain discretized with only solid elements. Keywords Nitsche’s method · Combined modeling · Shell element · Solid element

1 Introduction Structures are constructed with various shape of components which can be modeled as solids, plates, and beams. In the finite element analysis for such structures, the whole domain is often discretized with only one type of element. This approach does not seem to be appropriate in general cases, and it may exhibit difficulty in some cases. For example, in the structure composed of an assemblage of solids and plates, the whole domain is discretized with continuum elements, the

B

Takeki Yamamoto [email protected] Takahiro Yamada [email protected] Kazumi Matsui [email protected]

1

Department of Civil and Environmental Engineering, Tohoku University, Aza-Aoba 6-6-06, Aramaki, Aoba-ku, Sendai 980-8579, Japan

2

Graduate School of Environment and Information Science, Yokohama National University, Tokiwadai 79-7, Hodogaya-ku, Yokohama, Kanagawa 240-8501, Japan

local behavior can be evaluated effectively. However, huge computational cost may be required. On the other hand, if the entire region is discretized with structural elements, global bending behavior of thin-walled structures can be predicted efficiently, but they are not sufficient to evaluate the local behaviors. Thus, the effective procedure of discretizing the whole domain of the problem of interest is to select different types of elements according to the shape and the behavior. For example, the thin-walled portion can be discretized with shell elements, whereas the region in which details of the structural response needs to be should be modeled with solid elements. For the purpose of modeling such structures flexibly in the finite element analysis, it is necessary to connect shell and solid elements properly. In numerical procedures to connect structural and continuum elements, it is important to combine the properties of two types of elements on the connecting interface. On the interface, the nodal degrees of freedom of structural and continuum elements are incompatible with each other, and these elements cannot be connected directly. Several approaches that focus on combining the different types of elements are roughly classified into three groups. One is a class of intro-

123

Computational Mechanics

ducing the transition element that combines the properties on the interface. Surana presented isoparametric transition elements [1], which was developed for connecting thin-walled cross section and solid-like portion, and developed for linear elastic axisymmetric analysis [2] and three dimensional stress analysis [3]. Bathe and Ho [4] extended the classical C0 compatible displacement-rotation curved shell element to model the intersections of both shell and solid elements. In this approach, the formulation is restricted to static analysis and based on a stress-strain relationship of the shell element. Surana also proposed transition elements for nonlinear analyses [5,6]. In order to avoid the locking phenomena when transition elements are applied to the modeling of the interface, transition elements were improved. Cofer and Will [7] presented the transition element, which is connected quadratic shell and solid elements. Gmür and Schorderet [8] also developed the element for coupling superparametric shell elements and isoparametric solid elements. The transition element was applied to thermo-elastic problems by Chavan and Wrigger [9]. Furthermore, Gmür and Kauten [10] used transition elements, which were constructed from the approach [8], in the application of dynamic problems. Attempts have been made to connect hexahedral and tetrahedral elements for uniform strain by Dohrman and Key [11], and to dissimilar three dimensional finite element meshes by Dohrmann et al. [12]. However, transition elements, in which the displacement field is expressed by a combination of displacements of nodes of structural and continuum elements, show the unreasonable deformation, because the discrepancy of the displacement in the thickness direction between structural and continuum elements. Hence, Garusi and Tralli [13] presented transition elements constructed by the hybrid stress method, which assumes a stress state of the interface, for the modeling solid-to-beam and plate-to-beam structures. Ahn and Basu applied a p-convergent transition element [14] to the modeling of patch repaired plates [15], and to connect shell and solid elements [16]. In these transition elements, it is necessary to verify the domain modeled by transition elements in addition to regions discretized with structural and continuum elements. Thus, in the evaluation of simulation models using transition elements, it is important to assess the relationship of the response of structures and the behavior of transition elements. The second group is a class of coupling methodology based on the multipoint constraint (MPC) equations [17]. The multipoint constraint equations are more simpler approach for connecting different types of elements than introducing transition elements. Curiskis and Valliappan [18] proposed a technique for incorporating constraint equations, which define the relationship between various degrees of freedom, into linear algebraic systems. Abel and Shephard [19] introduced general constraint equations into finite element analysis, and Shephard [20] presented a procedure for the

123

application of linear multipoint constraints derived from the transformation approach, which reduces the number of equations to be solved by the number of constraints. When the multipoint constraint equations are employed, the unreasonable deformation can be observed for connecting structural and continuum elements in large strain ranges [21]. Hence, Jiain et al. [22] developed a connecting procedure for shell and solid elements by imposing the multipoint constraint equations, in which the additional degrees of freedom to represent the change in thickness of shell element is considered. Numerical procedures imposing constraint conditions of displacement continuity on the interface by the Lagrange multiplier method or the penalty method are widely used [23]. The Lagrange multiplier method, in which constraint equations and unknown variables are added to equilibrium equations, increases the number of unknowns in the derived linear system of equations. On the other hand, the penalty method, in which constraint equations are approximately satisfied without introducing additional variables, works when the penalty parameter is set to a sufficiently large value. In order to obtain the reasonable deformation on the interface connecting structural and continuum elements by using these approaches, additional techniques, in which represents the deformation in the thickness direction of structural elements, must be introduced. The final group of connecting different types of elements is to utilize the equilibrium equation for the work on the interface. Monaghan et al. [24], McCune et al. [25], and Shim et al. [26] presented the modeling to combine a reduced or lower dimensional element with a higher dimensional element. These approaches are based on the multipoint constraint equations to achieve compatibility of displacements and stress equilibrium on the interface between different types of elements. Robinson et al. [27] proposed a numerical procedure to create mixed dimensional models automatically. Aminpour and Krishnamurthy [28] and Aminpour et al. [29] developed the interface elements, which are coupled finite element meshes without any restrictions of the size and fitted or unfitted geometry. The formulation of interface elements is simpler than that of transition elements. However, it is inevitable to increase the number of unknowns in the derived linear system of equations, because these techniques are based on the Lagrange multiplier method. Osawa et al. [30] also proposed the technique for coupling shell to solid elements by introducing a fictitious shell plane perpendicular to the original shell plane. Numerical results obtained by this approach are in good agreement with those obtained from the multipoint constraint equations. In this approach, the degrees of freedom on the interface are increased, similar to the interface element [28,29]. Mata et al. [31] developed a two-scale approach (global and local) corresponding to the nonlinear response of structures. This approach evaluates the overall behavior of simulation models at global

Computational Mechanics

scale and the local behavior at local scale. Hence, neither the multipoint constraint equations nor Lagrange multipliers are employed. However, it is necessary to apply additional iterative schemes, which considers the interaction between scales, in addition to an iterative procedure at each scale. Recently, the Nitsche’s method [32], which is a technique to connect numerical solutions calculated in divided domains, has been focused. This technique was originally proposed to enforce Dirichlet boundary conditions weakly as an approach to equivalent pointwise constraints. The method was applied to connect the solutions for interface problems [33,34] and overlapping meshes [35,36], for domain decomposition methods [37,38], for coupling nonconforming patches in isogeometric analysis [39], and for meshless analysis [40]. The Nitsche’s method was also used in the interior discontinuous Galerkin method [41], and shows its effectiveness of the modeling the discontinuity problems [42] and the contact problems [43]. Furthermore, numerical procedures for connecting different types of elements were proposed by applying the Nitsche’s method to the domain decomposition methods [39,44]. In this paper, we propose a numerical procedure to couple shell to solid elements by using the Nitsche’s method. The proposed formulation is based on the domain decomposition methods, in which the whole domain of the problem of interest is discretized with shell and solid elements. On the interface, the penalty method is applied to satisfy the continuity of displacements, and the continuity of stress vectors is imposed by employing the Nitsche’s method. Hence, no additional degree of freedom on the interface is needed in this approach. The authors propose two types of shell elements coupled with solid elements in this paper. One of them is the conventional MITC4 shell element [45], which is one of the most popular elements in engineering applications. This approach shows the capability of discretizing the domain of the structure with the different types of elements. The other is the modified MITC4 shell element with additional degrees of freedom to represent thickness–stretch, which was developed by the authors [46], and is designated as MITC4ts shell element in this paper. In this approach, the continuity of displacements including the deformation in the thickness direction on the interface can be considered. As a key technique of this study, the MITC4ts shell element has a capability of evaluating the change in the thickness and using the three dimensional constitutive equations without any additional assumptions. Due to the formulation of the MITC4ts shell element, the thickness–stretch can be discretized at each element. Therefore, the discontinuity of the deformation in the thickness direction on the interface is allowed in the proposed connecting procedure. Thus, static condensation can be employed in constructing the stiffness matrix on the connecting interface, in which the degrees of freedom

to represent thickness–stretch in the MITC4ts shell element are condensed out. In the case of connecting the MITC4ts shell and the hexahedral elements, it is important to verify the effectiveness for considering the deformation in the thickness direction on the interface and applying the three dimensional constitutive equations to the whole domain. Some representative numerical examples are presented to verify the numerical procedure to couple shell to solid elements in elastic deformation problems. In small and large strain ranges, the behavior of the proposed simulation model is compared with that of the simulation model discretized with only solid elements.

2 Formulation This section describes the governing equations of the decomposed domains, which are discretized with shell and solid elements. On the connecting interface between shell and solid elements, the continuity conditions for displacement and stress vectors must be considered. In this approach, the continuity condition of displacement vector is satisfied approximately by the penalty method. On the other hand, the continuity condition of stress vector can be imposed by employing the Nitsche’s method.

2.1 Governing equations As shown in Fig. 1, the whole domain of the problem can be decomposed two domains discretized with both shell and solid elements. The whole domain  with boundary  can be divided into two non-overlapping regions, which are the combination of the domain of the continuum part (c) with boundary  (c) and that of the structural part (s) with boundary  (s) . Here, the superscripts (c) and (s) denote quantities defined in the domain discretized with solid elements and shell elements, respectively. On the interface connecting (c) and (s) , the boundary is denoted by  (cs) =  (c) ∩  (s) . In addition, the Dirichlet and Neumann boundaries are represented as u and t , respectively. Hence, the definitions of the boundaries in each region are expressed as

Domain of solid

Domain of shell

Fig. 1 Domain decomposition for connecting shell and solid elements

123

Computational Mechanics (c)

energy  associated with the displacement fields u(c) in (c) and u(s) in (s) is derived from (1) to (6) as

u(c) ∪ t ∪  (cs) =  (c) , (c)

u(c) ∩ t = ∅, u(c) ∩  (cs) = ∅, (c) t ∩  (cs) (s) u(s) ∪ t ∪  (cs) (s) u(s) ∩ t u(s) ∩  (cs) (s) t ∩  (cs)



 u , u

= ∅,

(c)

(s)



 =

=  (s) ,

0 (c)



−

= ∅,

  W (c) F (c) d

0 (c)



 u(c) · b(c) d −

  W (s) F (s) d 0 (s)   (s) (s) − u · b d −

= ∅, = ∅.

∇·σ

+b

∇·σ

(s)

+b

σ

=0

in  ,

(1)

(s)

=0

(s)

in  ,

(2)

(c)

= u¯

(c) t , (s) t , (c) t , (s) t ,

(3)

(c)

(c)

on

u(s) = u¯ (s)

on

¯ (c)

on

¯ (s)

on

u

(c) (c)

n

=t

σ (s) n(s) = t

(4) (5) (6)

where u denotes the displacement vector, σ is the Cauchy stress tensor calculated from the displacement field u, and b and t are the body force and the surface traction vectors, respectively. In this approach, the equilibrium equations (1), (2) are considered in each domain. The prescribed displacement and traction vectors are expressed as u¯ and ¯t , respectively. The outward unit normal vectors to the interfaces  (c) and  (s) are also represented as n(c) and n(s) , respectively. On the connecting interface  (cs) , the outward unit normal vector n can be defined as n = n(c) = −n(s) .

(7)

In addition, the continuity conditions on  (cs) can be written as u(c) = u(s) σ

(c)

n=σ n (s)

on  (cs) ,

(8)

on 

(9)

(cs)

.

For the purpose of solving the equilibrium equations (1), (2), the continuity conditions (8), (9) must be imposed on  (cs) , in addition to the Dirichlet boundary conditions (3), (4) and the Neumann boundary conditions (5), (6).

2.2 Weak formulation In this work, we consider the structure of which material is hyperelastic. For a hyperelastic material, the total potential

123

(c)

d

0  (s)

(s) u(s) · ¯t d,

(10)

The governing equations in  can be described using the domain decomposition method as (c)

u(c) · ¯t

+

0 (s)

(c)

0  (c)

where 0  is the physical region with the boundary 0  in the initial configuration, F is the deformation gradient tensor calculated from the displacement field u, and W is the stored energy function defined per unit reference volume. In this approach, the continuity condition (8) is imposed approximately by the penalty method. On the interface  (cs) , a gap vector for the displacement fields is expressed as g = u(c) − u(s) .

(11)

By using the gap vector (11) as a penalty term, the potential can be written as    W P u(c) , u(s) =

0  (cs)

p g · g d, 2

(12)

where p is the penalty parameter, which is a suitably large positive number having the physical interpretation of a stiff spring constant. In the penalty method, the constraint (8) is enforced by adding the potential W P in (12) to the function  in (10). Hence, the total potential energy function in the penalty method can be obtained as       P u(c) , u(s) =  u(c) , u(s) + W P u(c) , u(s) . (13) If p is large enough, the minimum of  P is attained with satisfying the constraint condition (8) approximately. The Nitsche’s method can be employed to impose the constraint condition (9). When the continuity conditions for both displacement fields and stress vectors are satisfied on the interface  (cs) , the work applied in terms of the surface traction on the connecting interface is equal to zero. From (7), (8), and (9), this relation can be expressed as 

 t  (cs)

u(c) · σ (c) n(c) d +

t  (cs)

u(s) · σ (s) n(s) d = 0. (14)

Computational Mechanics

By pulling back the relation (14) to the initial configuration, the work on the interface can be rewritten as   (c) (c) (c) u · P N d + u(s) · P (s) N (s) d = 0, (15) 0  (cs)

0  (cs)

where P is the first Piola–Kirchhoff stress tensor computed using the displacement vector u, N is the outward unit normal vector to the interface  (cs) in the initial configuration. Similarly to the relation (7) in the current configuration, the outward unit normal vector to the interface  (cs) are described as

the variational equation can be derived as   G (c) u(c) , u(s) ; δu(c)  = S(c) : δ E (c) d 0 (c)   ¯t (c) · δu(c) d − b(c) · δu(c) d − t (c) t  (c)      + δu(c) · p u(c) − u(s) − [[ P]] N d 0  (cs)    −α u(c) − u(s) · δ P (c) N d, 0  (cs)

N = N (c) = −N (s) .

(16)

From the work (15) and the normal vector (16) in the initial configuration, the potential using the Nitsche’s method can be derived as    W N u(c) , u(s) = −

 0  (cs)

 u(c) − u(s) · ([[ P]] N) d, (17)

where [[ P]] = α P (c) + (1 − α) P (s) , and α is the Nitsche’s parameter, where 0 ≤ α ≤ 1. Hence, the total potential energy function using the Nitsche’s method can be obtained adding the potential W N in (17) to the function P in (13) as       N u(c) , u(s) =P u(c) , u(s) + W N u(c) , u(s)     = u(c) , u(s) + W P u(c) , u(s)   + W N u(c) , u(s)    = W (c) F (c) d 0 (c)   (c) (c) (c) − u · b d − u(c) · ¯t d 0 (c) 0  (c)    + W (s) F (s) d 0 (s)   (s) − u(s) · b(s) d − u(s) · ¯t d 0 (s) 0  (s)  p + g · g d 0  (cs) 2    − u(c) − u(s) · ([[ P]] N) d. 0  (cs)

(18) By taking the variation of the potential energy (18) with respect to the admissible displacement field of u(c) in (c) ,

(19)

where S and E are the second Piola–Kirchhoff stress tensor and the Green–Lagrange strain tensor, respectively, and δu(c) is the virtual displacement vector in (c) . Similarly, the variational equation for the admissible displacement field of u(s) in (s) can be expressed as   G (s) u(c) , u(s) ; δu(s)  = S(s) : δ E (s) d 0 (s)   (s) (s) ¯t (s) · δu(s) d − b · δu d − t (s) t  (s)      − δu(s) · p u(c) − u(s) − [[ P]] N d 0  (cs)    − (1 − α) u(c) − u(s) · δ P (s) N d, 0  (cs)

(20)

where δu(s) is the virtual displacement vector in (s) . The incremental form of the variational equations (19) and (20) at time t + t are shown to be   G (c) u(c) + u(c) , u(s) + u(s) ; δu(c)      S(c) + S(c) : δ E (c) + E (c) d = 0 (c)   ¯t (c) · δu(c) d − b(c) · δu(c) d − t (c) t  (c)      +p δu(c) · u(c) + u(c) − u(s) + u(s) d 0 (cs)   δu(c) · [[ P +  P]] N d − 0  (cs)      u(c) + u(c) − u(s) + u(s) · δ P (c) N d −α 0  (cs)

=0

(21)

and   G (s) u(c) + u(c) , u(s) + u(s) ; δu(s)      S(s) + S(s) : δ E (s) + E (s) d = 0 (s)   ¯t (s) · δu(s) d − b(s) · δu(s) d − t (s)

t  (s)

123

Computational Mechanics  −p  +

0  (cs)

δu(s) ·



   u(c) + u(c) − u(s) + u(s) d

δu(s) · [[ P +  P]] N d      − (1 − α) u(c) + u(c) − u(s) + u(s) · δ P (s) N d 0  (cs)

In the above equations, the terms of the stiffness and the equivalent nodal force corresponding to the potentials for both the penalty method and the Nitsche’s method are added to the variational formulations for the displacement method.

0  (cs)

= 0,

(22)

where u is the incremental displacement vector from time t to t +t. Arranging the known terms and the unknown terms of the weak forms (21) and (22), the incremental variational formulations can be rewritten as     S(c) : δ E (c) d + S(c) : δ E (c) d 0 (c) 0 (c)    +p δu(c) · u(c) − u(s) d 0 (cs)   − δu(c) · [[ P]] N d 0  (cs)    −α u(c) − u(s) · δ P (c) N d 0  (cs)   (c) (c) (c) = δu · b d + δu(c) · ¯t d t (c) 0  (c)  − S(c) : δ E (c) d 0 (c)    −p δu(c) · u(c) − u(s) d 0 (cs)   + δu(c) · [[ P]] N d 0  (cs)    +α u(c) − u(s) · δ P (c) N d (23)

3 Discretization This section describes the numerical procedures to couple shell to solid elements. In this study, we propose two types of shell elements coupled with solid elements. One of them is the MITC4 shell element [45], in which the thickness is constant and a plane stress state is assumed in the transverse direction. The other is the MITC4ts shell element, in which the thickness–stretch can be evaluated using the additional degrees of freedom [46]. Note that the thickness–stretch is discretized at each element. Hence, the discontinuities for the thickness on the connecting interface can be allowed in this approach. Thus, static condensation can be employed in constructing the stiffness matrix on the connecting interface. As a key technique of this approach, the MITC4ts shell element has a capability of using the three dimensional constitutive equations without any assumption. Thanks to such a property of the MITC4ts shell element, the same constitutive equation can be applied uniformly to the whole domain. After discretizing the displacement fields with each element, we present the stiffness equation for connecting shell and hexahedral elements.

3.1 Displacement fields

0  (cs)

and 

In this work, the geometry and kinematics are described in a Cartesian coordinate system. The local coordinates ξ , η, and ζ are defined in an element. Note that (ξ, η) are the in-plane coordinates, while ζ denotes the transverse coordinate.



  S : δ E d + S(s) : δ E (s) d 0 (s) 0 (s)    −p δu(s) · u(c) − u(s) d 0 (cs)   + δu(s) · [[ P]] N d 0  (cs)    − (1 − α) u(c) − u(s) · δ P (s) N d 0  (cs)   (s) = δu(s) · b(s) d + δu(s) · ¯t d t (s) 0  (s)  − S(s) : δ E (s) d 0 (s)    +p δu(s) · u(c) − u(s) d 0 (cs)   − δu(s) · [[ P]] N d 0  (cs)    + (1 − α) u(c) − u(s) · δ P (s) N d. (s)

(s)

0  (cs)

123

3.1.1 Domain of continuum part A hexahedral element is employed to discretize the domain of the continuum part of the simulation model. The displacement fields of the hexahedral element can be expressed as 8      Na(c) ξ (c) , η(c) , ζ (c) ua(c) u(c) ξ (c) , η(c) , ζ (c) = a=1

  = Nc ξ (c) , η(c) , ζ (c) uc ,

(24)

(25)

where ua denotes the nodal displacement vector at node a, Na is the shape function described in local coordinates, and N and u represent the matrix of shape functions and the vector of nodal displacements, respectively. Here, the subscript c denotes quantities defined in solid elements.

Computational Mechanics

Deformed configuration

Deformed configuration

Initial configuration

Initial configuration

Hexahedral solid

MITC4 shell

Hexahedral solid

Fig. 2 Displacement field on connecting surface both MITC4 shell and hexahedral solid

3.1.2 Domain of structural part In this study, two approaches for modeling the domain of the structural part discretized with shell elements are presented. One is to use the MITC4 shell element, which is proposed to avoid the locking phenomena of the degenerated shell element [23]. The displacement fields of the conventional quadrilateral shell element are defined by u (ξ, η, ζ ) =

4 

Na (ξ, η) ua +ζ

a=1

4  ha a=1

2

Na (ξ, η) (v a −V a ) ,

MITC4ts shell

Fig. 3 Displacement field on connecting surface both MITC4ts shell and hexahedral solid

introduced to represent thickness–stretch. The formulation of the MITC4ts shell element is based on the MITC4 shell element, and assumed to be constant of the thickness in each element. The geometry, other than for the thickness–stretch, is represented with midsurface nodes, while the thickness– stretch is described with additional nodes, which have only the displacement variations in the thickness direction. Therefore, the displacement fields of the MITC4ts shell element can be expressed as 4      u(s) ξ (s) , η(s) , ζ (s) = Na(s) ξ (s) , η(s) ua(s) a=1

where ua denotes the nodal displacements at the midsurface node a, V a and v a are the director vectors at the midsurface node a in the initial configuration and the deformed configuration, respectively, and h a is the nodal thickness. In this study, we assume a constant shell thickness for each element, hence the displacement fields of the quadrilateral shell element can be described as   u(s) ξ (s) , η(s) , ζ (s) =

4 

  Na(s) ξ (s) , η(s) ua(s)

a=1

   he Na(s) ξ (s) , η(c) v a(s) − V a(s) 2 a=1   = Ns ξ (s) , η(s) , ζ (s) us , (26) +ζ

4 

where h e is the initial thickness of the element, which is constant in the deformed configuration. Here, the subscript s denotes quantities defined in shell elements. On the other hand, the other is to employ the MITC4ts shell element, in which additional degrees of freedom are

+ζ

4  h e  (s)  (s) (c)   (s) Na ξ , η v a − V a(s) 2 a=1



+ w ζ (s)

4 

Na(s) (0, 0)v a(s)

a=1

    (s) = Ns ξ (s) , η(s) , ζ (s) us + w ζ (s) v td     = Ns ξ (s) , η(s) , ζ (s) us + Nt ζ (s) ut

(27)

where w is the thickness variation assumed to be constant on the ξ (s) −η(s) plane in each element, and v td is the vector in the thickness direction, which is defined at each element. Note that the subscript t denotes quantities involving the thickness– stretch of the MITC4ts shell element. 3.1.3 Continuity condition on connecting interface The continuity condition (8) can be expressed by using the displacement fields (25), (26) for connecting the MITC4 shell element and the hexahedral element (see Fig. 2) as

123

Computational Mechanics     Nc ξ (c) , η(c) , ζ (c) uc = Ns ξ (s) , η(s) , ζ (s) us

 (cs) ,

on

(28)

while by using the displacement fields (25), (27) for coupling the MITC4ts shell element to the hexahedral element (see Fig. 3) as       Nc ξ (c) , η(c) , ζ (c) uc = Ns ξ (s) , η(s) , ζ (s) us + Nt ζ (s) ut

S(c) : δ E (c) d = {δuc }t kcc uc ,

0 (c)

(29)

In this study, the continuity conditions (28), (29) are satisfied approximately by the penalty method.

 p

0  (cs)

−p

0  (cs)

and that of the MITC4 shell element is given by     u(s) ξ (s) , η(s) , ζ (s) = Ns ξ (s) , η(s) , ζ (s) us .

=

rs + fsP + fsN rc + fcP + fcN

,

P δu(s) · u(c) d = {δus }t ksc uc ,

 P δu(s) · u(s) d = {δus }t kss us ,

δu(s) · b(s) d +

0  (cs)

δu

· P

(c)

N d

0  (cs)

N u(c) · δ P (s) N d = {δus }t ksc uc ,

 (1 − α)

0 (s)

 p

0  (cs)

δu(s) ·  P (s) N d

+ (1 − α)

0  (cs)

u

· δP

(s)

N d = {δus }

t

N kss us ,

(32)

123

(s)

d

0  (cs)

S(s) : δ E (s) d = {δus }t rs ,

  δu(s) · u(c) − u(s) d = {δus }t fsP ,

0  (cs)

0  (cs)

  δu(c) · u(c) − u(s) d = {δuc }t fcP ,

(37)

 δu(c) · [[ P]] N d    u(c) − u(s) · δ P (c) N d = {δuc }t fcN , +α

0  (cs)

0  (cs)

 (s)

δu(s) · ¯t



−p

 − (1 − α)

0  (s)

 −

 (s)

N u(s) · δ P (c) N d = {δuc }t kcs us ,

0 (c)

 α

0  (cs)

δu(s) · [[ P]] N d    u(c) − u(s) · δ P (s) N d = {δus }t fsN , + (1 − α) 0  (cs)   (c) (c) (c) δu · b d + δu(c) · ¯t d t (c) 0  (c)  S(c) : δ E (c) d = {δuc }t rc , (36) −



0  (cs)

δu(c) ·  P (s) N d



−

0 (s)

p

(35)





where  S(s) : δ E (s) d = {δus }t kss us ,

0  (cs)

0  (cs)



(31)

−p

N u(c) · δ P (c) N d = {δuc }t kcc uc ,



The discretized equation corresponding to the incremental variational formulations (23), (24) can be described as 

δu(c) ·  P (c) N d

0  (cs)

− (1 − α)

t (s)

us uc

(34)



+α



P δu(c) · u(s) d = {δuc }t kcs us ,



−α

From above definitions of displacement fields, the incremental displacement vector of the hexahedral element can be expressed as     (30) u(c) ξ (c) , η(c) , ζ (c) = Nc ξ (c) , η(c) , ζ (c) uc ,

P + kN P + kN ksc kss + kss ss sc P + kN P + kN kcs k + kcc cc cs cc

P δu(c) · u(c) d = {δuc }t kcc uc ,

0  (cs)

−α

3.2 Coupling of MITC4 shell element and solid element



(33)



 (cs) .

on



(38)

in which the superscript t denotes the transpose. Note that the stiffness matrix in the discretized equation (31) may be symmetry only when the Nitsche’s parameter α is equal to 0.5.

Computational Mechanics



3.3 Connecting MITC4ts shell element and solid element

t (s)

The incremental displacement vector of the hexahedral element is defined in (30). From the definition of the displacement fields (27), the incremental displacement vector of the MITC4ts shell element can be expressed as u

(s)

 δu(s) · b(s) d + 

−  p

0 (s)

0  (cs)

0  (s)

δu(s) · ¯t (s) d

S(s) : δ E (s) d = {δus }t rs + {δut }t rt ,

  δu(s) · u(c) − u(s) d = {δus }t fsP + {δut }t ftP ,



    ξ (s) , η(s) , ζ (s) = Ns ξ (s) , η(s) , ζ (s) us   + Nt ζ (s) ut .

−

0  (cs)

δu(s) · [[ P]] N d

+ (1 − α)

The discretized equation corresponding to the incremental variational formulations (23), (24) can be described as ⎫ ⎤⎧ P + kN k + kP + kN P + kN ksc kss + kss ⎨ us ⎬ ss st sc st st P + kN ⎣ k(s) + kP + kN ktt + kP + kN ⎦ ut ktc ts ts ts tt tt tc ⎩ ⎭ P + kN P + kN P + kN uc kcs k k + k cc cs cc cc ct ct ⎫ ⎧ ⎨ rs + fsP + fsN ⎬ = rt + ftP + ftN , (39) ⎭ ⎩ rc + fcP + fcN



 0  (cs)

 u(c) − u(s) · δ P (s) N d

= {δus }t fsN + {δut }t ftN .

⎡

where



0 (s)

S(s) : δ E (s) d = {δus }t kss us + {δus }t kst ut + {δut }t kts us + {δut }t ktt ut ,

 −p

P P δu(s) · u(c) d= {δus }t ksc uc +{δut }t ktc uc ,

0  (cs)

 p

0  (cs)

 α

0  (cs)

P δu(s) · u(s) d = {δus }t kss us + {δus }t kstP ut

+ {δut }t ktsP us + {δut }t kttP ut ,  δu(s) ·  P (c) N d−(1−α) u(c) · δ P (s) N d 0  (cs)

N N = {δus }t ksc uc + {δut }t ktc uc ,

(40)

 (1 − α)

0  (cs)

+ (1 − α)

 u(s) · δ P (s) N d

N u + {δu }t kN u , + {δut }t kts s t t tt

 0  (cs)

P u +{δu }t kP u , δu(c) · u(s) d= {δuc }t kcs s c t ct

 − (1 − α)

0  (cs)

 δu(c) ·  P (s) N d + α

N u +{δu }t kN u , = {δuc }t kcs s c t ct

(41)

where

k¯ ts = kts + ktsP + ktsN , k¯ tt = ktt + kttP + kttN , P N k¯ tc = ktc + ktc ,

N u + {δu }t kN u = {δus }t kss s s t st

−p

⎫ ⎧ ⎫ ⎤⎧ k¯ ss k¯ st k¯ sc ⎨ us ⎬ ⎨ r¯ s ⎬ ⎣ k¯ ts k¯ tt k¯ tc ⎦ ut = r¯ t , ⎩ ⎭ ⎩ ⎭ r¯ c uc k¯ cs k¯ ct k¯ cc ⎡

P N k¯ ss = kss + kss + kss , P N k¯ st = kst + kst + kst , P N + ksc , k¯ sc = ksc

δu(s) ·  P (s) N d

0  (cs)

P , kN , r , f P , and f N are defined in (33), (34), Here, kcc , kcc cc c c c (35), (36), (37), and (38), respectively. Note that the difference between kss in (32) and kss in (40) is the constitutive equation, which is employed to construct the stiffness matrices. The constitutive equation, in which the stiffness matrix (32) is constructed by using the MITC4 shell element, is modified by the assumption of a plane stress state in the transverse direction. On the other hand, the MITC4ts shell element has a capability of employing the three dimensional constitutive equations without any assumption. One of the most important features is that the same constitutive equation can be used in the whole domain, which is discretized with both shell and solid elements. Note that the stiffness matrix in the discretized equation (39) may be symmetry only when the Nitsche’s parameter α is equal to 0.5. Hence, the discretized equation (39) is rewritten as

u(s) · δ P (c) N d

0  (cs)

P N + kcs , k¯ cs = kcs P N ¯kct = kct + kct , P N + kcc , k¯ cc = kcc + kcc

r¯ s = rs + fsP + fsN , r¯ t = rt + ftP + ftN , r¯ c = rc + fcP + fcN .

123

Computational Mechanics

Global coordinates

Parametric coordinates

Parametric coordinates

Layer coordinates

Integration point for penalty terms

Integration point for penalty terms

Integration point for Nitsche’s terms

Integration point for Nitsche’s terms

Fig. 4 Numerical integration on connecting surface both shell and solid elements

Due to the formulation of the MITC4ts shell element, the thickness–stretch can be discretized at each element. Hence, the discontinuity of the deformation in the thickness direction on the interface is allowed in the proposed connecting procedure. The displacement variations, which represent thickness–stretch in the shell element, are independent of the displacement fields described by the hexahedral element. Thus, the stiffness matrices k¯ tc and k¯ ct are omitted from the discretized equation (41) as ⎡

k¯ ss ⎣ k¯ ts k¯ cs

k¯ st k¯ tt 0

⎫ ⎧ ⎫ ⎤⎧ k¯ sc ⎨ us ⎬ ⎨ r¯ s ⎬ 0 ⎦ ut = r¯ t . ⎩ ⎭ ⎩ ⎭ r¯ c uc k¯ cc

(42)

Note that the above modification of the discretized equation is employed only in the stiffness matrix. In the discretized equation (42), the equivalent nodal force vectors, which are derived from the potentials of both the penalty and the Nitsche’s methods, are calculated using the change in thickness and the transverse normal stress of the shell element.

123

If the equivalent nodal force vectors are computed in this approach with the change in thickness and the transverse normal stress ignored, the proposed approach may be regarded as the same procedure for connecting the MITC4 shell and the hexahedral elements. Hence, the discretized equation is modified with the minimal change of the stiffness matrix. Such a modification may cause numerical instabilities in computing the numerical solution of the discretized equation (42). In addition, the numerical solution obtained by (42) has not necessarily coincided with that by (41). Several numerical examples are presented to examine these fundamental features of the proposed approach. Static condensation can be employed in constructing the stiffness matrix on the connecting interface, and the degrees of freedom to represent thickness–stretch in the MITC4ts shell element are condensed out. From the second equation of (42), the displacement variations ut can be rewritten as   ut = k¯ tt−1 r¯ t − k¯ ts us .

(43)

Computational Mechanics

Substituting (43) into the first equation of (42), the discretized equation can be expressed as 

 k¯ ss ¯kcs

k¯ sc k¯ cc



us uc



=

 r¯ s , r¯ c

(44)

where  k¯ ss = k¯ ss − k¯ st k¯ tt−1 k¯ ts , r¯ s = r¯ s − k¯ st k¯ tt−1 r¯ t ,

which is the stiffness equation condensed with the degrees of freedom to represent thickness–stretch of the shell element on the interface. For the purpose of deriving the stiffness equation (44), it is necessary to treat ut independently of uc , which imposes the assumption on (41) to obtain (42). Note that the stiffness matrix in (44) may also be symmetry only when the Nitsche’s parameter α is equal to 0.5.

4 Numerical integration In this section, the numerical implementation to calculate the stiffness matrices and the nodal equivalent force vectors for connecting shell and solid elements is presented. In the discretized equations for connecting shell and hexahedral elements, the surface integration on the interface must be evaluated numerically. For simplicity in numerical integration, nodes of the shell element are located at the same position of those of the hexahedral element. Due to this restriction of the mesh density, the relationship of the positions of integration points in each domain becomes clear. In addition, the number of layers of the shell element has been taken as the same of elements in the transverse direction of the hexahedral element. Since the numerical integration on the connecting interface has not been restricted in constructing the stiffness equation, the non-matching meshes [39] or embedded meshes [36] can be used on the interface. In constructing the stiffness equations (31), (44) for connecting shell and hexahedral elements, the one-point integration for the penalty terms and 2 × 2 Gauss integration for the Nitsche’s terms are employed to evaluate the stiffness matrices and the nodal equivalent force vectors, as shown in Fig. 4. Hence, it is necessary to associate the positions of the integration pointat the local coordinates of the shell ele ment ξ (s) , η(s) , ζ (s) with those of the hexahedral element  (c) (c) (c)  ξ ,η ,ζ .

terms of the stiffness equations (31), (44) can be calculated using the reduced integration technique [47,48] as  f (s) f (c) d 0  (cs)     (s) (s) (s) (c) (c) (c) (cs) f (c) ξ0 , η0 , ζ0 w0 js(cs) , = f (s) ξˆ0 , ηˆ 0 , ζˆ0  (s) (s) (s)   (c) (c) (c)  where ξˆ0 , ηˆ 0 , ζˆ0 and ξ0 , η0 , ζ0 are the positions of the integration point for the reduced integration technique in the layer coordinate of the shell element and the local coordinate of the hexahedral element, respectively (see Fig. (cs) (cs) 4), w0 is the weight of the integration point, and js is the surface Jacobian. Note that the reduced integration technique should be employed for penalty terms to avoid traction oscillations [48,49] on the connecting interface.

4.2 Full integration technique for Nitsche’s terms In a similar manner to calculate the penalty terms of the stiffness equations (see Sect. 4.1), the Nitsche’s terms of the stiffness equations (31), (44) can be computed by employing the full integration technique as  f (s) f (c) d 0  (cs)

=

n gp 

    (s) (s) (s) (c) (c) (c) (cs) f (c) ξi , ηi , ζi wi js(cs) , f (s) ξˆi , ηˆ i , ζˆi

i=1

 (c) (c) (c)   (s) (s) (s)  where ξˆi , ηˆ i , ζˆi and ξi , ηi , ζi are the positions of the Gauss integration points for the full integration technique in the layer coordinate of the shell element and the local coordinate of the hexahedral element, respectively (see Fig. (cs) 4), wi is the weight of the integration point, and n gp is the number of the integration points.

5 Numerical examples Several numerical examples are chosen to verify the proposed procedure to couple shell to solid elements. Its performance is demonstrated in small and large deformation ranges by comparing the behavior of the proposed simulation model with that of the whole domain discretized with only solid elements. The material is described by the elastic material based on St. Venant–Kirchhoff model with Young’s modulus E and Poisson’s ratio ν. In the visualization of calculation results, we plot layers of shell elements as if they were hexahedral elements.

4.1 One-point integration technique for penalty terms

5.1 Tensile deformation of a bar

Let f (s) and f (c) be the functions (e.g. shape functions), which are defined in (s) and (c) , respectively. The penalty

To evaluate the behavior near the interface for connecting shell and solid elements, the simple model incorporating

123

Computational Mechanics

48

Midsurface 38

(a) Plane strain condition in y direction Solid Geometry:

Shell

Material:

48

Fig. 5 Geometry and material data of a bar for in-plane tensile deformation

38

Solid

in-plane tensile deformation is calculated. In addition, by changing the parameters, we observe the numerical stability of the proposed approach influenced by the model parameter. The geometric and material properties are given in Fig. 5. This model is discretized with 800 hexahedral elements (l (c) = 20 [mm], l (s) = 0 [mm]); 20 shell elements with 10 layers (l (c) = 0 [mm], l (s) = 20 [mm]); 40 hexahedral elements and 19 shell elements with 10 layers (l (c) = 1 [mm], l (s) = 19 [mm]). The stress distributions in the deformed configuration near the fixed edge is shown in Fig. 6. In the numerical results obtained by connecting the MITC4 shell and the hexahedral solid elements (see Fig. 6b), the unreasonable deformation, which caused by the discrepancy of the displacement in the thickness direction between shell and solid elements, can be observed near the connecting interface. On the other hand, the Nitsche’s method gives the reasonable deformation. Further, the reasonable deformation on the connecting interface can be obtained with the penalty parameter small enough by employing the Nitsche’s method. Thus, the influence of artificial numerical parameters on the numerical simulation results can be suppressed by using the proposed procedure. On the other hand, the numerical results obtained by connecting the MITC4ts shell and solid elements (see Fig. 6 c) illustrate that the continuities of the deformation in the thickness direction and the stress distribution are satisfied near the connecting interface. As seen in this figure, the reasonable stress distribution can be obtained from the three dimensional constitutive equations applied to the whole domain. Note that the solutions can be unreliable and inaccurate if the penalty parameter is set large in the connection of the MITC4ts shell and solid elements, due to the assumption for constructing a stiffness equation (see Sect. 3.3). However, this numerical instability can be avoided by using the Nitsche’s method with the small penalty parameter. Another numerical instability to connect shell and solid elements by using the Nitsche’s method can be observed

123

Shell 48

38

(b) Solid

Shell 48

38

Solid

Shell 48

38

(c) Fig. 6 Stress distributions of σx x of a bar for in-plane tensile deformation; a hexahedral solids, b connection of MITC4 shell and hexahedral solid, c connection of MITC4ts shell and hexahedral solid

in Fig. 7. The reasonable deformation can be obtained only when the Nitsche’s parameter is equal to zero. In α = 0.0, the slope of the interface is uniquely defined by using the director vector of the shell elements. On the other hand, in α = 1.0, the slope of the interface is expressed by using the position of nodes of solid elements near the connecting interface. In this situation, since the slope of the interface is

Computational Mechanics Table 1 Tip displacement of a cantilever beam for in-plane tensile deformation Discretization

l (c) (mm)

l (s) (mm)

Hexahedral solid 20 MITC4 shell and Hexa. solid MITC4ts shell and Hexa. solid

1

1

Parameter

0 19

19

x-dir. (mm) 5.1247

α = 0.0, p = 100 N/mm3

5.1247

p = 102 N/mm3

5.5224

α = 0.0, p = 100 N/mm3

5.1247

p = 102 N/mm3

5.5241

MITC4 shell

0

20

5.1340

MITC4ts shell

0

20

5.1340

Solid

[mm] (a)

Shell

(a) Solid

[mm] Shell

(b) Fig. 8 Norm of gap vectors on connecting surface of a bar for in-plane tensile deformation; a connection of MITC4 shell and hexahedral solid, b connection of MITC4ts shell and hexahedral solid

Solid

Shell 48

(b) Solid

Shell 38

Solid

Shell 48

(c) Solid

Shell 38

Fig. 9 Stress distributions of σx x of a bar for in-plane tensile deformation; connection of MITC4 shell and hexahedral solid

(d) Fig. 7 Deformation diagrams for connecting MITC4 shell and hexahedral solid of a bar for in-plane tensile deformation; a α = 0.5, p = 100 N/mm3 , b α = 1.0, p = 100 N/mm3 , c α = 0.5, p = 103 N/mm3 , d α = 1.0, p = 103 N/mm3

determined by the relative position of the nodes of solid elements, those nodes are not necessarily aligned on a straight line, and the shape of the interface is not uniquely defined, even if the penalty parameter is set large to enforce the con-

123

Computational Mechanics Table 2 Tip displacement of a cantilever beam with point load q = 4.0 × 10−3 N Discretization

l (c) (mm)

Hexahedral solid

20

0

1

19

MITC4 shell and Hexa. solid MITC4ts shell and Hexa. solid

1

l (s) (mm)

19

Parameter

x-dir. (mm)

z-dir. (mm)

− 3.9158 × 10−2

− 1.1418

α = 0.0, p = 101 N/mm3

− 3.9897 × 10−2

− 1.1511

p=

N/mm3

− 5.1574 × 10−2

− 1.3305

α = 0.0, p = 101 N/mm3

− 4.0253 × 10−2

− 1.1561

p = 102 N/mm3

− 5.1986 × 10−2

− 1.3357 − 1.1625 − 1.1683

102

MITC4 shell

0

20

− 4.0571 × 10−2

MITC4ts shell

0

20

− 4.0980 × 10−2

Table 3 Tip displacement of a cantilever beam with point load q = 0.2 N Discretization

l (c) (mm)

Hexahedral solid

20

MITC4 shell and Hexa. solid

1

MITC4ts shell and Hexa. solid

1

l (s) (mm)

Parameter

0 19 19

α = 0.0, p =

x-dir. (mm)

z-dir. (mm)

− 10.465

− 15.931

− 10.502

− 15.926

p = 102 N/mm3

− 11.402

− 16.647

α = 0.0, p = 101 N/mm3

− 10.543

− 15.951

101 N/mm3

− 11.412

− 16.650

MITC4 shell

0

20

− 10.557

− 15.983

MITC4ts shell

0

20

− 10.578

− 15.990

p=

Midsurface

Plane strain condition in y direction Geometry:

Material:

Fig. 10 Geometry and material data of a cantilever beam with point load

necting harder. Thus, the instability becomes remarkable as the Nitsche’s parameter approaches 1.0 (see Fig. 7). This feature is also observed when the MITC4ts shell element and hexahedral solid element are connected. The norm of gap vectors on the connecting interface is depicted in Fig. 8, where ζ is the thickness coordinate and ’+1’ and ’−1’ for ζ denote the position of top and bottom surfaces, respectively. Figure 8 shows that the norm decreases as the penalty parameter increase, when shell and solid elements are connected by only the penalty method. On the other hand, by the Nitsche’s method, the change in the norm is small when the penalty parameter is set small, and the norm is sufficiently

123

102

N/mm3

small. As the penalty parameter is set larger in employing the Nitsche’s method, the norm gets smaller, see Fig. 8. However, the deformation in the thickness direction and the stress distribution near the connecting interface seem to be unreasonable in Fig. 9, similar to the results obtained by only the penalty method (see Fig. 6 b). Thus, the penalty parameter is an important parameter that influences the determination of priority for the continuity condition of displacement vector or that of stress vector in the proposed approach. If the penalty parameter is set large, the continuity of the displacement vector is given priority over that of the stress vector. Note that a part of the numerical solution cannot be shown in Fig. 8 b since the numerical instabilities can be observed, due to the condensation of the stiffness matrices involving the deformation in the thickness direction (see Sect. 3.3). Considering these properties of the proposed approach, we verify the proposed procedure by employing the Nitsche’s method when the penalty parameter is set small in the subsequent numerical examples. We compare the displacement of the tip obtained by different discretizations to evaluate the effectiveness of the proposed procedure. The tip displacements of this model are shown in Table 1. The result obtained by the Nitsche’s method indicates a good performance, similar to solid elements. On the other hand, when discretizing with shell elements, since only the midsurface is constrained, an appropriate calculation result cannot be obtained.

Computational Mechanics

(a)

(a)

Solid

Shell

Solid

Shell

Solid

Shell

Solid

Shell

(b)

(b)

Solid

Shell

Solid

Shell

Solid

Shell

Solid

Shell

(c) Fig. 11 Stress distributions of σx x of a cantilever beam with point load q = 4.0 × 10−3 N; a Hexahedral solids, b connection of MITC4 shell and hexahedral solid, c connection of MITC4ts shell and hexahedral solid

(c) Fig. 12 Stress distributions of σzx of a cantilever beam with point load q = 4.0 × 10−3 N; a hexahedral solids, b connection of MITC4 shell and hexahedral solid, c connection of MITC4ts shell and hexahedral solid

5.2 Cantilever beam with point load A cantilever beam in plane strain subjected to a point load is chosen to demonstrate the effectiveness of the proposed procedure to couple shell to solid elements in small and large bending deformation. In this example, when uniform strain state does not occur on the connecting interface, and hence

we can verify that the shell and solid elements are properly connected by using the proposed approach in such situations. The geometry, material data, and boundary conditions are given in Fig. 10, and the model is discretized with 800 hexahedral elements (l (c) = 20 [mm], l (s) = 0 [mm]); 20 shell elements with 10 layers (l (c) = 0 [mm], l (s) = 20 [mm]);

123

Computational Mechanics

[mm] (a)

[mm] (b) Fig. 13 Norm of gap vectors on connecting surface of a cantilever beam with point load q = 4.0 × 10−3 N; a connection of MITC4 shell and hexahedral solid, b connection of MITC4ts shell and hexahedral solid

40 hexahedral elements and 19 shell elements with 10 layers (l (c) = 1 [mm], l (s) = 19 [mm]). Figures 11 and 12 illustrate the stress distributions near the fixed edge in small deformation range. The proposed procedure for connecting shell and solid elements gives the stress distributions of σx x similar to that in which whole domain is discretized with solid elements. On the other hand, Fig. 12 shows the discrepancy of the stress distributions of σzx near the connecting interface between the models using only the penalty method and that discretized with only solid elements. Thus, this discrepancy affects the behavior of the whole structure, as shown in Table 2. These results confirm that the proposed approach employing the Nitsche’s method is capable of connecting shell and solid elements reasonably. The norm of gap vectors on the connecting interface is shown in Fig. 13. The results obtained using the Nitsche’s method show that the norm is sufficiently small if the penalty parameter was changed. Figures 14 and 15 indicate the stress distributions in large deformation range. The stress distributions obtained by the Nitsche’s method are in good agreement with those of the model discretized with only solid elements. In addition, the

123

norm of gap vectors on the connecting interface becomes small enough by using the Nitsche’s method, as illustrated in Fig. 16. Furthermore, the Nitsche’s method improves the tip displacements of this model, similar to the model discretized with only solid elements, as shown in Table 3. Thus, these results illustrate that the proposed procedure using the Nitsche’s method can be more efficient than that employing only the penalty method. The cantilever beam or plate can be a suitable example for the verification of the combined modeling. In the conventional procedure of introducing transition elements[16], it is important to assess the behavior of the models by changing the size of domains discretized with shell, solid, and transition elements. On the other hand, when using the proposed approach, it is necessary to verify only the size of domains discretized with shell and solid elements. Note that it is also inevitable to assess the numerical parameters set on the connecting interface. This task is similar to select the constitutive equations of the region discretized with transition elements for adjusting the stress state between shell and solid elements[5]. Thus, no transition elements on the interface is required in the proposed approach, and hence the simulation models can be clearly distinguished between the regions discretized with shell and solid elements.

5.3 Cantilever beam with surface traction We examine a cantilever beam in plane strain applied a distributed load on the top surface, as shown in Fig. 17. This model is the most attractive example to demonstrate the predictive capabilities of the proposed procedure by comparing the stress distributions near the connecting interface. The geometric and material properties are also given in Fig. 17. This model is discretized with 800 hexahedral elements (l (c) = 20 [mm], l (s) = 0 [mm]); 20 shell elements with 10 layers (l (c) = 0 [mm], l (s) = 20 [mm]); 40 hexahedral elements and 19 shell elements with 10 layers (l (c) = 1 [mm], l (s) = 19 [mm]). Figures 18, 19, and 20 illustrate the stress distributions near the fixed edge in small deformation range. The proposed approach provides sufficiently accurate stress distributions of σx x and σzx , similar to that in which whole domain is discretized with the hexahedral solid elements. On the other hand, Fig. 20 shows the discontinuities of the stress distributions of σzz on the interface connecting the MITC4 shell element and the hexahedral solid element. Within the formulation of the MITC4 shell element, the surface traction is evaluated on the midsurface, because the plane stress assumption is enforced in the thickness direction. Thus, it is inevitable that the discontinuities of the stress distributions on the interface connecting the MITC4 shell and solid elements in this situation. It should be noted that this feature can be observed for the connection of conventional shell

Computational Mechanics

(a)

Solid

(a)

Solid Shell

Shell

Solid

Solid Shell

Shell (b)

Solid

(b)

Solid Shell

Shell

Solid

Solid Shell

Shell (c)

(c)

Fig. 14 Stress distributions of σx x of a cantilever beam with point load q = 0.2N; a Hexahedral solids, b connection of MITC4 shell and hexahedral solid, c connection of MITC4ts shell and hexahedral solid

Fig. 15 Stress distributions of σzx of a cantilever beam with point load q = 0.2 N; a hexahedral solids, b connection of MITC4 shell and hexahedral solid, c connection of MITC4ts shell and hexahedral solid

123

Computational Mechanics

(a) Solid

Shell

Solid

Shell

[mm] (a)

(b) [mm]

Solid

Shell

Solid

Shell

(b) Fig. 16 Norm of gap vectors on connecting surface of a cantilever beam with point load q = 0.2 N; a connection of MITC4 shell and hexahedral solid, b connection of MITC4ts shell and hexahedral solid

Midsurface

Plane strain condition in y direction Geometry:

Material:

(c)

Fig. 17 Geometry and material data of a cantilever beam with distributed load

Fig. 18 Stress distributions of σx x of a cantilever beam with distributed load q = 2.0 × 10−4 N/mm2 ; a hexahedral solids, b connection of MITC4 shell and hexahedral solid, c connection of MITC4ts shell and hexahedral solid

and solid elements. On the other hand, the MITC4ts shell element is capable of capturing the surface traction at the top surface. Hence, the reasonable stress distributions of σzz near the interface can be obtained when the MITC4ts shell and solid elements are connected, as shown in Fig. 20. Figures 21, 22, and 23 indicate the stress distributions in large deformation range. The stress distributions of σx x and

σzx obtained by the Nitsche’s method are in good agreement with those calculated with the hexahedral solid elements. In addition, Fig. 23 shows that the proposed procedure for connecting the MITC4ts shell and solid elements provides sufficiently accurate stress distributions of σzz , similar to the model discretized with only solid elements. Thus, the approach connecting the MITC4ts shell and solid elements

123

Computational Mechanics

(a)

(a)

Solid

Shell

Solid

Shell

Solid

Shell

Solid

Shell

(b)

(b)

Solid

Shell

Solid

Shell

Solid

Shell

Solid

Shell

(c)

(c)

Fig. 19 Stress distributions of σzx of a cantilever beam with distributed load q = 2.0 × 10−4 N/mm2 ; a hexahedral solids, b connection of MITC4 shell and hexahedral solid, c connection of MITC4ts shell and hexahedral solid

Fig. 20 Stress distributions of σzz of a cantilever beam with distributed load q = 2.0 × 10−4 N/mm2 ; a hexahedral solids, b connection of MITC4 shell and hexahedral solid, c connection of MITC4ts shell and hexahedral solid

can be more efficient than that coupling the MITC4 shell to solid elements in this situation. The tip displacements of this model are depicted in Tables 4 and 5. These results show that the proposed procedure is effective in connecting shell and solid elements.

5.4 Tensile deformation of a stepped bar A stepped bar in plane strain incorporating in-plane tensile deformation is investigated, as shown in Fig. 24. In this model, the geometric singularity can be observed at the portion where the thickness is discontinuous. Hence, we verify

123

Computational Mechanics

(a)

(a)

Solid

Solid

Shell

Shell

Solid

Solid

Shell

Shell

(b)

(b)

Solid

Solid

Shell

Shell

Solid

Solid

Shell

Shell (c)

Fig. 21 Stress distributions of σx x of a cantilever beam with distributed load q = 1.0 × 10−2 N/mm2 ; a hexahedral solids, b connection of MITC4 shell and hexahedral solid, c connection of MITC4ts shell and hexahedral solid

123

(c)

Fig. 22 Stress distributions of σzx of a cantilever beam with distributed load q = 1.0 × 10−2 N/mm2 ; a hexahedral solids, b connection of MITC4 shell and hexahedral solid, c connection of MITC4ts shell and hexahedral solid

Computational Mechanics

Midsurface

(a)

Plane strain condition in y direction Material:

Geometry:

Fig. 24 Geometry and material data of a stepped bar for in-plane tensile deformation

Solid Shell

Solid

(a)

Shell

(b)

Solid

(b)

Shell

Fig. 25 Stress distributions of σx x of a stepped bar discretized with one kind of element for in-plane tensile deformation; a hexahedral solid element, b MITC4ts shell element

Solid Shell (c)

Fig. 23 Stress distributions of σzz of a cantilever beam with distributed load q = 1.0 × 10−2 N/mm2 ; a hexahedral solids, b connection of MITC4 shell and hexahedral solid, c connection of MITC4ts shell and hexahedral solid

the effectiveness of the proposed approach by changing the size of domains discretized shell and solid elements. The geometric and material properties are also depicted in Fig. 24, and the model is discretized with 2, 200 hexahedral elements (l (c) = 20 [mm], l (s) = 0 [mm]); 20 shell elements with 10 layers (l (c) = 0 [mm], l (s) = 20 [mm]); 300 hexahedral elements and 19 shell elements with 10 layers (l (c) = 1 [mm], l (s) = 19 [mm]); 400 hexahedral elements and 18 shell elements with 10 layers (l (c) = 2 [mm], l (s) = 18 [mm]). In this example, it is not possible to assess the local behavior at the corner using only shell elements, as shown in Fig. 25. In contrast, the proposed procedure can be an appropriate approach of evaluating the local behaviors.

123

Computational Mechanics Table 4 Tip displacement of a cantilever beam with distributed load q = 2.0 × 10−4 N/mm2 Discretization

l (c) (mm)

Hexahedral solid

20

0

1

19

MITC4 shell and Hexa. solid MITC4ts shell and Hexa. solid

1

l (s) (mm)

19

Parameter

x-dir. (mm)

z-dir. (mm)

− 5.2803 × 10−3

− 0.43034

α = 0.0, p = 101 N/mm3

− 5.3477 × 10−3

− 0.43223

p=

N/mm3

− 7.4858 × 10−3

− 0.51827

α = 0.0, p = 101 N/mm3

− 5.3981 × 10−3

− 0.43450

p = 102 N/mm3

− 7.5379 × 10−3

− 0.52028

× 10−3

− 0.43788 − 0.44017

102

MITC4 shell

0

20

− 5.4730

MITC4ts shell

0

20

− 5.5242 × 10−3

Table 5 Tip displacement of a cantilever beam with distributed load q = 1.0 × 10−2 N/mm2 Discretization

l (c) (mm)

Hexahedral solid

20

MITC4 shell and Hexa. solid

1

MITC4ts shell and Hexa. solid

1

l (s) (mm)

Parameter

0 19 19

α = 0.0, p =

x-dir. (mm)

z-dir. (mm)

− 6.4411

− 13.632

− 6.0118

− 13.275

p = 102 N/mm3

− 7.1531

− 14.429

α = 0.0, p = 101 N/mm3

− 6.0615

− 13.321

101 N/mm3

− 7.1758

− 14.445

MITC4 shell

0

20

− 6.0837

− 13.362

MITC4ts shell

0

20

− 6.1165

− 13.388

p=

Solid

Shell

(a)

Solid

102

N/mm3

Solid

Shell

(a)

Shell

Solid

Shell

(b)

(b)

Fig. 26 Stress distributions of σx x of a stepped bar discretized with hexahedral solid (l (c) = 1 mm) and MITC4ts shell (l (s) = 19 mm) for in-plane tensile deformation; a α = 0.0, p = 101 N/mm3 , b p = 102 N/mm3

Fig. 27 Stress distributions of σx x of a stepped bar discretized with hexahedral solid (l (c) = 2 mm) and MITC4ts shell (l (s) = 18 mm) for in-plane tensile deformation; a α = 0.0, p = 101 N/mm3 , b p = 102 N/mm3

Figures 26 and 27 indicate the stress distributions of σx x in the deformed configuration. It is not possible to capture the

stress concentration when the corner is discretized with shell elements (see Fig. 26). On the other hand, when the corner is

123

Computational Mechanics Table 6 Tip displacement of a stepped beam for in-plane tensile deformation Discretization

l (c) (mm)

Hexahedral solid

20

0

MITC4ts shell and Hexa. solid

1

19

2

MITC4ts shell

0

l (s) (mm)

Parameter

x-dir. (mm)

Solid Shell

5.9594

18

α = 0.0, p = 100 N/mm3

5.9726

p = 102 N/mm3

6.4639

α = 0.0, p = 100 N/mm3

5.9593

p = 102 N/mm3

6.4585

20

(a)

5.8871

Solid Shell

(b)

Midsurface

Fig. 30 Stress distributions of σzx of a stepped beam with point load discretized with one kind of element; a hexahedral solid elements, b MITC4ts shell element

Plane strain condition in y direction Geometry:

Material: Solid Shell

Fig. 28 Geometry and material data of a stepped beam with point load

(a)

Solid Shell

(a)

(b) Fig. 31 Stress distributions of σx x of a stepped beam with point load discretized with hexahedral solid (l (c) = 1 mm) and MITC4ts shell (l (s) = 19 mm); a α = 0.0, p = 101 N/mm3 , b p = 102 N/mm3

(b) Fig. 29 Stress distributions of σx x of a stepped beam with point load discretized with one kind of element; a hexahedral solid element, b MITC4ts shell element

modeled with solid elements (see Fig. 27), the stress concentration at the corner can be represented appropriately, and the reasonable deformation can be obtained using the Nitsche’s method. In the proposed procedure, it should be noted that the portion, which the local behavior should be assessed, is discretized with solid elements. This important feature can be confirmed the numerical results of the tip displacements, as illustrated in Table 6.

123

Computational Mechanics Table 7 Tip displacement of a stepped beam with point load q = 0.2 N Discretization

l (c) l (s) Parameter (mm) (mm)

x-dir. (mm)

Hexahedral solid

20

0

− 9.9140 − 15.293

MITC4ts shell 1 and Hexa. solid

19

− 10.015 − 15.375 α = 0.0, p = 101 N/mm3 p = 102 N/mm3

2

18

0

Solid Shell

− 10.917 − 16.187

(a)

− 9.9095 − 15.275 α = 0.0, p = 101 N/mm3 p = 102 N/mm3

MITC4ts shell

z-dir. (mm)

− 10.736 − 15.903 − 10.024 − 15.194

20

Solid Shell

(b) Fig. 33 Stress distributions of σx x of a stepped beam with point load discretized with hexahedral solid (l (c) = 2 mm) and MITC4ts shell (l (s) = 18 mm); a α = 0.0, p = 101 N/mm3 , b p = 102 N/mm3

(a)

Solid Shell

(a) (b) Fig. 32 Stress distributions of σzx of a stepped beam with point load discretized with hexahedral solid (l (c) = 1 mm) and MITC4ts shell (l (s) = 19 mm); a α = 0.0, p = 101 N/mm3 , b p = 102 N/mm3 .

Solid Shell

The structures, in which the thickness is discontinuous, can be modeled by using the conventional procedure of introducing transition elements[8]. However, it is necessary to set the size of domains discretized with shell, solid, and transition elements suitably for obtaining reasonable numerical results. In addition, the constitutive equation of the region discretized with transition elements should be determined appropriately. On the other hand, in the proposed approach of connecting the MITC4ts shell and solid elements, the same constitutive equation can be applied uniformly to the whole domain. Thus, the proposed procedure is effective in modeling the complex structures, in which the thickness is not uniform.

123

(b) Fig. 34 Stress distributions of σzx of a stepped beam with point load discretized with hexahedral solid (l (c) = 2 mm) and MITC4ts shell (l (s) = 18 mm); a α = 0.0, p = 101 N/mm3 , b p = 102 N/mm3

5.5 Bending of a stepped beam A stepped beam in plane strain subjected to a point load is employed to examine the behavior of the model discretized with the proposed procedure, see Fig. 28. In this example, we verify the behavior in bending deformation of the structure including the geometric singularity to evaluate the effec-

Computational Mechanics

(a)

Midsurface Geometry: Material:

Fig. 35 Geometry and material data of a plate with distributed load

Solid Shell

tiveness of the numerical procedure to couple shell to solid elements. The geometry and material data for the numerical simulation are also given in Fig. 28. This model is discretized with 2, 200 hexahedral elements (l (c) = 20 [mm], l (s) = 0 [mm]); 20 shell elements with 10 layers (l (c) = 0 [mm], l (s) = 20 [mm]); 300 hexahedral elements and 19 shell elements with 10 layers (l (c) = 1 [mm], l (s) = 19 [mm]); 400 hexahedral elements and 18 shell elements with 10 layers (l (c) = 2 [mm], l (s) = 18 [mm]). In this example, it is not possible to evaluate the local behavior using only shell elements, as shown in Figs. 29 and 30. In the proposed approach, when the corner is discretized with shell elements, it is not sufficient for assessing the stress concentration near the corner, as illustrated in Figs. 31 and 32. On the other hand, Figs. 33 and 34 show that the local behaviors can be resolved properly, when the corner is discretized with solid elements. Further, the validity of the proposed procedure can also be confirmed from the results of the tip displacements, as shown in Table 7.

Solid Shell

(b)

Solid Shell

5.6 Plate bending with surface traction We examine a plate applied a distributed load on the top surface, as shown in Fig. 35. This example is chosen to demonstrate the predictive capabilities of the proposed procedure in three dimensional analysis. The geometry and material data for the numerical simulation are also given in Fig. 35. This model is discretized with 4, 000 hexahedral elements (l (c) = 20 [mm], l (s) = 0 [mm]); 200 shell elements with 10 layers (l (c) = 20 [mm], l (s) = 20 [mm]); 2, 000 hexahedral elements and 100 shell elements with 10 layers (l (c) = 10 [mm], l (s) = 10 [mm]). Figures 36, 37, and 38 show the stress distributions of the plate in the deformed configuration. The non-physical oscillation in stress distributions near the connecting interface in Figs. 36, 37, and 38 come from the discrepancy of

Solid Shell

(c) Fig. 36 Stress distributions of σx x of a plate with distributed load; a hexahedral solids, b connection of MITC4 shell and hexahedral solid, c connection of MITC4ts shell and hexahedral solid

123

Computational Mechanics

(a)

(a)

Solid

Solid

Shell

Shell

Solid

Solid

Shell

Shell

(b)

(b)

Solid

Solid

Shell

Shell

Solid

Solid Shell

Shell

(c) Fig. 37 Stress distributions of σzx of a plate with distributed load; a hexahedral solids, b connection of MITC4 shell and hexahedral solid, c connection of MITC4ts shell and hexahedral solid

123

(c) Fig. 38 Stress distributions of σzz of a plate with distributed load; a hexahedral solids, b connection of MITC4 shell and hexahedral solid, c connection of MITC4ts shell and hexahedral solid

Computational Mechanics Table 8 Tip displacement of a plate at center with distributed load Discretization

l (c) (mm)

l (s) (mm)

Hexahedral solid MITC4 shell and Hexa. solid

20 10

0 10

MITC4ts shell and Hexa. solid

10

10

0 0

20 20

MITC4 shell MITC4ts shell

Parameter

α = 0.0, p = 100 N/mm3 p = 102 N/mm3 α = 0.0, p = 100 N/mm3 p = 102 N/mm3

x-dir. (mm)

z-dir. (mm)

− 0.48820 − 0.46285 − 0.51460 − 0.46897 − 0.51474 − 0.56795 − 0.57186

− 4.0998 − 3.9989 − 4.1889 − 4.0215 − 4.1895 − 4.4167 − 4.4314

x-dir. (mm)

z-dir. (mm)

− 0.48684 − 0.46095 − 0.51218 − 0.46648 − 0.51233 − 0.56512 − 0.56904

− 4.0947 − 3.9918 − 4.1804 − 4.0122 − 4.1811 − 4.4072 − 4.4219

Table 9 Tip displacement of a plate at corner with distributed load Discretization

l (c) (mm)

l (s) (mm)

Hexahedral solid MITC4 shell and Hexa. solid

20 10

0 10

MITC4ts shell and Hexa. solid

10

10

0 0

20 20

MITC4 shell MITC4ts shell

the transverse shear components of stress. In the domain discretized with shell elements, the transverse shear components of stress are assumed to be constant with respect to the transverse direction in the Reissner–Mindlin plate theory. Thus, it is difficult to alleviate the non-physical oscillation in stress distributions in the proposed approach. However, the proposed approach employing the Nitsche’s method accurately predict the stress distributions of the whole structure except for the connecting interface, similar to those of the whole domain is discretized with solid elements, see Figs. 37 and 38. Furthermore, the connection of the MITC4ts shell and solid elements with the Nitsche’s method provides a good performance to capture the tip displacements at the center and the corner, similar to the model discretized with only solid elements, as shown in Tables 8 and 9.

5.7 Inflation of a cylinder subjected to internal pressure A simulation model of a cylinder is chosen to demonstrate the performance of the proposed procedure in the analysis of a thick-walled curved structure. The geometric and material properties are given in Fig. 39. On account of symmetry, one eighth of the cylinder is modeled using 12, 500 hexahedral elements (l (c) = 25 [mm], l (s) = 0 [mm]); 1, 250 shell elements with 10 layers (l (c) = 0 [mm], l (s) = 25 [mm]); 1, 500 hexahedral elements and 100 shell elements with 10 layers (l (c) = 3 [mm], l (s) = 22 [mm]). Figure 40 indicates the von Mises stress distributions of the cylinder in the deformed configuration. As seen in this

Parameter

α = 0.0, p = 100 N/mm3 p = 102 N/mm3 α = 0.0, p = 100 N/mm3 p = 102 N/mm3

Geometry: Material:

Fig. 39 Geometry and material data of a cylinder

figure, the MITC4ts shell element has the capability to assess the pressure that is applied to the internal surface. However, when the whole domain is discretized with shell elements, a boundary condition at the clamped edge cannot be considered adequately. Hence, the difference between the calculated results of the MITC4ts shell element and the hexahedral solid element becomes larger near the clamped edge. On the other hand, the connection of the MITC4ts shell and solid elements gives the stress distributions similar to that in which whole domain is discretized with solid elements. Thus, the proposed procedure for connecting the MITC4ts shell and solid elements is more effective than the simulation model discretized with only shell elements in this situation. It should be noted that the numerical results obtained by the proposed approach indicate the non-physical oscillation in stress distribution near the connecting interface due to the discrepancy

123

Computational Mechanics

: cross section

(a)

(b) Shell

Solid

Shell

Solid

(c) Fig. 40 von Mises stress distributions of a cylinder; a hexahedral solids, b MITC4ts shells, c connection of MITC4ts shell and hexahedral solid

of the transverse shear components of stress, as shown in Fig. 40.

6 Concluding remarks A numerical procedure to couple shell to solid elements by using the Nitsche’s method has been developed. In the proposed approach, the domain of the problem of interest is decomposed into two domains discretized with shell and solid elements. On the connecting interface between shell and solid elements, the continuity of the displacement vector is satisfied approximately by the penalty method, and that of

123

the stress vector is imposed by the Nitsche’s method. Two types of shell elements coupled with solid elements have been presented in this paper. One of them is the conventional MITC4 shell element, and the other is the modified MITC4 shell element with additional degrees of freedom to represent thickness–stretch. The connection of the MITC4 shell and solid elements by only the penalty method shows the unreasonable deformation on the connecting interface, due to the discrepancy of the displacement in the thickness direction between shell and solid elements. In this study, by virtue of considering the continuity of the stress vector on the connecting interface, the continuity of deformation can be improved. In addition, the reasonable deformation can be obtained with the penalty parameter small enough by using the Nitsche’s method. However, the unreasonable stress distribution can be observed near the connecting interface when the surface traction is subjected, because the plane stress condition is assumed in the transverse direction for conventional shell elements. In this approach, the surface traction is evaluated on the surface where the traction is applied in the domain discretized with solid elements, whereas it is assessed on the midsurface in the domain discretized with the MITC4 shell element. Hence, it is difficult to satisfy the continuity of the stress vector on the connecting interface by using the coupling of the MITC4 shell and solid elements in this situation. On the other hand, the coupling of the MITC4ts shell and solid elements has the capability to represent the stress state in a full three dimension for the whole domain. The most attractive property of this approach is that the MITC4ts shell element has a capability of employing the three dimensional constitutive equations without any assumptions. Hence, the same constitutive equation can be applied uniformly to the whole domain, which is discretized with both shell and solid elements. Thus, the continuity of the stress vectors including the transverse normal component can be satisfied on the interface. It should be noted that the numerical results indicate the discrepancy of the transverse shear components of stress, which are assumed to be constant in the domain discretized with shell elements with respect to the transverse direction in the Reissner–Mindlin plate theory. In addition, due to the formulation of the MITC4ts shell element, the discontinuities of the deformation in the thickness direction can be allowed at each element and on the connecting interface. Thus, static condensation can be employed in constructing the stiffness matrix on the connecting interface. Note that this relaxation causes numerical instabilities in computing the numerical solution. No transition elements on the interface in the connection of shell and solid elements is required in the proposed approach, and hence the simulation models can be clearly distinguished between the regions discretized with different types of elements. On the other hand, the conventional pro-

Computational Mechanics

cedure of introducing transition elements[15,16] should be verified the domain modeled by transition elements in addition to the regions discretized with shell and solid elements. In order to assess the validity of the conventional procedure, it is inevitable to repeat the pre-processing for mesh generation to change the size of the regions discretized with shell, solid, and transition elements. It is also noted that the numerical parameter, which is defined the stress state in the transition elements, should be determined appropriately[5]. Therefore, the proposed approach can be more efficient for coupling shell to solid elements than the conventional procedures. Representative numerical examples demonstrate the applicability and validity of the proposed procedure to couple shell to solid elements. The results reveal that the solutions obtained by the proposed approach provides more accurate than those by the approach using only the penalty method. In future work, we will apply the proposed procedure to practical problems in more complex structures. In such situations, the effectiveness of the proposed procedure and the validity of decomposed domains should be verified in detail.

References 1. Surana KS (1979) Isoparametric elements for cross-sectional properties and stress analysis of beams. Int J Numer Methods Eng 14:475–497 2. Surana KS (1980) Transition finite elements for axisymmetric stress analysis. Int J Numer Methods Eng 15:809–832 3. Surana KS (1980) Transition finite elements for three-dimensional stress analysis. Int J Numer Methods Eng 15:991–1020 4. Bathe KJ, Ho LW (1981) Some results in the analysis of thin shell structures. In: Wunderlich W, Stein E, Bathe KJ (eds) Nonlinear finite element analysis in structural mechanics. Springer, Berlin, pp 122–150 5. Surana KS (1982) Geometrically nonlinear formulation for the three dimensional solid-shell transition finite elements. Comput Struct 15:549–566 6. Surana KS (1983) Geometrically non-linear formulation for the axi-symmetric transition finite elements. Comput Struct 17:243– 255 7. Cofer WF, Will KM (1991) A three-dimensional, shell-solid transition element for general nonlinear analysis. Comput Struct 38:449–462 8. Gmür TC, Schorderet AM (1993) A set of three-dimensional solid to shell transition elements for structural dynamics. Comput Struct 46:583–591 9. Chavan KS, Wriggers P (2004) Consistent coupling of beam and shell models for thermo-elastic analysis. Int J Numer Methods Eng 59:1861–1878 10. Gmür TC, Kauten RH (1993) Three-dimensional solid to beam transition elements for structural dynamics analysis. Int J Numer Methods Eng 36:1429–1444 11. Dohrmann CR, Key SW (1999) A transition element for uniform strain hexahedral and tetrahedral finite elements. Int J Numer Methods Eng 44:1933–1950 12. Dohrmann CR, Key SW, Heinstein MW (2000) Methods for connecting dissimilar three-dimensional finite element meshes. Int J Numer Methods Eng 47:1057–1080

13. Garusi E, Tralli A (2002) A hybrid stress-assumed transition element for solid-to-beam and plate-to-beam connections. Comput Struct 80:105–115 14. Ahn JS, Basu PK, Woo KS (2010) Analysis of cracked aluminum plates with one-sided patch repair using p-convergent layered model. Finite Elem Anal Des 46:438–448 15. Ahn JS, Basu PK (2011) Locally refined p-FEM modeling of patch repaired plates. Compos Struct 93:1704–1716 16. Ahn JS, Woo KS (2011) Analysis of Cantilever plates with stepped section using p-convergent transition element for solid-to-shell connections. Adv Struct Eng 14:1167–1183 17. Liu GR, Quek SS (2013) The finite element method: a practical course, 2nd edn. Butterworth-Heinemann, Oxford 18. Curiskis JI, Valliappan S (1978) A solution algorithm for linear constraint equations in finite element analysis. Comput Struct 8:117–124 19. Abel JF, Shephard MS (1979) An algorithm for multipoint constraints in finite element analysis. Int J Numer Methods Eng 14:464–467 20. Shephard MS (1984) Linear multipoint constraints applied via transformation as part of a direct stiffness assembly process. Int J Numer Methods Eng 20:2107–2112 21. Bournival S, Cuilliere JC, Francois V (2010) A meshgeometry based method for coupling 1D and 3D elements. Adv Eng Softw 41:838–858 22. Jialin W, Zhiwen L, Xing M, Jianzheng Z (1996) A DOF expanding method for connecting solid and shell element. Commun Numer Methods Eng 12:321–330 23. Bathe KJ (1996) Finite element procedure. Prentice-Hall Inc, Upper Saddle River 24. Monaghan DJ, Doherty IW, Court DM, Armstrong CG (1998) Coupling 1D beams to 3D bodies. In: 7th International meshing roundtable. Sandia National Laboratories, Dearborn, Michigan, pp 285–293 25. McCune RW, Armstrong CG, Robinson DJ (2000) Mixeddimensional coupling in finite element models. Int J Numer Methods Eng 49:725–750 26. Shim KW, Monaghan DJ, Armstrong CG (2002) Mixed dimensional coupling in finite element stress analysis. Eng Comput 18:241–252 27. Robinson TT, Armstrong CG, Fairey R (2011) Automated mixed dimensional modelling from 2D and 3D cad models. Finite Elem Anal Des 47:151–165 28. Aminpour MA, Krishnamurthy T (1997) A two-dimensional interface element for multi-domain analysis of independently modeled three-dimensional finite element meshes. AIAA, No.97-1297, 1853–1861 29. Aminpour MA, Krishnamurthy T, Fadale TD (1998) Coupling of independently modeled three-dimensional finite element meshes with arbitrary shape interface boundaries. AIAA, No.98-2060, 3014–3024 30. Osawa N, Hashimoto K, Sawamura J, Nakai T, Suzuki S (2007) Study on shell-solid coupling FE analysis for fatigue assessment of ship structure. Mar Struct 20:143–163 31. Mata P, Barbat AH, Oller S (2008) Two-scale approach for the nonlinear dynamic analysis of RC structures with local non-prismatic parts. Eng Struct 30:3667–3680 32. Nitsche J (1971) Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind. Abhandlungen aus dem Mathematischen Seminar der Universitat Hamburg 36:9– 15 33. Hansbo A, Hansbo P (2002) An unfitted finite element method, based on Nitsche’s method, for elliptic interface problems. Comput Methods Appl Mech Eng 191:5537–5552

123

Computational Mechanics 34. Hansbo P (2005) Nitsche’s method for interface problems in computational mechanics. GAMM-Mitt 28:183–206 35. Sanders JD, Puso MA (2012) An embedded mesh method for treating overlapping finite element meshes. Int J Numer Methods Eng 91:289–305 36. Sanders JD, Laursen TA, Puso MA (2012) A Nitsche embedded mesh method. Comput Mech 49:243–257 37. Fritz A, Hüeber S, Wohlmuth BI (2004) A comparison of mortar and Nitsche techniques. CALCOLO 41:115–137 38. Apostolatos A, Schmidt R, Wüchner R, Bletzinger KU (2014) A Nitsche-type formulation and comparison of the most common domain decomposition methods in isogeometric analysis. Int J Numer Methods Eng 97:473–504 39. Nguyen VP, Kerfriden P, Brino M, Bordas SPA, Bonisoli E (2014) Nitsche’s method for two and three dimensional NURBS patch coupling. Comput Mech 53:1163–1182 40. Costa JC, Pimenta PM, Wriggers P (2016) Meshless analysis of shear deformable shells: boundary and interface constraints. Comput Mech 57:679–700 41. Arnold DN, Brezzi F, Cockburn B, Marini LD (2002) Unified analysis of discontinuous galerkin methods for elliptic problems. SIAM J Numer Anal 39:1749–1779 42. Mergheim J, Steinmann P (2006) A geometrically nonlinear FE approach for the simulation of strong and weak discontinuities. Comput Methods Appl Mech Eng 195:5037–5052

123

43. Wriggers P, Zavarise G (2008) A formulation for frictionless contact problems using a weak form introduced by Nitsche. Comput Mech 41:407–420 44. Guo Y, Ruess M (2015) Nitsche’s method for a coupling of isogeometric thin shells and blended shell structures. Comput Methods Appl Mech Eng 284:881–905 45. Dvorkin EN, Bathe KJ (1984) A continuum mechanics based four node shell element for general non-linear analysis. Eng Comput 1:77–88 46. Yamamoto T, Yamada T, Matsui K (2017) A quadrilateral shell element with degree of freedom to represent thickness–stretch. Comput Mech 59:625–646 47. Zienkiewicz OC, Taylor RL, Too JM (1971) Reduced integration techniques in general analysis of plates and shells. Int J Numer Methods Eng 3:275–290 48. Kikuchi N (1982) A smoothing technique for reduced integration penalty methods in contact problems. Int J Numer Methods Eng 18:343–350 49. Oden JT, Kikuchi N (1982) Finite element methods for constrained problems in elasticity. Int J Numer Methods Eng 18:701–725

Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.