Large deformation and wrinkling analyses of bimodular structures and membranes based on a peridynamic computational framework

H. Li·Y. G. Zheng·Y. X. Zhang·H. F. Ye·H. W. Zhang

Abstract In this paper, the quasi-static large deformation, wrinkling and fracture behaviors of bimodular structures and membranes are studied with an implicit bond-based peridynamic computational framework.Firstly,the constant and tangential stiffness matrices of the implicit peridynamic formulations for the nonlinear problems are derived, respectively. The former is constructed from the linearization of the bond strain on the basis of the geometric approximation while the latter is established according to the linearization of the pairwise force by using first-order Taylor’s expansion.Then,a bimodular material model in peridynamics is developed, in which the tensile or compressive behavior of the material at each point is conveniently described by the tensile or compressive states of the bonds in its neighborhood.Moreover,the bimodular material model is extended to deal with the wrinkling and fracture problems of membranes by setting the compressive micro-modulus to be zero.In addition,the incremental-iterative algorithm is adopted to obtain the convergent solutions of the nonlinear problems.Finally,several representative numerical examples are presented and the results demonstrate the accuracy and efficiency of the proposed method for the large deformation,wrinkling and fracture analyses of bimodular structures and membranes.

Keywords Bimodular structures·Wrinkling membranes·Fracture problems·Peridynamics·Implicit numerical algorithm

1 Introduction

The asymmetric properties under tension and compression are ubiquitous in many engineering materials, such as ceramics,polymers,rubber,and composite materials.These materials often exhibit different mechanical behaviors in tension and compression because they have different tensile and compressive moduli[1–4].This type of material is called bimodular material and has been extensively used in the aerospace engineering area due to the superior mechanical characteristics arising from its asymmetric mechanical behaviors.For instance,tensegrity structures,which are the self-stressed structures existing prestresses in the cables and struts to maintain initial geometrical shapes,have been widely utilized in space engineering because they are flexible and light weight [5]. Thin membranes, a special type of bimodular material by removing the corresponding compressive resistance of the materials (i.e., the compressive modulus in the bimodular material model is set to be zero),can endure large tensile loads,whereas it is wrinkled when a small compressive or shear load is applied[6,7].These have been extensively used as foldable components in aerospace engineering.Because of these special features,it is of great importance to deeply understand the deformation and actuating mechanisms of bimodular materials for their practical engineering applications.

During the past several decades,the mechanical analyses of bimodular materials have evolved into a mature subject and extensively related studies have been conducted from three aspects.From the theoretical aspect,Jones[8],Vijayakumar and Rao[9]and Liu et al.[10]discussed the stress–strain relations of the general materials, the composites and the NiTi shape memory alloys with different moduli in tension and compression,respectively.Then,Yao and Ye[11]performed analytical studies on the bending responses of columns and beams with different tensile and compressive moduli.From the experimental side, Blandino et al. [12] and Wong and Pellegrino [13] investigated the wrinkling phenomena of a square membrane under a symmetrical load and a rectangle membrane under a shear load.Moreover,from a numerical simulation aspect, the analyses of the wrinkled membranes were conducted by adopting many approaches, such as the finite element method [14–22], a series of variational principles and the related bounding theorems[23,24],the finite element method with parametric variational principle and co-rotational approach[25–28],the peridynamic method of trusses and tensegrity structures [29]. Among these works,the peridynamic method has many natural advantages to deal with the geometrical nonlinear problems and it becomes a great promising method to simulate the bimodular materials with large deformation.

The peridynamic theory proposed by Silling[30]is a nonlocal theory and is governed by an integral equation. Since the proposal of the peridynamics, the fundamental theories of the peridynamics have been gradually developed from different perspectives in recent years[31–38].Meanwhile,the peridynamic method has been widely applied to the numerical simulation in various fields. For instance, a morphing approach was developed by Han et al. [39, 40] to analyze damage problems based on the peridynamics.An improved peridynamic approach was presented by Huang et al. [41]for the brittle fracture analysis. The peridynamic formulation for the transient heat conduction in bodies involving discontinuities was developed by Bobaru and Duangpanya[42,43].The peridynamic model of thermal diffusion problem was proposed by Oterkus et al.[44].The 3D peridynamic model for pitting corrosion damage was investigated by Jafarzadeh et al. [45]. The coupled peridynamic model for porous flow in porous media were simulated by Katiyar et al.[46],Ouchi et al.[47]and Lai et al.[48].The peridynamic simulation for fatigue life and residual strength prediction of composite laminates was presented by Hu and Madenci[49]. In addition, the peridynamic method was applied to perform numerical simulations of the stretching,wrinkling,rupture, and dynamic tearing of membranes and biomembranes [50–52]. These works show that the peridynamic method is an effective route to deal with the problems of both the geometrical and material nonlinear characteristics.However, so far there have been few studies to address the peridynamic model for the nonlinear mechanical analysis of bimodular structures and membranes. Moreover, the computational efficiency is low when the explicit peridynamic formulation is used to analyze the quasi-static problems and the convergent behaviors are unsatisfied when the linearized theory [33] is applied to deal with the nonlinear problems.

Therefore, to improve the computational efficiency and generate good convergent behaviors, an efficient implicit peridynamic computational framework for the quasi-static large deformation and wrinkling analyses of bimodular structures and membranes is presented in this paper. Firstly, the constant and tangential stiffness matrices of the implicit peridynamic formulations for the nonlinear problems are derived, respectively. The former is constructed from the linearization of the bond strain on the basis of the geometric approximation, that is,the bond strain is considered the linear function of the displacements of the material points approximately and the current displacements of the material points are decomposed as the sum of the initial displacements and the incremental displacements. Otherwise, the latter is established according to the linearization of the pairwise force by using first-order Taylor’s expansion. As far as the stiffness matrix of each bond is obtained, the global stiffness matrix can be assembled according to the stiffness matrix of each bond. Then, the bimodular material model in peridynamics is proposed in this paper. The tensile or compressive state of the material at each material point is described on the basis of the corresponding state of the bonds in its neighborhood. Thus, we only need to determine the state of each bond in each iteration step and to update the corresponding stiffness matrix of each bond in the simulating processes.Meanwhile,the bimodular material model is also used to analyze the wrinkling and fracture problems of the membranes by setting the compressive micro-modulus to be zero. Moreover, load and displacement control based incremental-iterative methods are utilized to solve the nonlinear problems. It is worth mentioning that the prescribed boundary and loading conditions can be accurately imposed on the actual material points in the present implicit peridynamic method.

The remaining sections of this paper are organized as follows. Section 2 briefly describes the fundamental concepts of the bond-based peridynamic formulation.Section 3 shows the derivations of the constant and tangential stiffness matrices for the implicit peridynamic formulations in detail.The bimodular material and membrane models in peridynamics are presented in Sect.4.The specific computational procedure of the present method for solving the nonlinear problems involving material and geometrical nonlinearities is then given in Sect. 5. In Sect. 6, several representative numerical examples are presented and the results verify the accuracy and efficiency of the proposed method by comparing it with the reference solutions.Finally,some conclusions are given in Sect.7.

2 Fundamental bond-based peridynamic formulation

The fundamental integral equation of motion for any material point based on the bond-based peridynamic model is given by[30]

The pairwise force f satisfies the basic mechanical considerations of the conservation oflinear and angular momenta and can be expressed as

in which ξ = xj-xiis the relative position of these two material points in the initial configuration and called a bond,is the relative displacement at time t,c(‖ξ‖)is the micro-modulus that represents the elastic stiffness of the corresponding bond,is the unit vector along the direction of the relative position in the current configuration.The bond strain s is defined as

It is noted that the micro-modulus c(‖ξ‖)in this paper is obtained based on the equivalence between the strain energy densities of the peridynamic and the classical continuum models.Moreover,to reduce the influence of the geometric boundary on the computational accuracy,the micro-modulus c(‖ξ‖) is calculated numerically according to the actual neighborhood for every discretized material point[41],that is

where E is the Young’s modulus, ν is the Poisson’s ratio,Nxiis the number of material points among the neighborhood of material point xiand Vxjis the volume of material point xj.

3 Implicit incremental formulations

In the present work, the quasi-static problem under large deformation with small strain situation is considered. Usually, the external loads are applied incrementally for this nonlinear problem. It is assumed that the system at the nth load step is in equilibrium and the displacementis known.Then,based on the fundamental equation of motion(i.e.,Eq.(1)),the static equilibrium equation at the(n+1)-th load step can be expressed as

The meshless method[31]can be used to solve the above equation. The domain of interesting is discretized into N material points and then the equilibrium equation corresponding to each material point xican be written as

3.1 Constant stiffness matrix method:linearization of the bond strain via geometric approximation

Based on the geometric relationship of the bonds as shown in Fig.1(for a two-dimensional case),the bond strain given by Eq.(3)can be rewritten as a linear function of displacements.

Fig.1 Schematic diagrams of the bond with a initial,b deformed,c current configurations

Then the corresponding bond strain for two-dimensional problems can be calculated by

It can be seen from Eq. (9) that the pairwise force is a nonlinear function of displacement,and thus the equilibrium Eq.(6)is a nonlinear equation.To solve this problem,Eq.(9)should be linearized.Since the deformation with small rotation in each load step is considered here,dθ is small enough in each iteration step, and cos(θn+dθ) ≈cos θnand sin(θn+dθ) ≈sin θn. Then the linearized pairwise force can be written as

For simplicity,the pairwise force can be rearranged to

Substituting Eq. (11) into Eq. (6), the discretized incremental formulation of the static equilibrium equation for material point xibased on the initial configuration can be expressed as

with

3.2 Tangential stiffness matrix method:linearization of the pairwise force via Taylor’s expansion

On the other hand,to obtain the general nonlinear algebraic equations of geometrical nonlinear problems,Eq.(6)can be approximated as follows according to the first-order Taylor’s expansions,that is

where δη is the same as the incremental displacements δdij.

According to Eqs. (11), (16) and (17), Eq. (18) can be written as

In addition,based on Eqs.(2)and(3),the derivative of the pairwise force fijcan be expressed as

Then, Eq. (19) can be expressed in the form of Eq. (15)and the corresponding stiffness matrixcan be rewritten as

Comparing with Eq. (16), it can be seen that Eq. (21) is variable in pace with the deformation of the structure,and is related to the current configuration of the deformed structure and has an additional term which represents the “geometric”or“nonlinear”contribution to the stiffness matrix.Thus,Eq. (21) can be called as the tangential stiffness matrix.It is noted that the form of the tangential stiffness matrixis similar to that of the tangential stiffness matrix of a truss element on condition of geometrically nonlinear analysis. Furthermore, Eq. (21) can be reduced to Eq. (16) if‖ξ+η‖→‖ξ‖on condition of small deformation.

3.3 Assembling of the global incremental equation

Finally, the discrete equation for the whole domain can be obtained by assembling Eq.(15)for each material point and can be written as

where Kniis a 2 by N×2 matrix that denotes a portion of Knand can be assembled by knij, A represents an assembling operator that is similar to the operator in the conventional finite element method[26].Moreover,ifcan be calculated by Eq.(16)or Eq.(21);otherwise,

Actually,it is more convenient to assemble the global stiffness matrix Knvia the stiffness matrix Knbcorresponding to each bond,that is

with

in which the b-th bond is formed by material points xiand xj.From Eqs.(16)and(21),it can be seen that ifthenandis symmetric for each bond,which will result in a symmetric and banded global stiffness matrix Kn.In addition,the total number of the bond Nbondcan be obtained by

It is worth noting that the prescribed boundary and loading conditions for the final discretized governing equation(i.e.,Eq.(22))are imposed on the corresponding discretized material points of the actual material region.

4 Bimodular material and membrane models in peridynamics

In this paper,the structures composed of bimodular materials are considered.As is known,the tensile or compressive state of two-dimensional bimodular materials is usually described by the stress state of a local point in the continuum mechanics model.However,in the present peridynamic model,the tensile or compressive state of the bimodular materials at every material points can be simply presented by the force–elongation relations of the bonds in their neighborhoods as shown in Fig. 2a,b, in which c+and c-are the micro-moduli that denote the bonds are under tensile and compressive states,respectively.Therefore,the corresponding pairwise force of a bond for the bimodular material together with crack propagation can be expressed as

Fig.2 Peridynamic constitutive models.a Bimodular materialb Bimodular material c+ ＜c-.c Membrane c- =0

where the micro-moduli c+(ξ) and c-(ξ) are related to the tensile and compressive Young’s moduli and Poisson’s ratios,respectively,is used to denote the failure state of the bond with s0the critical stretch refer to the critical failure state of the bond.

Consequently,according to Eq.(21),the tangential stiffness matrix of the tensile bond of the bimodular material can be written in the following form

and the tangential stiffness matrix of the compressive bond of the bimodular material can be written in the following form

In addition,the bimodular material model in peridynamics can be easily extended to deal with the wrinkling and fracture analyses of membranes with no compressive resistance by setting the compressive micro-modulus to be zero as shown in Fig.2c.Therefore,the corresponding pairwise force of a bond of the membrane can be expressed as

Moreover,to measure the wrinkled level of the deformed membrane,we define a scalar field of the wrinkled level as the fraction of the compressive bonds at a material point in its neighborhood. Thus, the scalar L of the wrinkled level can be expressed as

where sign(s) is a sign function and it is defined as sign.It can be seen that,withdenoting that the membrane at the point xiis in the tensile state in all directions andrepresenting that the membrane at the point xiis in the slack state in all directions.In addition,the virial stress of each material point is defined as[53]

5 Implementation of the algorithm

In this paper, the bimodular structures and membranes that experience large deformation and small strain are discussed.In other words,both the material and geometrical nonlinearities are considered in the present work.To solve the nonlinear equation, the external load is usually applied through several load steps and the load control or displacement control based incremental-iterative method[54]is utilized to obtain the convergent solution in each load step. Therefore, the linearized incremental formulation of Eq. (22) should be rewritten as

where the superscripts n and n+1 denote the n-and(n+1)-th iteration step,respectively,the superscript k+1 represents the(k+1)-th load step.

Then, the detailed computational procedure of the proposed numerical algorithm is listed as follows.

1. Input the material and geometrical parameters of the bimodular structures,discretize the system in space with material points and setup the other information,such as the horizon radius δ,the maximum number of the iteration step nmaxand the maximum number of the load step kmax,etc.

2. Loop for the load step(initialize k =0).

2.1. Set k =k+1.

2.2. Loop for the iteration step (initialize n = 0 and ΔUk+1=0).

2.2.1. Set n =n+1.

2.2.2. Form the stiffness matrixof the bond,

2.2.3. Form the global stiffness matrix Kn,k+1and the global internal force vector

2.2.4. Calculate the incremental displacement δUn+1,k+1of the current iteration step,update the total increment displacement ΔUk+1of the current load step according to ΔUk+1= ΔUk+1+ δUn+1,k+1and the total displacement Uk+1according to U k+1 =Uk+ΔUk+1.

2.2.5. Determine whether the solution satisfies the convergent condition er r ≤φ.If YES,finish the calculation of the current load step and go to step(3);else,loop for the iteration step and go to step(2.2.1)when n ＜nmax,or go to step(2.2)to recalculate the current load step by decreasing the assigned load increment when n = nmax(If the value of the load increment is smaller than a limitation,finish the computations).

3. If k =kmax,go to step(4);else,go to step(2.1).

4. Finish the computations and output the results.

The flow chart as shown in Fig. 3 depicts the detailed computational procedure mentioned above,in which the convergence criterion is defined as

where φ is the convergence tolerance and φ =1.0×10-5is used in the following simulation examples.

6 Numerical examples

In this section, two examples of bimodular structures, two examples of wrinkled membranes and one example of fractured membrane are taken into account to demonstrate the accuracy and efficiency of the proposed implicit peridynamic method for quasi-static large deformation, wrinkling and fracture analyses of bimodular structures and membranes.The solutions obtained by the explicit peridynamic method with an adaptive dynamic relaxation technique are taken as the reference solutions and are labeled as EXPD while the results obtained by the present method are labeled as IMPD.

Fig.3 Flow chart of the computational procedure

In addition,the errors of the displacement fields are defined as

where UIMPDis the solution obtained by the present method and URstands for the reference solutions.

6.1 Tensile bimodular column with self-weight

First,a tensile bimodular column with self-weight,which is equivalent to a one-dimensional analytical example is considered[27].The geometrical information and the boundary and loading conditions are depicted in Fig. 4a. The length of the column is L = 10 m and the cross sectional area is A=1 m2.The external load is P =6 Pa and the self-weight is γ = 2.0 Pa/m. Furthermore, the compressive Young’s modulus is fixed as E-= 5000 Pa and both the Poisson’s ratios under tension and compression are assumed to be zero.Figure 4b shows that the column can be divided into two parts that experience tension and compression,respectively,under the condition of considering both the external load and the self-weight.The analytical solution of the displacement along axis z can be found by[27]

Fig.4 Schematic diagrams of a bimodular bar with gravity

in which h = L - P/(γ A)is a critical value of tension and compression.The upper part(z ＞h)is in tension while the lower part(z ＜h)is in compression.Moreover,the column is discretized into 200 material points with the same grid size of Δz =0.05 m.The horizon of the material point is chose to be 3Δz.Three different types of c-/c+ratios(i.e.,c-/c+=1, 5, 10) are considered and all loads are applied by one load step for each case.It is emphasized that the tangential stiffness matrix is reduced to the constant stiffness matrix for the one-dimensional problem.Figure 5 shows the variations of the displacements along z direction for different kinds of c-/c+ratios, which demonstrates that the results obtained by the present method fit well with the analytical solutions.Further, Fig. 6 depicts the errors of the displacement fields obtained by the present method and the analytical method changing with ratio of c-/c+.It can be seen from the slope of the curve that the results are satisfied when the c-/c+ratio is very large.Consequently,the present method is valid and accurate for the mechanical analysis of the bimodular structures.

6.2 Bimodular cantilever beam with large deformation

A cantilever beam as shown in Fig. 7 is considered in this example.The sizes of the beam are 2 m×0.21 m×0.01 m.

Fig.5 Distributions of the vertical displacement along the bar for different types of c-/c+ ratios

Fig.6 L2 error of the displacement fields changing with the c-/c+ratio

Fig.7 Schematic diagram of a cantilever beam subjected to a concentrated force

Fig.8 Distributions of the horizontal displacement along the middle line obtained by the present method and the explicit peridynamic method

Fig.9 Distributions of the vertical displacement along the middle line obtained by the present method and the explicit peridynamic method

Table 1 Iteration steps at each load step with the constant stiffness matrix method for different types of c+/c- ratios

Table 2 Iteration steps at each load step with the tangential stiffness matrix method for different types of c+/c- ratios

Plain stress condition is assumed. For the boundary and loading conditions, the left side of the beam is fixed and a concentrated force F = 1000 N is imposed on point A.The tensile modulus is E+= 3.84×108Pa and the Poisson’s ratio is ν+= ν-= 1/3. The beam is discretized into 200 × 21 material points with the same grid size of Δx = Δy = 0.01 m. The horizon of the material point is δ = 3Δx. Four different types of c+/c-ratios (i.e.,c-/c+= 1, 5, 10, 20) are considered and all loads are applied by ten load steps for each case.Figures 8 and 9 show the distributions of the displacements along the middle line obtained by the present method(Here,the tangential stiffness matrix is adopted.)and the explicit peridynamic method with an adaptive dynamic relaxation technique.It can be seen from Figs.8 and 9 that the results obtained by the present method are in good agreement with those obtained by the explicit method.Tables 1 and 2 give the number of iteration at each load step for different types of c+/c-ratios when the constant and tangential stiffness matrices are applied, respectively,(i.e.,Eqs.(16)and(20)).Table 1 declares that the convergent errors do not satisfy the convergent criteria because they are oscillated periodically when the deformation of structure and the c+/c-ratios are large for the constant stiffness matrix method. However, Table 2 illustrates that the convergent behaviors are very good even the deformation of structure and the c+/c-ratios are large for the tangential stiffness matrix method. The reason for this is that the tangential stiffness matrix obtained by the Taylor’s expansion is on the basis of the current configuration and contains the “geometric”or“nonlinear”term.Therefore,only the tangential stiffness matrix(i.e.,Eq.(21))is used in the present method in the following simulations.In addition,the curves of the convergent behaviors as shown in Fig.10 investigate that the convergent solutions can be obtained by the explicit peridynamic method through the adaptive dynamic relaxation technique for these nonlinear problems.As for the computational efficiency,the present method and the explicit peridynamic method take 2.31 s and 0.01 s for each iteration step,respectively.However, the former only needs about 20 iteration steps while the latter method need about 9000–13,000 iteration steps to obtain the convergent solutions.Consequently,the computational efficiency of the present implicit peridynamic method is higher than that of the explicit peridynamic method.Moreover, Fig. 11 shows the deformed shapes of the beam with different c+/c-ratios,which represents that the deformation of the bimodular beam is asymmetric and neutral layer shifts up comparing with the classical beam.

6.3 Wrinkling of a square membrane

Fig.10 Convergence behaviors of the vertical displacement of the point A in the explicit peridynamic method for different types of c+/c-ratios

A classical example of the wrinkling analysis of a square membrane is considered. The square membrane is shown in Fig. 12a, which has been experimentally investigated by Blandino et al. [12]. Figure 12b shows the schematic of the numerical model of the square membrane that is discretized into 100×100 material points with the same grid size of Δx = Δy = 5.0×10-3m. The thickness of the membrane is 1.0 × 10-3m. In order to reduce the nonlocal effects, the horizon of the material point is chosen as δ =5.0×10-3m.The tensile Young’s modulus of the membrane is E+= 1×105Pa and the corresponding Poisson’s ratio is ν+= 1/3,while the compressive Young’s modulus E-and Poisson’s ratio ν-of the membrane are set to be zero.Four concentrated forces with F1= 5 N are symmetrically applied on the four corners of the membrane. Moreover, to prevent rigid body motion, the four material points located at the center of the square membrane are fixed.Figure 13a,b plots the x-directional and y-directional displacement contours of the membrane, respectively. Figure 14a, b depicts contours of the maximumprincipal stress and the corresponding magnitude/direction, respectively. Figure 15a,b shows the wrinkling patterns obtained by the present method and the finite element method in Ref.[12],respectively,in which the blue and white regions are wrinkled while the red and gray regions are taut.From these figures,it can be seen that the wrinkling phenomena of the square membrane obtained by the present method fit well with the reference solutions.In addition,it is noted that all of the external loads are applied by only one load step and the corresponding number of the iteration step is 6.

Fig.11 Deformed configurations of the bimodular beam contoured with the vertical displacement. a c+/c- = 1. b c+/c- = 5. c c+/c- = 10.d c+/c- =20

Fig.12 Square wrinkled membrane. a Experimental model [12].b Numerical model

Fig.13 Displacement contours of the square membrane obtained by the present method.a x-directional displacement.b y-directional displacement

Fig.14 Stress patterns.a Maximum principal stress.b Magnitude and direction of the maximum principal stress obtained by the present method

6.4 Rectangle Kapton membrane in shear

Fig.15 Wrinkling patterns. a Obtained by the present method.b Obtained by Blandino et al.[12]using the finite element method

Fig.16 Rectangle Kapton membrane under shear loading. a Experimental model[13].b Numerical model

Another benchmark problem of a rectangle membrane under simple shear that was experimentally investigated by Wong and Pellegrino[13]is taken into account.Figure 16 depicts the corresponding experimental model and the numerical model. As shown in Fig. 16b, the length and width of the Kapton membrane are 0.380 m and 0.128 m, respectively,and its thickness is 2.5 × 10-5m. The tensile Young’s modulus is E+= 3.5×109Pa and the Poisson’s ratio is ν+=1/3,while the compressive Young’s modulus E-and Poisson’s ratio ν-of the membrane are set to be zero. For the boundary and loading conditions,the lower edge of the membrane is fixed while a uniform horizontal displacement Ux=1.6×10-3m is applied on the upper edge and keeps the other degrees constrained.In addition,the upper edge is initially imposed to a uniform pre-stress of 1.5×106Pa in the y direction.To solve this example,the membrane is discretized into 380×128 material points with the same grid size of Δx =Δy =1.0×10-3m.Figure 17a,b shows the x-directional and y-directional displacement contours of the rectangle membrane with the horizon of δ =3.0×10-3m,respectively. Figure 18 depicts the curves of the maximum and minimum principal stresses of the material points along the middle line of the rectangle membrane(i.e.,the A-A section shown in Fig.16b)obtained by the present method with different horizon sizes and the finite element method [27].It can be seen that the maximum and minimum principal stresses obtained by the present method are larger than those obtained by the finite element method due to the difference between the nonlocal and local models.Furthermore,the patterns of the maximum principal stress and the corresponding magnitude/direction obtained by the present method are plotted in Fig. 19a,b, respectively, which are similar to those obtained by Wong and Pellegrino[13]and Zhang et al.[27].In addition,it is noted that all of the external loads are applied by only one load step with 16 iteration steps.

Fig.17 Displacement contours of the membrane obtained by the present method with δ = 3.0×10-3 m. a x-directional displacement. b ydirectional displacement

Fig.18 Curves of the maximum and minimum principal stresses along the middle line of the membrane obtained by the present method with different horizon sizes and the finite element method[27]

6.5 Fracture membrane in tension

Fig.19 Stress patterns.a Maximum principal stress.b Magnitude and direction of the maximum principal stress obtained by the present method

Finally,a fracture problem of a rectangle membrane in tension is taken into account.As shown in Fig.20,the dimension of the membrane is 0.1 m × 0.04 m× 0.001 m together with a 0.006 m vertical notch at the bottom center. For the boundary conditions,a pair of equal and opposite horizontal displacement Uxis applied on the left and right sides. The Young’s moduli are E+= 1.0×105Pa and E-= 0 while the Poisson’s ratio are ν+= 1/3. Then, the membrane is discretized into 200×80 uniform material points.The horizon of the material point is δ =3Δx and the corresponding critical bond stretch is s0=6.0×10-3.The minimum displacement increment is 1.0×10-6m.Figure 21 shows the crack patterns of the fracture membrane at imposed displacements of Uxis 2.1×10-4m,2.4×10-4m,2.7×10-4m and 3.0×10-4m(scale=50).Moreover,it can be seen from the wrinkled level contours that the membrane located close to the crack is almost slack.Different from the definition of the slack,wrinkled and tensional states through the positive and negative states of the principal stresses in the local methods(e.g., the finite element method), it can be simply defined in the present nonlocal model that the membrane where the wrinkled levelandare generally under slack state,wrinkled state and tensional state,respectively.

Fig.20 Fracture problem of a rectangle membrane in tension

Fig.21 Crack and wrinkled level patterns of the fracture membrane at imposed displacements of a Ux =2.1×10-4 m,b Ux =2.4×10-4 m,c Ux =2.7×10-4 m and d Ux =3.0×10-4 m(scale=50)

7 Conclusions

This paper presents an implicit computational framework for the quasi-static mechanical analysis of the bimodular structures and membranes with large deformation based on the bond-based peridynamic theory. In this framework, the constant and tangential stiffness matrices of the implicit peridynamic formulations for the static equilibrium problems are presented.The former is derived from the geometric approximation by considering the bond strain as the linear function of the displacements of the material points approximately and decomposing the current displacements of the material points as the sum of the initial displacements and the incremental displacements. The latter is derived directly by the linearization of the pairwise force using first-order Taylor’s expansion.In addition,the bimodular material model in peridynamics is proposed, in which the tensile or compressive states of the material at each material point is dexterously described by the corresponding tensile or compressive states of the bonds in its neighborhood. Moreover, the bimodular material model in peridynamics is easily extended to simulate the wrinkling and fracture analyses of the membranes by setting the compressive modulus to be zero.In addition,the incremental-iterative algorithmwith displacement or load control strategy is adopted to obtain the convergent solutions and the computational framework is listed for the nonlinear problems.At last,the results of the examples show the accuracy and efficiency of the present method for the quasi-static large deformation, wrinkling and fracture analyses of the bimodular structures and membranes by comparing it with the reference solutions.The results also show that the convergent behaviors of the present method are very good when the tangential stiffness matrix is utilized for these nonlinear problems.However,the convergent errors of the present method could be oscillated periodically when the constant stiffness matrix is applied for these nonlinear problems.For the computational efficiency,in spite of that the computational time for each iteration step in the present method is longer than that in the explicit peridynamic method,the total number of the iteration step in the present method is much less than that in the explicit peridynamic method. Consequently, the computational cost of the present method is less than that of the latter in general.Finally,the present method in this paper provides an effective numerical simulating route to study the mechanical characteristics of the bimodular materials and structures in aerospace engineering.

AcknowledgementsThe work was supported by the National Natural Science Foundation of China (Grants 11672062, 11772082, and 11672063), the 111 Project (Grant B08014), and the Fundamental Research Funds for the Central Universities.