Concurrently coupled solid shell-based adaptive multiscale method for fracture

Available online at www.sciencedirect.com ScienceDirect Comput. Methods Appl. Mech. Engrg. 319 (2017) 338–365 www.elsevier.com/locate/cma Concurrently coupled solid shell-based adaptive multiscale method for fracture P.R. Budarapua, J. Reinosob,∗, M. Paggia aMulti-scale Analysis of Materials Research Unit, IMT School for Advanced Studies Lucca, Piazza San Francesco 19, 55100 Lucca, Italy bElasticity and Strength of Materials Group, School of Engineering, University of Seville, Camino de los Descubrimientos s/n, 41092, Seville, Spain Received 18 October 2016; received in revised form 23 January 2017; accepted 16 February 2017 Available online 1 March 2017 Highlights •Continuum-based phantom node method combined with an enhanced strain-based solid shell element. •Coupling with a molecular statics to generate a multiscale framework for the simulation of cracks in thin structures. •Coupling of the continuum and atomistic models performed through the use of a bridging scale method. •Simulations of crack propagation in three scenarios for Silicon specimens with a cubic diamond lattice structure containing propagating cracks. Abstract A solid shell-based adaptive atomistic–continuum numerical method is herein proposed to simulate complex crack growth patterns in thin-walled structures. A hybrid solid shell formulation relying on the combined use of the enhanced assumed strain (EAS) and the assumed natural strain (ANS) methods has been considered to efficiently model the material in thin structures at the continuum level. The phantom node method (PNM) is employed to model the discontinuities in the bulk. The discontinuous solid shell element is then concurrently coupled with a molecular statics model placed around the crack tip. The coupling between the coarse scale and the fine scale is realized through the use of ghost atoms, whose positions are interpolated from the coarse scale solution and enforced as boundary conditions to the fine scale model. In the proposed numerical scheme, the fine scale region is adaptively enlarged as the crack propagates and the region behind the crack tip is adaptively coarsened in order to reduce the computation costs. An energy criterion is used to detect the crack tip location. All the atomistic simulations are carried out using the LAMMPS software. A computational framework has been developed in MATLAB to trigger LAMMPS through system command. This allows a two way interaction between the coarse and fine scales in MATLAB platform, where the boundary conditions to the fine region are extracted from the coarse scale, and the crack tip location from the atomistic model is transferred back to the continuum scale. The developed framework has been applied to study crack growth in the energy minimization problems. Inspired by the influence of fracture on current–voltage characteristics of thin Silicon photovoltaic cells, the cubic diamond lattice structure of Silicon is used to model the material in the fine scale region, whilst the Tersoff potential function is employed to model the ∗Corresponding author. E-mail addresses: [email protected] (P.R. Budarapu), [email protected] (J. Reinoso), [email protected] (M. Paggi). http://dx.doi.org/10.1016/j.cma.2017.02.023 0045-7825/ c 2017 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http:// creativecommons.org/licenses/by-nc-nd/4.0/). P.R. Budarapu et al. / Comput. Methods Appl. Mech. Engrg. 319 (2017) 338–365 339 atom–atom interactions. The versatility and robustness of the proposed methodology is demonstrated by means of several fracture applications. c 2017 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http:// creativecommons.org/licenses/by-nc-nd/4.0/). Keywords: Multiscale methods; Solid shell finite element; Phantom node method for fracture; Atomistic simulations; Adaptivity; Silicon solar cells 1. Introduction In engineering applications, the global response of the system is often governed by the material behavior at small length scales. For example, the macroscopic properties of a material such as toughness, strength and ductility are strongly influenced by small scale defects like cracks and dislocations, which are initiated and evolve at the micro and nano scales. Hence, in the ambitious aim to derive the overall full-scale global mechanical response using a bottom-up approach, the sub-scale behavior has to be accurately computed. Although molecular dynamics (MD) simulations promise to reveal the fundamental mechanics of material failure by modeling the atom interactions, they are still prohibitively expensive to be employed in industrial applications [1,2]. Therefore, a plausible alternative to reduce the computational demand is to couple the continuum scale with the discrete scale using a multiscale approach. In this concern, the Quasi-Continuum Method (QCM) developed in [3] constitutes a new frontier for the formulation of novel multiscale methods coupling atomistic and continuum domains. In the QCM, the continuum degrees of freedom need to be located at the positions of the atoms at the interface, requiring a very fine grading of the continuum mesh around the defects. In two scale coupling, concurrent multiscale methods are mainly distinguished based on the coupling domain as the ‘Interface’ or the ‘Handshake’ coupling. The coupling is achieved across the boundary in the former case, whereas the regions are coupled over a finite overlapping domain in the latter. Classical examples of the ‘Interface’ and ‘Handshake’ couplings are the bridging scale method (BSM) and the bridging domain method (BDM), respectively. In particular, the BSM is based on the projection of the molecular dynamics (MD) solution onto the coarse scale shape functions to effectively address the spurious wave reflections in dynamic settings [4,5]. Through the use of the BSM and the Virtual Atom Cluster (VAC) model, bending of carbon nanotubes has been simulated in a two scale framework in [6]. A variation of the previous technique is the so-called adaptive multiscale method (AMM) for quasi-static crack growth which combines the VAC model [6], the phantom node method (PNM) [7–9], and the enhanced BSM [10,11]. The AMM approach considered that the coarse scale domain occupies the whole domain in the BSM, whereby the coupling between both scales is performed by enforcing displacement boundary conditions on the ghost atoms due to the fact that they follow the motion of the continuum. Therefore, using this approach, the coarse scale and fine scale problems can be solved independently in distinct computations. Regarding the BDM [12], this technique is based on a domain decomposition and a linear energy weighting procedure in the bridging domain. One advantage of the BDM over other methods relies on the fact that the nodes on the “continuum–atomistic” region do not need to be coincident with the atoms. Due to its versatility, the BDM has been also applied to dynamic problems in [13]. In this setting, Gracie et al., [14,15] have extended the bridging domain method (XBDM) to effectively account for dislocations and cracks. Further extensions of the XBDM to model cracks and dislocations in three dimensions can be found in [16,17]. A computational library of multiscale modeling of material failure has been proposed in [18]. Apart from the previous methodologies, an alternative multiscale strategy is proposed in [19], which is based on the principles of variational multiscale method (VMM) [20] through the exploitation of the concept of splitting the displacement field into large and small scale components. The small scale displacements were locally supported by assuming appropriate constraint conditions. An embedded statistical coupling method to couple MD atoms with finite element (FE) nodes with a statistical averaging of atomistic displacements in local atomic volumes associated with each FE node in an interface region has also been proposed in [21]. In this context, the finite element method (FEM) and MD computational systems are defined independent from each other and the interaction was included via an iterative update of the boundary conditions. A heterogeneous multiscale method by explicitly coupling the atomistic/continuum interface multiscale model to study the dynamics of brittle cracks in crystalline solids has been developed in [22]. Most of the previous multiscale methods have been also applied to modeling physical phenomena different from fracture in solids. This is the case of the technique developed by Molinari and coauthors within the framework of 340 P.R. Budarapu et al. / Comput. Methods Appl. Mech. Engrg. 319 (2017) 338–365 BDM [23], whereby a direct multiscale method coupling MD and FE simulations to investigate the contact area evolution of rough surfaces under normal loading has been proposed. Particularly, this approach has been considered to address the difficulties with regard to dealing with higher temperatures in the bridging domain. Many of the aforementioned multiscale methods do not adaptively adjust the fine scale domain as the defects propagate. Adaptive multiscale methods have been significantly improved in the last few years following the numerical procedures developed in [24–26,15], to quote a few of them. An adaptive multiscale technique to simulate crack propagation and crack coalescence, based on the extended finite element method (XFEM), has been proposed in [27]. A three-dimensional automatic adaptive mesh refinement for crack propagation based on the modified super convergent patch recovery technique by applying the asymptotic crack tip solution and using the collapsed quarter-point singular tetrahedral elements at the crack tip region has been developed in [28,29]. This latter exploits the use of aposteriori error estimator, and therefore the error of fracture parameters can be assessed and the crack path pattern can be accurately predicted. Recently, an efficient coarse graining (CG) technique to efficiently convert a given atomistic region to an equivalent coarse scale region has been developed in [30]. The vast majority of the previous methodologies are mostly suitable for 2D applications. However, over the years, modeling of three-dimensional complex fracture patterns which can undergo coalescence and branching in the continuum remains a challenging task as a consequence of the arising difficulties, especially from geometric signature. These problems are indeed relevant for complex technological applications such as Silicon photovoltaics, where crack branching and complex crack patterns due to impacts are observable via the electroluminescence technique in thin solar cells embedded into photovoltaic modules, see e.g. [31–34]. These applications usually regard the use of thin-walled structures (shells), which endow additional difficulties due to their slender and curved shape. Analyzing the different modeling options for fracture without enrichment, on the one hand, the popular interface cohesive fracture method (CFM), which is based on the incorporation of interface elements in the FE mesh [35–39], usually requires the a-priori knowledge of the crack path in its implicit version. In particular, as was addressed in [33], the CFM can be especially useful in order to assess the effect of crack opening on the electric response of solar cells. On the other hand, the strong discontinuity approach (SDA) [40] could be used for triggering fracture, but the main problem of crack nucleation would remain unsolved. An analogous drawback affects several formulations based on enriched FE techniques using the partition of unity (PU) concept such as XFEM [41,42] and the phantom node method (PNM), whose extension for the analysis of fracture in plates and shells can be found in [43–46]. In this context, the formulation of the PNM for thin shells has been proposed in [47,48], whilst a combination of XFEM and solid-like shell (solid shells) element topologies to study the delamination in composite structures has been developed in [49]. Complying with this technique, a discontinuous solid shell element for fracture in composites, accounting for the thickness stretching when fracture capabilities are embodied, has been proposed in [50]. In addition to the previous strategies, alternative techniques to model fracture in shells for various applications can be found in the related literature. In particular, a FEM-based computational method for the fracture of plates and shells on the basis of edge rotation and load control has been proposed in [51]. In this latter approach, the authors considered the crack front nodes as rotation axes, which results in each crack front edge in surface discretizations affecting the position of only one or two nodes. A meshfree method for thin shells with finite strains and arbitrary evolving cracks has been described in [52], where the authors eliminated the membrane locking by the use of a cubic or fourth-order polynomial basis. However, third-order completeness was necessary to remove membrane locking [53], which resulted in the use of very large domains of influence that made the method computationally expensive. A modified method using an extrinsic basis to increase the order of completeness of the approximation to reduce the computational cost has been developed in [53]. Tackling different applications, the fluid–structure interaction of fracturing structures under impulsive loads is described in [54], where a Kirchhoff–Love shell theory is adopted to model the structure and the cracks are treated by either discrete or continuous discontinuities. Recently, a novel phase-field model for fracture in Kirchhoff–Love thin shells using the local maximum-entropy (LME) meshfree method is described in [55] allowing complex fracture patterns via a smeared crack representation. An extended isogeometric element formulation (XIGA) for analysis of through-the-thickness cracks in thin shell structures based on Non-Uniform Rational B-Splines (NURBS) is proposed in [56], in which the singular field near the crack tip and the discontinuities across the crack are simulated based on the Kirchhoff–Love theory. According to the previous overview on the current modeling procedure in shells, in the current investigation a novel solid shell-based adaptive multiscale numerical method coupled with molecular statics to resolve the fine scale nonlinear phenomena at the crack tips is developed. Regions around crack tips are explicitly modeled in the atomistic P.R. Budarapu et al. / Comput. Methods Appl. Mech. Engrg. 319 (2017) 338–365 341 scale, whilst a self-consistent continuum model is employed elsewhere. Cracks in the continuum scale are simulated based on the PNM. In the subsequent developments, with the aim of clarifying the concepts used throughout the article, the continuum domain is denoted as the “coarse scale region”, whilst the atomistic sub-domain is denominated as the “fine scale region”. On the theoretical side, the current discontinuous solid shell model is formulated through the postulation of the mixed Hu–Washizu variational principle [57]. Trapezoidal and transverse shear locking are removed through the assumed natural strain (ANS) method [58]. The solid shell concept herein exploited allows the use of unmodified three-dimensional constitutive formulations within the computations and is locking-free, tackling the following pathologies: (i) membrane locking, (ii) Poisson thickness locking, (iii) volumetric locking, (iv) trapezoidal locking, and (v) transverse shear locking. Special attention is devoted to the numerical treatment along with the corresponding finite element formulation. Integrating these aspects, the proposed formulation, which represents a progress in the line of multi-scale fracture methods for thin structural elements, is applied to a series of test problems to show its capability to predict crack nucleation, growth, and complex fracture patterns. The manuscript is organized as follows. The solid shell-based three dimensional multiscale method is introduced in Section 2. An overview of the mathematical formulation of the discontinuous solid shell exploiting the PNM is provided in Section 3. Coupling conditions and the solution algorithm are discussed in Section 4. The present multiscale method is validated through 3D numerical examples in Section 5. The key contributions are summarized in Section 6, along with perspective applications to the field of photovoltaics. 2. Solid shell-based three dimensional multiscale method: coupling procedure In this section, the central aspects of the solid shell-based three dimensional multiscale method for the adaptive simulation of quasi-static crack growth are outlined. In the current modeling framework, whereas the coarse scale is considered through the adoption of a three dimensional solid shell model that allows the use of fully three dimensional constitutive formulations, the energy minimization in the fine scale domain is carried out using the open source Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) software [59]. The atom–atom interactions of Silicon in the fine scale are modeled using the Tersoff potential function [60] since Silicon photovoltaics is the target application. Tersoff potential has been successfully applied to predict mechanical properties of Graphene [61–64]. However, it is worth mentioning that other potentials can be used for other materials without any loss of generality. The initial size of the fine scale domain is chosen such that all the mechanical characteristics of the crack growth process around the crack tip are captured. Therefore, the initial fine scale domain size should include sufficient region ahead and behind the crack tip. Thus, a very large initial domain can lead to higher computational costs, whereas a very small initial domain can originate jump of the crack tip out of small fine scale domains. Note however that the selection of the size of the fine scale region are influenced by some model parameters: (1) type of problem (static/dynamic), (2) geometry and boundary conditions and (3) rate of loading and the type fracture (brittle/ductile). Hence, we followed the guidelines mentioned in [10] and tested several sizes before finally arriving at the domain sizes used in the present work. Specific details about the preliminary studies regarding the size of the fine scale are omitted for the sake of brevity, though a comprehensive parametric study was carried out till achieving numerical convergence in terms of the estimation of crack growth and computational efficiency for the applications herein presented. With reference to characteristics of the interaction between both scales, the coupling procedure between the coarse and the fine scales is realized by enforcing the displacement boundary conditions on the ghost atoms, in line with the bridging scale method (BSM). Ghost atom positions are interpolated based on the coarse scale solution. The BSM has been enhanced in the present study in order to account for the presence of cracks. To this aim, a user-defined MATLAB interface has been developed to activate LAMMPS in each load/time step. The LAMMPS input file is suitably updated with the positions of atoms determined from the latest deformed configuration. Correspondingly, the ‘crack tip’ in an atomistic domain is identified as the intersection of atoms on the crack surface on either side of the crack. Since the atoms on the crack surface possess the highest energy, they are identified based on an energy criterion. Adaptive refinement and coarse graining schemes are activated depending on the location of the crack tip. The above steps led to an adaptive continuum–atomistic multiscale method in the framework of enhanced BSM for crack growth in three dimensions, using the phantom node method to model crack propagation in the continuum, based on the solid shell. The novelties of the present method include: (i) a novel multiscale method coupling the 342 P.R. Budarapu et al. / Comput. Methods Appl. Mech. Engrg. 319 (2017) 338–365 Fig. 1. (a) Schematic of a three dimensional coupled continuum–atomistic model. (b) Mechanics of coarse scale domain modeled with solid shell element. (c) Fine scale region showing the arrangement of atoms in the diamond cubic lattice structure of Silicon. Crack in the coarse scale region is modeled using the Phantom node method and the fine scale model is embedded at the crack tip. continuum based hybrid solid shell with molecular dynamics in the framework of extended bridging scale method, adopting the Phantom node method to model crack in the continuum; (ii) a novel multiscale interface developed in MATLAB triggering LAMMPS through the system command; and (iii) Application of the present method to study complex fracture observed in monocrystalline Silicon photovoltaic solar cells. To the best of the authors’ knowledge, the proposed method is a key novel technique as compared to the several multiscale strategies reviewed in Section 1. In order to illustrate the proposed methodology, a simple representation is discussed in the sequel. Consider a three dimensional multiscale model with an initial crack as shown in Fig. 1(a). The shaded area in Fig. 1(a) corresponds to the solid shell-based coarse scale approximation. Squares denote the finite element nodes in the continuum discretization. The solid shell elements used in modeling the continuum domain along with the governing equations are shown in Fig. 1(b), whereas the particular arrangement of atoms in the diamond cubic lattice structure of Silicon used in the fine scale model (see Fig. 1(a)) is depicted in Fig. 1(c). The solid circles in Fig. 1(a) and (c) represent the atoms in the atomistic model, where the rose and pink atoms in Fig. 1(a) indicate the ghost atoms and the atoms on the crack surface, respectively. The material response at the crack tip is expected to be highly non-linear and/or non-homogeneous, and away from the tip it is expected to be homogeneous. The initial crack in the fine scale region is created by deleting the bonds between the atoms on the crack surface and updating the neighbor list accordingly. The neighbor list is generated based on a radius of influence which can be provided as an input data. The bond potential corresponding to the atomistic model is calculated using the potential function, which depends on the distance between the atoms. The potential energy of an atom of interest is estimated from the interaction potential function, which is based on the distance to its neighbors within the domain of influence. Thus, in order to create a crack, atoms on one side of the crack surface are restricted to interact with the atoms laying on the other side of the crack surface. In the present methodology, this is achieved by creating regions on either side of the crack surface and restricting the interactions between them. Additionally, ghost atoms are located in the coarse region but within the cutoff radius of the atoms in the fine region. Their positions are interpolated from the coarse scale solution and enforced as the boundary conditions for the fine scale problem. The crack originates from the coarse scale region with the crack tip captured in the fine region. The fine scale region is adaptively adjusted as the crack propagates, following the method proposed in [10]. Finally, the centro symmetry parameter (CSP) is used in order to identify the atoms on the crack surface and hence to detect the crack tip location [30]. Specific details with regard to the estimation of the total displacement field of the coupled model and the internal forces in the fine scale region are addressed in Appendix A and Appendix B, respectively. P.R. Budarapu et al. / Comput. Methods Appl. Mech. Engrg. 319 (2017) 338–365 343 3. Coarse scale continuum formulation In this section, a brief introduction to the fundamentals of the shell formulation using a continuous and discontinuous kinematic description is given. The interest for the development of the proposed multiscale methodology incorporating the so-called solid shell concept relies on the fact that this parametrization is one of the most popular strategies for shells, see [57,65–68] for more details. In particular, this approach has several attractive aspects from the operative standpoint: (i) the complete kinematic compatibility with respect to standard continuum elements and contact/interface element formulations [69]; (ii) the avoidance of rotational degrees of freedom to update the shell director vector along the deformation process; and (iii) the use of three-dimensional constitutive laws without modifications. In the present investigation, fracture events in the continuum are modeled by means of a locking-free discontinuous solid shell formulation, whose kinematics relies on the PNM. The current element formulation includes the use of both the ANS and the EAS methods to remove locking deficiencies according to the formulation developed in [57]. The seminal idea for the development of EAS-based shell formulations regards tackling the so-called Poisson thickness locking in thin structures using three dimensional constitutive laws. Note also that this deficiency can be also circumvented by means of embodying a quadratic displacement interpolation over the shell thickness [70]. Nevertheless, due to the potential and adaptability of the EAS method, further common locking pathologies can be remedied using this numerical strategy. This is the case of formulations proposed by different authors [57,68,71], which alleviate the following pathologies: (i) transverse thickness locking, (ii) volumetric locking, and (iii) membrane locking. The previous techniques are subsequently combined with the PNM to account for arbitrary cracks within the shell body. In this context, Dolbow and Devan [72] integrated the EAS method with discontinuous enrichment to alleviate volumetric locking in finite strain applications using an enhanced deformation gradient in the form of that proposed in [73]: F:= Fu  enrich. with discontinuity +˜ F  incompatible ,(1) where Fuidentifies the displacement-derived deformation gradient including the discontinuous kinematics and ˜ F stands for the incompatible deformation gradient. Alternatively, the enhancement scheme relying on the additive decomposition of the Green–Lagrange strain tensor [57,71] is pursued, as: E:= Eu  enrich. with discontinuity +˜ E  incompatible and enriched ,(2) where, similarly, Euand ˜ Eaccount for the compatible with discontinuous displacements and the incompatible counterparts of the Green–Lagrange strain tensor, respectively. 3.1. Kinematic formulation Within the finite deformation setting, let X(ξ1, ξ2, ξ3)∈ΩC 0denote the position vector of an arbitrary material point in the reference configuration Ω0at time t0, and x(ξ1, ξ2, ξ3)∈ΩCindicating the corresponding position vector in the current configuration ΩCat time t∈R+. In the sequel, the superscript Creferred to the coarse scale is removed to alleviate the notation. The parametric curvilinear coordinates are denoted by ξ={ξ1, ξ2, ξ3}, with ξi∈[−1,1] (i=1,2,3), see Fig. 2. The co-variant basis in the reference (Gi) and current (gi) configurations is defined as Gi(ξ):= ∂X(ξ) ∂ξi;gi(ξ):= ∂x(ξ) ∂ξi,(3) where Gi·Gj=δj iand gi·gj=δj i. The metric tensors in the reference and current configuration are defined as: G=Gi j Gi⊗Gj=Gi j Gi⊗Gj;g=gi j gi⊗gj=gi j gi⊗gj.(4) 344 P.R. Budarapu et al. / Comput. Methods Appl. Mech. Engrg. 319 (2017) 338–365 Fig. 2. Definition of shell body, where {Ei}i=1,3,{ei}i=1,3denote the standard Cartesian settings in the reference Ω0and current Ωconfigurations, respectively. The coordinates ξ= {ξ1, ξ2, ξ3}define the parametric space. The parametrization of the shell body through the solid shell representation can be expressed in terms of the materials points on the top Xt(ξ1, ξ2) and the bottom surfaces Xb(ξ1, ξ2): X(ξ1, ξ2, ξ3)=1 2(1+ξ3)Xt(ξ1, ξ2)+1 2(1−ξ3)Xb(ξ1, ξ2).(5) Analogously, in the current configuration, this representation yields x(ξ1, ξ2, ξ3)=1 2(1+ξ3)xt(ξ1, ξ2)+1 2(1−ξ3)xb(ξ1, ξ2).(6) Therefore, the kinematic field, u, adopts the following form u(ξ1, ξ2, ξ3)=x(ξ1, ξ2, ξ3)−X(ξ1, ξ2, ξ3)=v(ξ1, ξ2)+ξ3w(ξ1, ξ2),(7) where v(ξ1, ξ2) denotes the displacement of the shell mid-surface and w(ξ1, ξ2) is the difference vector mapping the shell director vector between the reference and the current configurations as given below: v(ξ1, ξ2)=1 2[ut(ξ1, ξ2)+ub(ξ1, ξ2)],w(ξ1, ξ2)=1 2[ut(ξ1, ξ2)−ub(ξ1, ξ2)].(8) The displacement-derived deformation gradient (Fu) is defined as: Fu:= gi⊗Gi,(9) and the corresponding displacement-derived Green–Lagrange deformation tensor is expressed as: Eu:= 1 2[(Fu)TFu−G]=1 2[gi j −Gi j ]Gi⊗Gj. (10) The components of stress (S) and strain (E) tensors in the Voigt notation are arranged as follows: S=[S11,S12,S13,S22,S23,S33]T;E=[E11,2E12,2E13,E22,2E23,E33]T.(11) 3.2. The enhanced assumed strain (EAS) method The variational basis of the present strategy relies on the multi-field Hu–Washizu variational principle for the fields {u,˜ E,S}, where Sdenotes the second Piola–Kirchhoff stress tensor, which is energetically conjugated to the Green–Lagrange strain tensor. Hence, the corresponding boundary-value problem can be defined as: find {u,˜ E,S}, with Vu={δu∈[H1(Ω0)] :δu=0on ∂Ω0,u}identifying the space of admissible displacement variations, and P.R. Budarapu et al. / Comput. Methods Appl. Mech. Engrg. 319 (2017) 338–365 345 Fig. 3. Phantom node method (PNM): schematic representation of a cracked element containing (d) straight and (e) angled cracks. V˜ E,VS=[L2(Ω0)], respectively standing for the admissible space for the enhanced strain and stress fields, such that: R(u, δu,˜ E, δ ˜ E,S, δS)=∫Ω0[∂Ψ(E) ∂E:∂Eu ∂uδu]dΩ+∫Ω0[∂Ψ(E) ∂E:δ˜ E]dΩ+∫Ω0 S:δ˜ EdΩ +∫Ω0 δS:˜ EdΩ+δΠext(u)=0,∀(δu, δ ˜ E, δS)∈Vu×V˜ E×VS.(12) In Eq. (12),δΠext denotes the virtual contribution stemming from the prescribed external actions; {δu, δ ˜ E, δS} denote arbitrary variations of the corresponding fields u,˜ Eand S, whilst Ψ(E) identifies the Helmholtz free-energy function, which is a function of the Green–Lagrange strain tensor and can accommodate any material model. This formulation can be further simplified by invoking the orthogonality condition between the enhanced strain and the stress field, as explained in [57]. Through the use of standard arguments, the second Piola–Kirchhoff stress tensor can be defined as: S(E):= ∂Ψ/∂E=∂EΨ. Specific details concerning the FE discretization of the current enhanced-based solid shell element are given in Appendix C and Appendix D. 3.3. Phantom node method (PNM) for solid shells This section briefly outlines the main aspects with regard to the proposed PNM formulation for solid shells. In particular, the discontinuous formulation herein proposed relies on an eight node solid shell element. Note that differing from [50] the current model does not incorporate any internal geometrical node, since the different locking pathologies are alleviated through the EAS method. Let us consider an arbitrary shell body with a surface of discontinuity Γc, as shown in Fig. 3(a)–(b). According to PNM [7], the kinematics of a cracked element can be described by superimposing two separate displacement fields, which are active only in a determined region of the domain. Consequently, a completely cut element can be represented as a union Ωph. elem 0=Ωelem1 0∪Ωelem2 0, of two elements separated along the crack surface, see Fig. 3(d)–(e). The superscript ‘ph. elem’ refers to the considered phantom element and ‘elem1’ and ‘elem2’ denote the sub elements after splitting, see Fig. 3(d)–(e). This formalism is expressed by setting that the crack surface divides the shell domain into two sub-domains, viz. Ω0=Ω0(+)⋃Ω0(−). Correspondingly, two phantom domains are defined: Ωp 0=Ωp 0(+)⋃Ωp 0(−). Since the elements in the two sub-domains do not share any nodes in common, their displacements are independent, resulting in the expected discontinuity across the surface. In the corresponding FE discretization it is usual to use the terminology of real and phantom nodes in each of the subdomains to make reference to either active or inactive degrees of freedom. Through the definition of fas the signed distance measured from the crack surface, W+ 0,W− p,W− 0and W+ pas the nodes belonging to Ω0(+),Ωp 0(−),Ω0(−)and Ωp 0(+), respectively, the discontinuous interpolation of the displacement field is given by: u(X,t)=∑ I∈{W+ 0,W− p} uI(t)NI(X)H(f(X)) +∑ J∈{W− 0,W+ p} uJ(t)NJ(X)H(−f(X)) (13) 346 P.R. Budarapu et al. / Comput. Methods Appl. Mech. Engrg. 319 (2017) 338–365 where His the Heaviside function. In line with [46,47], the standard approximation of the displacements on each part of the cracked element Ω0(+)and Ω0(−), which are extended to their corresponding phantom domains Ωp 0(−)and Ωp 0(+) introduces the continuous displacement field. Derived from the previous results, the enriched kinematic field given in Eq. (13) leads to the consideration of the following discontinuous operators [50]: (i) the deformation mapping, (ii) the deformation gradient, (iii) the compatible Green–Lagrange strain tensor, and (iv) the incompatible strain tensor: ϕ={ϕ1∀X∈Ω0(+) ϕ2∀X∈Ω0(−) ;F={Fu 1∀X∈Ω0(+) Fu 2∀X∈Ω0(−) (14) Eu={Eu 1∀X∈Ω0(+) Eu 2∀X∈Ω0(−) ;˜ E={˜ E1∀X∈Ω0(+) ˜ E2∀X∈Ω0(−).(15) According to the previous definitions, it can be seen that the displacement jump between the two flanks of the crack can be computed by taking the difference of the displacement fields of the two domains of the cracked element. Thus, from Fig. 3(d)–(e), let us define the displacements of element 1 and element 2 as: u1(X,t)=∑ I∈{W+ 0,W− p} uI(t)NI(X)H(f(X)) (16a) u2(X,t)=∑ J∈{W− 0,W+ p} uJ(t)NJ(X)H(−f(X)).(16b) Therefore, based on Eq. (13), the total displacement field can be expressed as the summation of displacements of element 1 and element 2. However, from Eqs.(7) and (8) one obtains u1(X,t)=v1+ξ3w1=1 2[(ut1+ub1)+ξ3(ut1−ub1)],(17a) u2(X,t)=v2+ξ3w2=1 2[(ut2+ub2)+ξ3(ut2−ub2)].(17b) The substitution of Eqs. (17a) and (17b) into Eq. (13), yields: u(X,t)=(1 +ξ3)(ut1+ut2 2)+(1 −ξ3)(ub1+ub2 2).(18) To express the weak form of Eq. (12) by considering kinematic discontinuity due to the presence of cracks, we recall the additive property of integrals [50] and we exploit the minimization of this functional with respect to the independent fields, i.e. u1,u2,˜ E1and ˜ E2. Using the discontinuous kinematic definition introduced in Eq. (13), the variational form given in Eq. (12) can be expressed as: R(+)(u1, δu1,˜ E1, δ ˜ E1)=∫Ω0(+) S:δEu 1dΩ+∫Ω0(+) S:δ˜ E1dΩ+δΠext(+)(u1)=0,(19) R(−)(u2, δu2,˜ E2, δ ˜ E2)=∫Ω0(−) S:δEu 2dΩ+∫Ω0(−) S:δ˜ E2dΩ+δΠext(−)(u2)=0.(20) In Eqs.(19) and (20), the stress field has been already removed from the formulation through the aforementioned orthogonality condition between interpolation spaces associated with the stress and the incompatible strain fields. Furthermore, it is worth mentioning that following [57], the above approach requires the definition of the incompatible strains in both of the two-subdomains leading to a duplication of the incompatible strains. Moreover, the cracked elements are integrated over their corresponding active domains. After the insertion of the discretization schemes corresponding to the displacements and incompatible strains outlined above, and performing the consistent linearization of the residual equations at each of the domains Ω0(+) and Ω0(+), the final system of equations at the element level reads: [kdd kdς kςdkςς ][∆d ∆ς]=[fext 0]−[fint fEAS].(21) P.R. Budarapu et al. / Comput. Methods Appl. Mech. Engrg. 319 (2017) 338–365 353 Fig. 8. Adaptive refinement and coarse graining of the fine scale region as the crack grows. (a) Deformed configuration of the multiscale model after 38 load steps. A close up of the regions around the (b) left crack tip and the (c) right crack tip of the deformed configuration in (a), showing the direction the crack growth. (d) A three dimensional picture of (a), showing the atoms in the thickness direction. (e) Area around the vertical edges of the fine scale regions not containing the crack. High energy atoms (lying on the crack surface) are observed to be touching the boundaries of the fine scale region, which is the right time for adaptive refinement. (f) Multiscale model after an adaptive refinement after 39 load steps. Adaptive refinement and coarse graining algorithms (see [10]) are activated after 39 load steps as the cracks grow. As a result, the two fine scale regions are merged after 39 load steps and the combined fine scale region is adaptively adjusted, as plotted in (g)–(i). multiscale model at the end of the simulation is shown in Fig. 8(i). In this plot, almost a complete merging of cracks and hence the separation of the fine scale region into two parts can be noticed. 354 P.R. Budarapu et al. / Comput. Methods Appl. Mech. Engrg. 319 (2017) 338–365 Fig. 9. Schematic of a three dimensional coupled continuum–atomistic model to simulate the mode III crack propagation. An initial edge crack oriented along the x-direction, located in the middle of the domain length along the ydirection is created to study the growth. An atomistic region (ΩA) is created around the crack tip to capture the growth mechanics. Displacement boundary conditions are specified on the left and right edge nodes in the continuum region (ΩC) whereas, the nodes on the top and bottom edges are arrested along the y-direction. 5.2. Example 2: Out-of-plane crack growth The second application under consideration regards an out-of-plane (Mode III) crack propagation of a slab with an edge crack along the x-direction. The initial crack is located in the middle of the domain along the y-direction, see Fig. 9. An atomistic region ΩAis considered around the crack tip. Positive and negative displacements along the z-direction are specified to the left edge nodes on either side of the crack surface, in the continuum region (ΩC). Nodes with specified displacements on the left edge are allowed moving along the x-direction only. All the degrees of freedom of the right edge nodes are restrained, whereas the top and bottom edge nodes are allowed to move along the xand zdirections only, see Fig. 9. Consider a three dimensional coarse scale model with dimensions 120.0˚ A×80.0˚ A×14 ˚ A. An initial edge crack of length 40.0 ˚ A along the x-direction, located at 40.0 ˚ A in the middle of the domain along the y-direction, is created in the coarse region, see Fig. 9. Therefore, the initial crack tip is located at (40 ˚ A, 40 ˚ A). The model is discretized with 46 ×46 nodes along the xand ydirections, respectively. The PNM complying with the procedure outlined in Section 3.3 is employed to model the crack in the coarse region. In the current initial model there are 14 completely cracked elements and one tip element with 60 phantom nodes in total. Mode III crack propagation involves large displacements applied in the out-of-plane direction. Consequently, a single large atomistic region is considered in the initial multiscale model. An initial fine scale region measuring 84.0˚ A×33.94 ˚ A×14.0˚ A with 1787 active atoms and 468 ghost atoms is created including the crack tip position as shown in Fig. 10(a). This plot also shows the geometry and the location of the initial edge crack, apart from the highlighted ghost atoms. A uniform displacement load of 120 ˚ A along the z-direction is specified on the left edge nodes, in 210 equal pseudo-time steps. Displacements on either side of the crack surface on the left edge are specified in order to trigger Mode III crack growth (see Fig. 9), so that the crack opens and propagates. Fig. 10(b) shows the deformed configuration of the multiscale model after 96 load steps, whereas a three dimensional view showing the deformation of the coarse and fine scales along the thickness direction is provided in Fig. 10(c). Another isometric view of Fig. 10(b) is shown in Fig. 10(d). The distribution of the potential energy in an isometric view of the isolated atomistic region after 96 load steps is shown in Fig. 10(e), whilst Fig. 10(f) shows a closeup of the atoms around the crack tip, where the range of potential energy is shown in the color bar in Fig. 10(g). A close observation of the fine scale region in Fig. 10(e) and (f) indicate that the bonds of the atoms around the crack tip are not broken and hence the crack is not yet about to propagate. A close up view of the region in the multiscale model around the crack tip after 124 load steps is shown in Fig. 10(h), whereas Fig. 10(i) depicts an isolated picture of the deformed atomistic region around the crack tip after P.R. Budarapu et al. / Comput. Methods Appl. Mech. Engrg. 319 (2017) 338–365 355 Fig. 10. Mode III crack propagation in the multiscale frame work. (a) Initial configuration of the multiscale model showing the discretized coarse region, geometry of the crack and the atoms in the fine scale region along with the ghost atoms along the boundaries of the fine scale region. (b) Deformed configuration after 96 load steps. Two different views of the deformed configuration in (b) are shown in (c) and (d). Distribution of the potential energy in the (e) fine scale region and (f) a closeup of the atoms around the crack tip, extracted from the deformed configuration of the multiscale model in (b). (g) A color bar showing the range of potential energy plotted in (e) and (f). (h) A close up of the region around the crack tip after 124 load steps, where the isolated fine scale region of (h) is plotted in (i). A close up of the fine scale region around the crack tip is shown in (j). Two different views of the deformed configuration of the coarse region at the end of the simulation are plotted in (k) and (l), where (m) shows the isolated fine scale region in (k). 356 P.R. Budarapu et al. / Comput. Methods Appl. Mech. Engrg. 319 (2017) 338–365 Fig. 11. (a) Distribution of the potential energy with strain and (b) load–displacement diagram observed during the simulation of mode III crack growth in the second numerical example. 124 load steps. Analyzing in detail these results, it can be seen that Fig. 10(i) highlights the significant elongation of the bonds of the Silicon atoms around the crack tip that are about to break. The deformed atomistic region around the crack tip after 132 load steps is shown in Fig. 10(j). Performing a careful qualitative comparison of Fig. 10(i) and (j) it can be seen that after 132 load steps some bonds are observed to be broken and hence this marks the onset of crack propagation. Further application of the load leads to breakage of more bonds and hence crack growth. In comparison with the results regarding in-plane fracture, larger displacements are required to propagate the crack in Mode III. This is also confirmed from the potential energy distribution plotted in Fig. 11(a). According to this evolution it is observed that cumulative strain of the order of 1.48 is required to break the first bonds. However, due to the brittle character of Silicon, further propagation of the crack is very fast and unstable. Variation of the force with the total displacement, along the z-direction is plotted in Fig. 11(b). Note that, because of the selected boundary conditions in this numerical example (see Fig. 9), a complete material separation is not possible. In other words, the simulation was terminated when the crack tip reaches close to the right edge. In addition to the previous considerations, it is worth mentioning that based on the location of the crack tip in the fine scale model, the crack in the continuum is consistently grown by adding the phantom nodes to the newly identified split elements. The process is repeated until the end of the simulation. The deformed configuration of the continuum at the end of the simulation is shown in Fig. 10(k) and (l), respectively, whereas, a zoom around the crack tip of the deformed fine scale region is shown in Fig. 10(m). At the end of the simulation we noticed that the crack propagated horizontally, and the tip is located at ≈75% of the domain length along the x-direction. 5.3. Example 3: Crack growth due to indentation The last example regards a highly technological application, which is the case of Silicon solar cells embedded into photovoltaic laminates. These systems are prone to cracking at the Silicon (thin) layer, whilst the most preferential external actions lead to bending-dominated deformation patterns. This fact motivates the use of three-dimensional description of the crack pattern for thin structures. For instance, bending can be induced by a uniform snow pressure acting on the laminate simply supported along the edges, leading to a complex crack pattern as shown in Fig. 12. Alternatively, hail impacts can induce localized cracked areas under the punched region, see Fig. 13. In both cases, the crack pattern identified using the electroluminescence technique [32–34] influences the electric power output, since the black areas are electrically insulated areas of the solar cell not contributing to solar energy conversion. As experimentally shown in [32], cracks are not all fully insulated, and the degree of insulation is dependent on the crack opening, which can be suitably passed as input to an electric model to predict the current–voltage response of the photovoltaic laminate, see [33]. Therefore, an accurate assessment of crack opening is essential and a global–local finite element formulation has been proposed in [33] in the case of bending. Clearly, in the case of hail impacts, damage and crack growth is induced by relative out-of-plane (Mode III) displacements between the loaded punch P.R. Budarapu et al. / Comput. Methods Appl. Mech. Engrg. 319 (2017) 338–365 357 Fig. 12. Example of crack pattern in Silicon solar cells due to bending. Fig. 13. Example of crack pattern in Silicon solar cells due to hail impact. and the rest of the solar cell [34]. For such a case, the present framework offers an excellent possibility to accurately predict the evolution of material separation during impact by using the atomistic model. To show an example, a uniform out-of-plane displacement along the z-direction is applied to the upper side of the shell in a specific area. Relying on symmetry arguments, only a quarter of the shell is modeled. Symmetry boundary conditions are applied on the right and bottom edges, whereas all the degrees of freedom on the left and top edges are restrained. Consider a three dimensional atomistic model with dimensions 720.0˚ A×720.0˚ A×44 ˚ A to embody a circular indentation area of diameter 94.0 ˚ A. A quarter model of the shell has in plane dimensions 180.0˚ A×180.0˚ A×44 ˚ A. The quarter model consists of 76 296 atoms in total. A quarter circle with a radius of 47.0 ˚ A has been created with the lower right corner as the center. The quarter circle is further extruded to a quarter cylinder along the thickness direction. Atoms on the upper side of the quarter cylinder portion are subjected to uniform imposed displacements along the z-direction. The simulation is carried out in 1530 load steps in total, by specifying a displacement of 0.05 ˚ A in each step. Figs. 14(a) and (b) shows the deformed configuration after 100 load steps. From these graphs it can be observed that atoms in the first few layers are separated from the atoms on the cylinder, indicating the initiation and the growth of the separation between the material inside the punched area from the remainder of the shell. An isometric view 358 P.R. Budarapu et al. / Comput. Methods Appl. Mech. Engrg. 319 (2017) 338–365 Fig. 14. Simulation of indentation in a shell based on the atomistic model. Distribution of the potential energy at various instances of the punching process. (a) An isometric view and (b) a side view of the deformed configuration after 100 load steps. (c) and (e) An isometric view of the deformed configuration after 500 and 1530 steps. (d) and (f) A zoom of the region in (c) and (e) respectively, around the punching area. after 500 steps is shown in Fig. 14(c), where almost 50% of the material in the punched region is separated from the shell. A zoom of the region around the quarter cylinder in the top view is shown in Fig. 14(d), whereas a clear step separating the cylinder and remainder of the shell can be seen in Fig. 14(d). Progressing on the loading application, Fig. 14(e) corresponds to the deformed configuration after 1530 steps. Analyzing in detail these results, it can be seen that almost the entire material under the punch is separated from the plate material, as shown in the zoomed picture in Fig. 14(f). Finally, the evolution of the potential energy with strain is plotted in Fig. 15(a). In this graph it is observed that the potential energy seems to fluctuate after the strain is reaching 0.48. These fluctuations correspond the breaking of bonds on each layer of atoms in the plate with the corresponding atoms in the cylinder. Correspondingly, a drop in the force can be observed as plotted in Fig. 15(b). After breaking few initial bonds, lower forces are required for P.R. Budarapu et al. / Comput. Methods Appl. Mech. Engrg. 319 (2017) 338–365 359 Fig. 15. (a) Distribution of the potential energy with strain and (b) load–displacement diagram generated based on the MD simulations of hole punching in a rectangular panel. further material separation, see Fig. 15(b). A complete separation and further movement of the punched cylinder can be observed in Fig. 14(e) and (h) after 1530 steps. Thus, a continuous drop of the potential energy after the strain value of about 1.0 indicates the complete material separation and hence a rigid body motion occurs (though not a sudden drop is obtained due to the fact that there are still some atoms attached to). Similar behavior can be observed in the load–displacement diagram as well, see Fig. 15(b). The outcome of atomistic simulations can be profitably used as input of the coarse scale shell model in view of global/local simulations. The coarse scale model based on the solid shell can be used to estimate the electric response of the solar cell by using a generalization of the electric model in [33] accounting for a localized electric resistance dependent not only on the in-plane crack opening displacement, but also on the out-of-plane relative displacement along the crack front. 6. Conclusions In this work, a new continuum-based PNM combined with an enhanced strain-based solid shell element has been coupled with a molecular statics (MS) model to generate a multiscale framework for the simulation of cracks in thin structures. The use of the PNM over the standard XFEM formulation has been motivated due to its simpler numerical implementation. Atomistic simulations in the fine scale region are carried out by triggering the LAMMPS software through the system command, which is integrated into an in-house MATLAB code. Cracks have been incorporated into the fine scale model by deleting atomic bonds, whereas the PNM has been used to simulate crack propagation in the coarse scale. The coupling of the continuum and atomistic models has been realized through the use of a bridging scale method (BSM). The original BSM has been additionally enhanced with respect to its original formulation, so that arbitrary cracks are admissible at the coarse scale using the PNM. Therefore, the proposed numerical strategy can be categorized as a novel three-dimensional adaptive multiscale method (3DAMM) for thin structures. To demonstrate the modeling capacity of the current 3DAMM, it has been used to simulate crack propagation in three scenarios with regard to Silicon specimens with a cubic diamond lattice structure containing propagating cracks. In the first example, in plane crack propagation of two edge cracks has been modeled based on the proposed method in order to show an application featuring multiple cracks. Adaptive refinement and coarsening schemes have been implemented during crack propagation to improve the performance by reducing the related computation effort. In the second example, a problem displaying out-of-plane crack growth has been simulated. Geometrical nonlinear effects have been considered in the computations, so that tearing mode has been achieved. In the third example, a mechanical indentation problem of Silicon solar cell has been simulated based on the atomistic model and the obtained results open new perspectives for the integration with the coarse scale solid shell model for the computation of the electric 360 P.R. Budarapu et al. / Comput. Methods Appl. Mech. Engrg. 319 (2017) 338–365 power degradation of the solar cell. In particular, the computed out-of-plane relative displacements along the crack front are envisaged to be the essential quantity predicted by the atomistic model, which can be employed as input to a generalized electric model of the solar cell in the presence of cracks subject to bending-dominated deformation patterns. Savings in computational times are observed within 50%–400% for the numerical examples considered in this article in comparison with the numerical strategy in [33]. Acknowledgments The authors acknowledge funding from the European Research Council (ERC), Grant No. 306622 through the ERC Starting Grant “Multi-field and multi-scale Computational Approach to Design and Durability of PhotoVoltaic Modules”—CA2PVM. JR is also grateful to the support of the Spanish Ministry of Economy through the projects (DPI2012-37187, MAT2015-71036-P, and MAT2015-71309-P) and Andalusian Government (Projects of Excellence P11-TEP-7093 and P12-TEP-1050). Appendix A. Estimation of the total displacement field of the coupled model Complying with the current two scale approach, the total displacement field uαof an atom αis decomposed into coarse and fine scale components [10]: uα=uC α+uA α,(A.1) where uC αis the coarse scale component within the domain ΩC, whereas uA αis the fine scale component within the domain ΩAwhose projection onto the coarse scale is zero. The fine scale component uA αis the difference between the actual position of an atom αand the interpolated position on the coarse scale. In other words, uA αis vanishing in the regions away from the crack tip, and hence uC αis sufficient to model the deformation in the coarse scale region. On the other hand, in the fine scale region, both coarse and fine scales components have to be considered. Thus, let the coarse scale displacement uC αof an atom αbe represented by a set of FEM basis functions defined over a set of nC nodal points: uC α= nC ∑ I=1 NI(Xα)uC I,(A.2) where NI(Xα) is the shape function defined at node I, estimated at the αth atom with the material coordinate Xα,uC I is the continuum displacement vector at node I, and nCdenotes the number of coarse scale nodes within ΩC. Note that the coupling procedure herein described is performed at each incremental time step due to the fact that the position of the coarse scale (shell structure) is updated along the solution process. The reference configurations at the coarse and fine scales are denoted by ΩC 0and ΩA 0, respectively. Additionally, it is worth mentioning that the material points at the reference coarse scale are denoted by X∈ΩC 0, whose positions are transformed into the current positions x∈ΩCvia the nonlinear mapping operator: x=ϕ(X,t). Appendix B. Estimation of the internal forces of the fine scale model In molecular statics (MS), the objective is to determine the positions of the atoms for the given boundary conditions, by minimizing the potential energy of the system, which is given by Π=Wint −Wext (B.1) where Wint represents the internal energy of the system and Wext is the external contribution. Consider the simplest atom–atom interactions in which the potential energy is only a function of the distance between two atoms, the total internal energy of the system is given by summing the energies of all the atomic bonds over all the atoms, as given below: Wint =1 2 nA ∑ α=1 nA ∑ β=α V(rαβ ) (B.2) where V(rαβ ) is the bond potential between the atoms αand β, separated by distance rαβ . The system potential energy is minimized, when the first derivative of the potential function with respect to the positions of the atoms is equal to P.R. Budarapu et al. / Comput. Methods Appl. Mech. Engrg. 319 (2017) 338–365 361 zero. Therefore, for any given atom λ, the first derivative of the system potential energy with respect to the position vector rλis ∂(1 2∑nA α=1∑nA β=αV(rαβ )) ∂rλ −∂Wext ∂rλ =0 (B.3) where the internal forces acting on atom λare given by Fint λ=1 2 nA ∑ α=1 nA ∑ β=λ −∂V(rαβ ) ∂rαβ ∂rαβ ∂rλ (B.4) and the external forces acting on atom λare Fext λ= −∂Wext ∂rλ .(B.5) The residual forces on each atom are R=Fint −Fext.(B.6) The distance rαβ in (B.4) is defined as rαβ = |rα−rβ| =    √ 3 ∑ j=1 (rαj−rβj)2(B.7) where jis the free index. Substituting Eq. (B.7) into Eq. (B.4) yields Fint α= − nA ∑ β=α ∂V(rαβ ) ∂rαβ (rα−rβ rαβ ).(B.8) Details of the derivation of Eq. (B.8) are explained in the Appendix of [10]. Appendix C. Finite element formulation of the coarse scale model This appendix section outlines the finite element formulation of the coarse scale model relying on an enhancedbased solid shell model. Regarding the kinematic field, based on the isoparametric parametrization, the discretization of the reference and current position vectors of the present 8-node solid shell element can be expressed as: X≈ 4 ∑ A=1 1 2(1 −ξ3)NA(ξ1, ξ2)Xb,A+1 2(1 +ξ3)NA(ξ1, ξ2)Xt,A=NXe,(C.1a) x≈ 4 ∑ A=1 1 2(1 −ξ3)NA(ξ1, ξ2)xb,A+1 2(1 +ξ3)NA(ξ1, ξ2)xt,A=Nxe,(C.1b) where NA(ξ1, ξ2) are the in-plane shape functions in the natural space: NA=1 4(1+ξ1 Aξ1)(1+ξ2 Aξ2),with ξ1 A, ξ2 A= ±1.(C.2) In Eqs.(C.1a) and (C.1b), the discrete top and bottom position vectors of the node Ain the reference and current configuration are expressed by the pairs (Xb,A,Xt,A) and (xb,A,xt,A), respectively. These nodal component vectors are gathered at the element level by the vectors Xeand xe, respectively. After multiplying the in-plane shape functions by those associated with the linear interpolation over the thickness [67], we obtain: NA=1 8(1+ξ1 Aξ1)(1+ξ2 Aξ2)(1+ξ3 Aξ3),with ξ1 A, ξ2 A, ξ3 A= ±1.(C.3) Therefore, the interpolation of the displacement field, its variation and increment can be expressed as: u≈Nd;δu≈Nδd;∆u≈N∆d,(C.4) 362 P.R. Budarapu et al. / Comput. Methods Appl. Mech. Engrg. 319 (2017) 338–365 Fig. D.16. Collocation points in the element parametric space: use of the ANS method to tackle transverse shear and trapezoidal locking. where Ncorresponds to the matrix operator that collects the shape functions at the element level, and ddenotes the nodal displacements. Similarly, the interpolation of the variation and the increment of the displacement-derived strain field is estimated through the compatibility operator B(d) as: δEu≈B(d)δd,∆Eu≈B(d)∆d.(C.5) The enhanced strain field based on the EAS method is employed here to remedy membrane, volumetric and Poisson thickness locking deficiencies. Therefore, the interpolation of the incompatible strains, the corresponding variation and increment can be expressed as: ˜ E≈M(ξ)ς, δ ˜ E≈M(ξ)δς, ∆˜ E≈M(ξ)∆ς, (C.6) where ςstands for the incompatible strain vector. The design of the interpolation matrix Mis chosen based on the combination of schemes proposed in [68,71] leading to 7 EAS parameters: ˜ M= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ξ10 0 0 0 0 0 0ξ1ξ20 0 0 0 00000 0 0 000ξ20 0 0 00000 0 0 0000ξ3ξ1ξ3ξ2ξ3 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ .(C.7) Finally, it is noted that the relationship between the matrices ˜ Mand Mis established by a transformation mapping from the natural space to the global Cartesian setting: ˜ E=[det J0 det J]T0˜ Mς=Mς, (C.8) where J,J0identify the Jacobian and its evaluation at the element center, respectively, and T0is a transformation matrix [68]. Appendix D. The assumed natural strain (ANS) method The ANS method in shells is a collocation procedure which has been extensively employed to tackle the transverse shear [74] and trapezoidal locking [75]. According to [74], the alleviation of transverse shear locking can be performed by defining four sampling points A, B, C and D located at the four mid-points of the element edges on the shell mid-surface ξ3=0. The coordinates of the previous collocation points in the parametric element space are given by: ξA=(0,−1,0),ξB=(1,0,0), ξC=(0,1,0)and ξD=(−1,0,0), refer Fig. D.16. Hence, the modified interpolation of the transverse shear strain components, E13 and E23 obeys, {2EAN S 13 2EAN S 13 }={(1 −ξ2)2E13(ξA)+(1 +ξ2)2E13(ξC) (1 +ξ1)2E23(ξB)+(1 −ξ1)2E23(ξD).}(D.1)