Peridynamic Shell Model Based on Micro-Beam Bond

Guojun Zheng,Zhaomin Yan,Yang Xia,Ping Hu and Guozhe Shen,★

1School of Automotive Engineering,Dalian University of Technology,Dalian,116024,China

2State Key Laboratory of Structural Analysis for Industrial Equipment,Dalian University of Technology,Dalian,116024,China

ABSTRACT Peridynamics(PD)is a non-local mechanics theory that overcomes the limitations of classical continuum mechanics(CCM)in predicting the initiation and propagation of cracks.However,the calculation efficiency of PD models is generally lower than that of the traditional finite element method(FEM).Structural idealization can greatly improve the calculation efficiency of PD models for complex structures.This study presents a PD shell model based on the micro-beam bond via the homogenization assumption.First,the deformations of each endpoint of the micro-beam bond are calculated through the interpolation method.Second,the micro-potential energy of the axial,torsional,and bending deformations of the bond can be established from the deformations of endpoints.Finally,the micro moduli of the shell model can be obtained via the equivalence principle of strain energy density(SED).In addition,a new fracture criterion based on the SED of the micro-beam bond is adopted for crack simulation.Numerical examples of crack propagation are provided,and the results demonstrate the effectiveness of the proposed PD shell model.

KEYWORDS Peridynamics;shell;micro-beam bond;interpolation method;crack propagation;homogenization

1 Introduction

Peridynamic (PD) theory was proposed by Silling [1] as a new formulation of the classical continuum mechanics (CCM) to solve fracture and failure issues.The PD equations of motion appear as integral-differential equations [2],and the integral does not depend on the continuity of the displacement.Thus,PD theory is spontaneously applied to analyze and solve discontinuity problems,such as crack propagation.In recent years,PD theory has been rapidly developed and applied to fracture or damage analysis of solid structures,such as the fracture problems of isotropic and anisotropic materials[3–5],the problems of plasticity[6–8],the crack branch[9–11]and impact[12] of brittle materials,composite materials [13–15],crack bending problem [16],and nonlinear analysis[17,18].

However,the calculation efficiency of PD models is usually lower than that of other traditional numerical methods,such as the finite element method(FEM),especially for complex structures.The coupling method is a strategy to improve calculation efficiency.Han et al.[19]proposed a local/nonlocal coupling model via the Arlequin approach;this model couples energy equations of the PD and CCM models in an overlapping domain.Han et al.[20]presented a hybrid adaptive local-continuum damage/PD model by using the morphing method[21,22].Structural idealization is often employed in engineering as a useful method to reduce the calculation time of PD models.Silling et al.[23]proposed a two-dimensional(2D)PD model for membranes and used this model for the numerical simulation of membrane stretching and dynamic tearing.Gerstle et al.[24]proposed the micropolar PD model that has two stiffness parameters.This model is mainly used to solve 2D plane problems.O’Grady et al.[25]introduced a non-ordinary state-based PD(SBPD)Kirchhoff-Love plate model based on beam models to represent the bending deformation of the plate,which is derived from the concept of a rotating spring with the bond between two material points.Taylor et al.[26]presented a 2D bond-based PD(BBPD)plate model through the asymptotic analysis of thin plates.This model can effectively analyze thin plate structures,including fractures caused by bending.Diyaroglu et al.[27] proposed a new PD Mindlin plate model based on the BBPD formulation while accounting for transverse shear deformation.The PD equations of motion were established and simplified to the same form as those of the classical Mindlin plate equations.Chowdhury et al.[28]proposed an SBPD model for linear elastic shells by introducing a shear correction factor and a shear deformation correction considering shell thickness.Ta¸stan et al.[29] developed a PD model for orthotropic Kirchhoff thin plates based on the BBPD formulation by using the homogenized approach.Sarego et al.[30]proposed a 2D linearized ordinary SBPD model for elastic deformation and computed the stiffness matrix for 2D plane stress/strain conditions,which is suitable for constructing membrane structures.Yang et al.[31]implemented the PD Mindlin plate formulations presented by Diyaroglu et al.[27] in the FEM framework.Nguyen et al.[32]developed a new ordinary SBPD model to predict the thermomechanical behaviors of 3D shell structures with six degrees of freedom,providing numerical methods and computational steps to deal with complex shell structures.Yang et al.[33]proposed a new PD Kirchhoff plate model based on the SBPD formulation,which is computationally efficient because each material point of the plate model has only one degree of freedom.Yolum et al.[34]presented a PD formulation for orthotropic Mindlin plates with transverse shear deformation.Vazic et al.[35]established a PD model for Mindlin plates based on Winkler’s elastic theory,which is an extension of the Mindlin plate model proposed by Diyaroglu et al.[27].Oterkus et al.[36]developed a new PD formulation for shell membranes.PD equations of motion were established by using Euler-Lagrange equations,and the bond constant of the shell membrane was determined by comparing PD equations of motion with classical equations of motion under a special condition that horizon tends to zero.Zhang et al.[37]developed a SBPD theory for nonlinear Reissener-Mindlin shells to model and predict large deformations of thick-walled shell structures.The controlling equations of the SBPD shell theory were derived by the non-local balance laws and the kinematic assumptions of Reissner and Mindlin plate and shell theories.Yang et al.[38]proposed a new SBPD model for functionally graded Kirchhoff plates and obtained the equations of motion of the proposed model by using the Euler-Lagrange equations and Taylor expansion.Zhang et al.[39] developed a non-local stress resultant-based geometrically-exact shell theory,which can simulate fracture and crack extension in shell structures under finite deformations.Shen et al.[40,41]proposed a PD shell model based on the micro-beam bond via the interpolation method,which is the foundation of this study.

The structural idealization strategies for PD plate and shell models proposed above can be divided into two types:BBPD model and SBPD model.Compared with the SBPD model,the BBPD model has the advantages of easy programming,simple calculation,and high calculation efficiency and has a wide range of application prospects.In the study,the in-plane deformation part of the PD shell in our previous work[40]is improved,and a new fracture criterion based on the strain energy density(SED)of the micro-beam bond is proposed.The reconstruction process of the whole PD shell model is described in detail through the strain homogenization method.First,the deformations of each endpoint of the micro-beam bond are calculated by the interpolation method.Second,the micro-potential energy of the axial,torsional,and bending deformations of the bond can be established from the deformations of endpoints.Finally,the micro moduli of the shell model can be obtained via the equivalence principle of SED.

The rest of this paper is organized as follows.Section 2 reviews the basic formulations of the microbeam bond with axial,torsional,and bending deformations.Section 3 mainly describes the derivation process of the micro moduli of various deformations in the PD shell model and presents a new energy failure criterion of the micro-beam bond for PD shells.In Section 4,some numerical examples are given,including a V-shaped notch plate under pressure loading,a notched plate under stretching,a notched cylinder under internal pressure,and a notched S-shaped curved shell under stretching.Finally,Section 5 provides the conclusions.

2 Micro-Beam Bond

As shown in Fig.1,the micro-beam bond is assumed to be composed of point xiand its adjacent point xj,which is located in the neighborhood horizonHxi.The neighborhood horizon of point xiis a circle of radiusδcentered on it.For the local coordinate system of the neighborhood horizonHxi,the normal of the horizon is defined as theZ-axis,and theX-axis and theY-axis lie in the plane of the horizon and are orthogonal,respectively.For the local coordinate system of the micro-beam bond,thex-axis is defined as the vector from point xito point xj,thez-axis is in the same direction as theZ-axis,and they-axis is perpendicular to thex-axis.Each micro-beam bond is subjected not only to axial deformation but also to torsional and bending deformations.The micro moduli of various deformations for the micro-beam bond model can be solved independently.

Figure 1:Micro-beam bond with its(a)horizon and(b)local coordinate system

In the local coordinate system of the micro-beam bond,the relationship between the force/moment densities and displacements of endpoints of the micro-beam bond can be expressed as

F=Ku

where

in which B and DPDcan be expressed as follows

and

In Eqs.(1)–(3),χis the shape parameter whose value ranges from 0 to 1;ξis the length of the micro-beam bond;andcax,cby,cbz,andctorare the axial deformation parameter,the bending deformation parameter around they-axis,the bending deformation parameter around thez-axis,and the torsional deformation parameter,respectively.

3 PD Shell Model

As shown in Fig.2,the horizonHxiof point xiis assumed to be exactly enclosed by a square element with nodes 1,2,3,and 4,and the length of the square element isL=2δ.For the local coordinate systemXYZof the square element,theX-axis is defined as the vector from node 1 to node 2,theY-axis is the vector from node 2 to node 3,which is perpendicular to theX-axis,and theZ-axis is the normal of the square element.For the convenience of calculation,the homogenization assumption for the region within the neighborhood horizon of point xican be extended to any region within the square element,considering that the overlapping region between the neighborhood horizon of point xiand the square element occupies the majority of the square element.

Figure 2:Horizon of point xi is exactly enclosed by a square element

3.1 Micro Moduli of in-Plane Deformation for PD Shells

The in-plane deformation only includes the translational displacement of theX-axis and theY-axis and the rotation angle displacement around theZ-axis.Therefore,in the coordinate system of the square element,the translational displacement along theX-axis and theY-axis and the rotation angle around theZ-axis of each node can be defined asUI,VI,andΘIZ,respectively.The in-plane displacement vector of the square element can be expressed as

where

The in-plane displacement vector at any position in the square element can be expressed by the Allman interpolation function[42]and Eq.(4)as

According to Eq.(5),the in-plane displacement vector at endpoints xiand xjof the micro-beam bond in the local coordinate systemXYZcan be obtained as

In the local coordinate system of the micro-beam bond,the in-plane displacement vector at endpoints xiand xjof the micro-beam bond is defined as

Through coordinate transformation,the in-plane displacement vector at endpoints xiand xjof the micro-beam bond in its local coordinate system can be obtained as

whereθis the angle between thex-axis and theX-axis.

As shown in Fig.3,the micro-beam bond can only be subjected to axial deformation along thex-axis and transversal deformation around thez-axis for the in-plane deformation.Thus,the strain of the micro-beam bond can be expressed in the local coordinate system of the micro-beam bond as

whereεx,κzrepresent the axial strain and the curvature in thez-axis direction,respectively,and

whereχis the shape parameter,andξis the length of the bond.

Figure 3:Micro-beam bond subjected to axial and transversal deformation for the in-plane deformation(a)undeformed state and(b)deformed state

Substituting Eq.(6)into Eq.(7)can be obtained as

According to the homogenization assumption,the strain in the square element is assumed to be equal everywhere and is=[εX εY γXY]T.Thus,the in-plane displacement vector of the square element can be defined as

whereεX,εY,andγXYrepresent the positive strain in theX-axis direction,the positive strain in theY-axis direction,and the shear strain in theX-axis and theY-axis directions,respectively.

Substituting Eq.(9)into Eq.(8)can be obtained as

where

The micro moduli matrix of the micro-beam bond corresponding to axial and transversal deformations is defined as

Therefore,the micro-potential energy of the micro-beam bond can be expressed by Eqs.(10)and(11)as

Then,the in-plane SED of point xiin PD theory can be expressed by Eq.(12)as

The in-plane SED of point xiin CCM can be expressed as

where

whereEis Young’s modulus,andνis Poisson’s ratio.

By setting Eq.(13)equal to Eq.(14),the values ofcaxandcbzcan be obtained as

3.2 Micro Moduli of Bending Deformation for PD Shells

The bending deformation only includes the transversal displacement along theZ-axis and the rotational displacement around theX-axis and theY-axis.Therefore,in the coordinate system of the square element,the transversal displacement along theZ-axis and the rotation angles around theX-axis and theY-axis of each node can be defined asWI,ΘIX,andΘIY,respectively.The bending displacement vector of the square element can be expressed as

where

The bending displacement vector at any position of the square element can be expressed by a cubic interpolation function[40]and Eq.(16)as

According to Eq.(17),the bending displacement vector at endpoints xiand xjof the micro-beam bond in the local coordinate systemXYZcan be obtained as

In the local coordinate system of the micro-beam bond,the bending displacement vector at endpoints xiand xjof the micro-beam bond is defined as

Through coordinate transformation,the bending displacement vector at endpoints xiand xjof the micro-beam bond in its local coordinate system can be obtained as

As shown in Fig.4,the micro-beam bond can only be subjected to transversal deformation around they-axis and torsional deformation around thex-axis.Thus,the strain of the micro-beam bond can be expressed in the local coordinate system of the micro-beam bond as

whereκy,φxrepresent the curvature in they-axis direction and the torsional curvature,respectively,and

Figure 4:Micro-beam bond subjected to transversal and torsional deformation for the bending deformation(a)undeformed state and(b)deformed state

Substituting Eq.(18)into Eq.(19)can be obtained as

According to the homogenization assumption,the strain in the square element is assumed to be equal everywhere and is~εXY=[κY κX κXY]T.Thus,the bending displacement vector of the square element can be defined as

whereκX,κY,andκXYrepresent the curvature in theX-axis direction,the curvature in theY-axis direction,and the torsional curvature in theX-axis and theY-axis directions,respectively.

Substituting Eq.(21)into Eq.(20)can be obtained as

where

The micro moduli matrix of the micro-beam bond corresponding to bending and torsional deformations is defined as

Therefore,the micro-potential energy of the micro-beam bond can be expressed by Eqs.(22)and(23)as

Then,the bending SED of point xiin PD theory can be expressed by Eq.(24)as

The bending SED of point xiin CCM can be expressed as

where

By setting Eq.(25)equal to Eq.(26),the values ofcbyandctorcan be obtained as

So far,the four micro moduli of various deformations for the PD shell model in Eqs.(15)and(27)have all been solved,and the stiffness matrix of a micro-beam bond for PD shells can be formed easily.

3.3 Failure Criterion

By comparing the relationship between the elongation ratesand the critical elongation rates0of the bond between two material points,this single judgment criterion cannot accurately describe the damage degree of the micro-beam bond.In this study,the fracture criterion based on the idea of strain energy is adopted to simulate the cracks.The proposed PD shell model belongs to the 2D condition.As shown in Fig.5,similar to the 3D condition,the fracture energy[43,44],i.e.,the energy required to open a new fracture surface per unit area,can be determined by the energy released by breaking all bonds across the unit fracture surface.Assuming that the critical micro-potential energy of each bond isωc,the integral equation of the fracture energyGcfor the 2D condition can be defined as

wheretis the thickness of the shell andωccan be expressed as

where

whereεx,κy,κz,andφxrepresent the axial strain,the curvature in they-axis and thez-axis directions and the torsional curvature of the micro-beam bond,respectively.

Figure 5:Schematic of the fracture surface

Assuming that the micro-beam bond is uniformly deformed in the neighborhood horizon of point xi,all strains in the neighborhood horizon are equal.Therefore,Eq.(29)can be rewritten as

whereξis the length of the micro-beam bond.

Then,substituting Eq.(30)into Eq.(28),the integral fracture energy equation can be calculated as

Finally,the energy failure criterion of a micro-beam bond can be expressed by Eq.(31)as

4 Numerical Examples

4.1 V-Shaped Notch Plate under Pressure Loading

A soda-lime brittle glass plate is selected to demonstrate that the proposed shell model can simulate its crack propagation.The geometry and dimensions of the glass specimen are shown in Fig.6.Thec sample has the dimension of widthW=100 mm,heightH=150 mm,and thicknesst=4.65 mm.There is a V-shaped notch with an angle of 40°on the left side of the sample,and a pre-crack exists at the end of the V-shaped notch.The material properties of the specimen are taken from the literature[45],with Young’s modulus ofE=72 GPa,Poisson’s ratio ofν=0.23,density ofρ=2440 kg/m3,and fracture energy ofGc=3.8 J/m2.In the experiment[46],the impact effect on the V-shaped notch is caused by an impactor.Instead of simulating the impactor itself,the notch pressure load is chosen to be estimated in the computational tests to reduce the computational cost of the model.The loading conditions for the notch pressure load are taken from the literature[45],where it is assumed that the load is uniformly distributed along the loading area.Only one layer of material points of the V-shaped notch plate is discretized,with a discretization spacing ofΔ=0.5 mm.According to suggestions of the literature[47],the selected neighborhood horizon of the material points isδ=3Δ,and the time increment step isΔt=1×10-7s,which is sufficient to satisfy the stability condition of numerical integration.In this study,a precise volume correction algorithm[48]for 2D PD models is introduced to improve the computational accuracy of numerical integration.

Figure 6:Geometry and dimensions of the V-shaped notch plate

Numerical results are shown in Figs.7a–7c.As the loading time is prolonged,the crack propagates slowly in the pre-crack direction and the crack bifurcation appears at approximately 60μs,at a position 60 mm away from the left side of the sample,with a bifurcation angle of around 45°.Two bifurcated crack propagation paths can be observed at the front end of the plate.Experimental results provided in [46] are shown in Fig.7d.Numerical results in Fig.7c are in good agreement with experimental results of fractured specimens in Fig.7d,which shows that the PD shell model and fracture criterion proposed in this study can simulate the crack propagation path very accurately.

Figure 7:Numerical results of crack propagation for the V-shaped notch plate at loading time of(a)45 μs,(b) 81 μs,and (c) 112 μs;(d) experimental photographs of fractured specimens,reproduced with the permission of[46],copyright@Journal of the Mechanics and Physics of Solids,2018

4.2 Pre-Notched Plate under Stretching

As shown in Fig.8,the material characteristics,geometric shape,and boundary conditions of the rectangular thin plate are selected as described in literature [49].The rectangular thin plate has the dimension of lengthL=100 mm,widthW=40 mm,and thicknesst=1 mm.The pre-notch is located in the middle of the left side of the rectangular thin plate,whose length isl=50 mm.The Young’s modulus,density,Poisson’s ratio,and fracture energy of the thin plate areE=72 GPa,ρ=2440 kg/m3,ν=1/3,andGc=135 J/m2,respectively.The tensile loading with equal magnitude and the opposite direction is exerted to the upper and lower edges of the thin plate,whose magnitude is kept constant asσ=12 MPa throughout the simulation.The neighborhood horizon of the material point isδ=3Δ,and the same fracture criterion as in the previous example is selected.

Figure 8:Dimensions and boundary conditions of the pre-notched plate

Four discrete mesh sizes of 1.25,1,0.8,and 0.5 mm are selected to investigate the effect of mesh refinement on cracks,respectively.The crack propagation paths of the plate using four uniform grids at 47μs are shown in Fig.9.The tensile loading on the upper and lower edges of the thin plate increases abruptly and then remains constant after reaching the required value.As the loading time is prolonged,the crack propagation path slowly expands along the direction(X-axis)of the pre-crack until near 21μs,crack bifurcation occurs,and two bifurcated crack propagation paths can be observed at the front end of the plate.The crack propagation paths are highly similar for all four models.As the horizon decreases,the damage region becomes smaller,i.e.,the damage region becomes more precise.The fine mesh captures the degree of damage more accurately than the coarse mesh.Numerical results show that the crack propagation paths predicted by the proposed PD shell model are roughly similar to the numerical results in the references[48,49].

Figure 9:Crack propagation paths of the pre-notched plane under different discretizations.Mesh spacing Δ of(a)1.25 mm,(b)1.0 mm,(c)0.8 mm,and(d)0.5 mm

4.3 Pre-Notched Cylinder under Internal Pressure

A cylindrical shell[50]subjected to internal pressure in the axial direction and symmetrical prenotches with bending deformation and membrane deformation is selected to verify the applicability of the proposed PD shell model to the crack propagation of the ordinary shell model.As shown in Fig.10,the cylinder has a cross-sectional radius ofR=5 mm,an axial length ofL=30 mm,and a thickness oft=0.125 mm.Two symmetrical pre-cracks are located in the middle of the cylindrical shell axially(Z-axis),with a single crack length ofl=4 mm.Two cases (a) without and (b) with tear straps are considered.The purpose of installing the tear straps is to restrain the axially propagating crack.The mechanical parameters of the cylinder are taken from literature[50],where Young’s modulus isE=70 GPa,Poisson’s ratio isν=0.3,and the fracture energy isGc=1.5 J/m2,respectively.The displacements in theX-axis and theY-axis directions at both ends of the cylinder are constrained.Exploiting symmetry,only one-half of the geometry of the cylinder is modeled,and only one layer of material points of the thin cylinder is discretized with a discretization spacingΔof 0.25 mm.The selected neighborhood horizon isδ=3Δ,and the fracture criterion is as described in Section 3.3 above.The damage clouds of the cylinder in case (a) and case (b) at different force balance stages are shown in Figs.11 and 12.For the cylinder without tear straps,Fig.11 shows that a straight crack parallel to the axial direction(Z-axis)of the cylinder propagates toward its ends under the action of internal pressure.As expected,a straight crack propagation path is observed.For the cylinder with tear straps,Fig.12 shows that the crack propagates in the axial direction in the same way as case(a)until the crack approaches the tear straps,where the direction of crack propagation suddenly changes from axial to circumferential,and the crack branching occurs in the circumferential direction.The presence of tear straps significantly affects the crack propagation path.The damage numerical results of these two cases are the same as those in the literature[50].Thus,this phenomenon is successfully simulated by the proposed PD shell model.

Figure 10:Axially pre-notched thin cylinders under internal pressure(a)without tear straps and(b)with tear straps

Figure 11:Crack propagation path of the axially pre-notched cylinder without tear straps under internal pressure (a) 0.81 MPa,(b) 0.49 MPa,(c) 0.35 MPa,(d) 0.24 MPa,(e) 0.16 MPa,and (f)0.095 MPa

Figure 12:Crack propagation path of the axially pre-notched cylinder with tear straps under internal pressure(a)0.85 MPa,(b)0.50 MPa,(c)0.40 MPa,(d)0.47 MPa,(e)0.53 MPa,and(f)0.46 MPa

4.4 Pre-Notched S-Shaped Curved Shell under Stretching

In this example,an S-shaped curved shell with pre-notches is selected to verify that the proposed PD shell model is applicable to any shell structure and to deal with thin shells with arbitrary cracks.The geometry and dimensions of the long side of the S-shaped curved shell are shown in Fig.13,which is formed by two arcs tangent to each other with a radius of R=80 mm.An external load of equal magnitude in opposite directions is exerted on the entire long edge of the shell with a load magnitude ofσ=12 MPa.Two cases of pre-notches with an angle of(a)0°and(b)45°are considered,as shown in Fig.13,to observe the crack propagation.The length of the short side of the shell isL=50 mm,and a pre-notch with a length ofl=10 mm exists at the tangent in the middle of the two circular curved surfaces of the thin shell,which is uniformly distributed along the tangent of the curved surface.The material parameters of the thin shell are as follows:Young’s modulus ofE=72 GPa,Poisson’s ratio ofν=1/3,thickness oft=1 mm,and fracture energy ofGc=80 J/m2.Only one layer of material points of the S-shaped surface is discretized,with a discretization spacing ofΔ=1 mm.The neighborhood horizon of the material point isδ=3Δ,and the same fracture criterion as in the previous examples is selected.

Figure 13:S-shaped shell structure with pre-notch(a)parallel to YZ-plane and(b)45°with YZ-plane

The damage diagrams of the shell structure at different loading times for two cases are shown in Figs.14 and 15.Fig.14 shows that the crack parallel to theYZ-plane propagates in two opposite directions,and the propagation path is roughly parallel to theYZ-plane.However,at the end of the crack propagation,crack bifurcation occurs,and two slightly bifurcated crack propagation paths can be observed.As shown in Fig.15,the crack at an angle of 45°with theYZ-plane first propagates slightly on the shell surface along the pre-notched direction,the direction of crack propagation path changes,and the crack propagates roughly along theYZ-plane in two opposite directions.When the loading time is approximately 21μs,the first crack bifurcation occurs,and the new crack propagates roughly along theX-axis in two opposite directions.The original and bifurcated cracks propagate simultaneously in the established directions.Eventually,the bifurcation propagation path of the crack can be observed very clearly.

Figure 14:Propagation of crack parallel to the YZ-plane,and the time of(a)3 μs,(b)10 μs,(c)17 μs,(d)24 μs,(e)31 μs,and(f)38 μs

Figure 15:Propagation of the crack at an angle of 45°with the YZ-plane,and the time of(a)3 μs,(b)10 μs,(c)17 μs,(d)24 μs,(e)31 μs,and(f)38 μs

5 Conclusion

In this study,a PD shell model based on the micro-beam bond via the homogenization method is presented.First,the deformations of each endpoint of the micro-beam bond are calculated using the interpolation method.Second,the micro-potential energy of the axial,torsional,and bending deformations of the bond can be established from the deformations of endpoints.Finally,the micro moduli of the shell model can be obtained via the equivalence principle of SED.The most critical aspect is the establishment of the conversion relationship between strains of the micro-beam bond and strains of the square element,which can significantly simplify the process of deriving the micro moduli.Different deformations of the micro-beam bond are considered;thus,it is not restricted by material parameters.On the basis of elasticity and small deformations,the 2D plane stress model and the bending model are combined to form a general PD shell model,which is a combination of axial,torsional,and bending deformations.The effectiveness of the proposed PD shell model is demonstrated by several numerical examples.The strain homogenization method used to construct the PD shell model in this study will be applied to the construction of other models in the future,such as anisotropic and large deformation models.

Acknowledgement: The authors wish to express their appreciation to the reviewers for their helpful suggestions which greatly improved the presentation of this paper.

Funding Statement:This work was funded by Project of the National Natural Science Foundation of China(Grant No.11872017).The support is gratefully acknowledged.

Conflicts of Interest:The authors declare that they have no conflicts of interest to report regarding the present study.

Appendix A.Shape functions for in-plane deformation of the PD shell model

The parameters of the bottom right corner mentioned in Eq.(A1)can be defined as

where MOD(x,y)represents the remainder ofxdivided byy.

XcandYcare the coordinates of the square element centroid.The components of(p,q)can be expressed as

whereK,I,Jis an ordered triplet,and the order of value is(5,1,2),(6,2,3),(7,3,4),(8,4,1).

Appendix B.Shape functions for bending deformation of the PD shell model

whereXcandYcare the coordinates of the square element centroid.The components of(p,q)in Eq.(B1)can be expressed as