CN102930079B - Method for analyzing interlaminar damage of composite material laminate - Google Patents

Abstract

The invention discloses a method for analyzing interlaminar damage of a composite material laminate. The degree of freedom of an auxiliary node of a combined interface unit is condensed out by deducing a rigidity equation of the combined interface unit, so that the combined interface unit is a whole in a finite element model and only has eight main nodes. By the combined interface unit and the finite element model, a simulation result is relatively good; a curve obtained by calculation is basically overlapped with a test curve on a linear segment and is consistent with the tendency of the test curve on a falling segment behind a peak value; simulation of the start of the interlaminar damage and simulation of a damage evolution process accord with the reality; and a load peak value obtained by calculation is 62.8 N, a test value is 63.1 N, and the error is -0.4 percent. Compared with the prior art, the method has the characteristics of small calculation scale and high calculation efficiency; and furthermore, the combined interface unit has a thickness and is easy to implement compared with a cohesion unit with zero thickness during pretreatment modeling.

Description

A kind of method of COMPOSITE MATERIALS laminate interlayer damage

Technical field

The invention belongs to Computational Mechanics field, be specifically related to a kind of for calculating the method for analog composite material layer plywood interlayer damage.

Background technology

Composite laminated plate is widely used in a lot of fields, especially very common in aeronautic structure.Laminate is bonded together and is formed whole structural slab by multilayer lamina, and each lamina is called a sublayer.The material principal direction of lamina can be pressed as required to different directions and different order and lay, obtain the laminate of different mechanical properties, so the mechanical property of laminate have designability.Interlayer layering is one of main damage type of laminate, and the relative in-plane strength of interlaminar strength makes a little less than laminate easily produce interlayer damage, especially more obvious under percussive action.The appearance of interlayer damage reduces structural stability and residual intensity greatly, and may cause catastrophic consequence.The how interlayer of Accurate Prediction and COMPOSITE MATERIALS laminate damage, is used composite structure very important to appropriate design.

Cohesion elements method (Cohesive Zone Model s, CZM) is that present analysis laminate interlayer damages highly effective Finite Element Method.CZM method is carried out the interlayer glue-line of simulation layer plywood by introducing cohesion boundary element, between each layer of laminate, or in structure, may there is the location arrangements cohesion boundary element (Cohesive Element) of layering, when cohesion boundary element when the stress value under other unit effects reaches boundary strength around interlayer damage initial; When the further bearing load of cohesion boundary element, and energy release rate interlayer damage expansion while reaching critical value.CZM method is easy to finite element program and realizes, and does not need to provide the initial sum of damaging between just can prediction interval initial damage position and expands.Yet one of shortcoming of CZM method is in order to obtain rationally result accurately, just need finite element model to divide very meticulous grid, can make result of calculation enough reasonable.At present, during most researchers analysis layer plywood interlayer performance, all use solid element and cohesion boundary element to carry out modeling: with solid element simulation, between the solid element of neighbouring sublayer, to insert the interlayer glue-line that cohesion boundary element that zero thickness or thickness are very little carrys out simulation layer plywood in the face of sublayer.It is just large that solid element is compared the calculated amount of plate shell unit own, makes to calculate to be on a grand scale when grid division is very thin, and this has just limited the future in engineering applications of CZM method.

Summary of the invention

For overcoming the be on a grand scale deficiency of restriction future in engineering applications of the calculating existing in prior art, the present invention proposes the method for a kind of COMPOSITE MATERIALS laminate interlayer damage.

Detailed process of the present invention is:

Step 1, sets up finite element model;

Step 2, extracts displacement of joint and the coordinate of each combined interface unit; The process of extracting comprises:

A. extract displacement of joint; From the finite element model having established, extract the displacement of each combined interface unit major node; This major node displacement is expressed as by vector form:

$u_M{(u_M^{1-},v_M^{1-},w_M^{1-},{\mathrm{\θ}}_{\mathrm{Mx}}^{1-},{\mathrm{\θ}}_{\mathrm{My}}^{1-}...,u_M^{4+},v_M^{4+},w_M^{4+},{\mathrm{\θ}}_{\mathrm{Mx}}^{4+},{\mathrm{\θ}}_{\mathrm{My}}^{4+})}_{40\×1}^T---\left(1\right)$

In formula (1), u _mthe major node displacement of combined interface unit, represent that respectively i major node is to along x, y, the translation displacement of z direction (i=1,2,3,4); represent that respectively i major node is to around x, the rotation displacement of y axle, subscript negative sign represents the displacement component of lower major node, and plus sige represents the displacement component of major node, and node numbering is as shown in Figure 2; Described major node is to referring to x in same combined interface unit, and y coordinate is identical, two major nodes that z coordinate is different;

B. extract node coordinate; From the finite element model having established, extract the coordinate of each combined interface unit major node; This major node coordinate is expressed as by vector form:

$x_M={(x_M^{1-},y_M^{1-},z_M^{1-},...,x_M^{4+},y_M^{4+},z_M^{4+})}_{24\×1}^T$

In formula, x _mthe major node coordinate of combined interface unit, respectively i the x that major node is right, y, z coordinate (i=1,2,3,4), subscript negative sign represents lower major node coordinate components, plus sige represents major node coordinate components;

The node of described combined interface unit refers to the major node of combined interface unit; Described major node refers to the node that combined interface unit is connected with plate shell unit;

Step 3, calculates the geometric matrix of each combined interface unit; This process comprises the following steps:

A. calculate the right coordinate of the secondary node in each combined interface unit; According to extracting the coordinate of each combined interface unit major node obtaining in formula (2), by formula (3):

$x_S=\frac{1}{2}\left[I_{12\×12}\right|I_{12\×12}]x_M---\left(3\right)$

Obtain the right coordinate array x of the secondary node in combined interface unit _s:

$x_S={(x_S^1,y_S^1,z_S^1,...,x_S^4,y_S^4,z_S^4)}_{12\×1}^T---\left(4\right)$

In formula (3), x _sthe right coordinate array of secondary node of combined interface unit, I _{12 * 12}it is unit diagonal matrix; In formula (4), respectively i the right x of secondary node, y, z coordinate (i=1,2,3,4);

Described secondary node is in the secondary node that refers to inside, combined interface unit and comprise, a pair of secondary node that position overlaps mutually; Described secondary node refers to the node that rigidity unit is connected with cohesion unit;

B. calculate the shape function matrix of each combined interface unit; By formula (5), calculate shape function matrix

$H(\mathrm{\ξ},\mathrm{\η})=\left[\begin{array}{ccccccc}{N}_1& 0& 0& & N_4& 0& 0\\ 1& N_1& 0& |...|& 0& N_4& 0\\ 0& 0& N_1& & 0& 0& N_4\end{array}\right]---\left(5\right)$

In formula (5), H (ξ, η) is shape function matrix; N _ibe shape function (i=1,2,3,4), be specially:

N ₁=(1-ξ)(1-η)/4，N ₂=(1+ξ)(1-η)/4（6）

N ₃=(1+ξ)(1+η)/4，N ₄=(1-ξ)(1+η)/4

In formula (6), ξ and η are the coordinates under natural system of coordinates, ξ, η ∈ [1,1];

C. calculate the transition matrix of each combined interface unit; By formula (7) compute vector v _ξand v _ηx, y, z component:

${v_{\mathrm{\ξ}}}_i=\underset{k=1}{\overset{4}{\mathrm{\Σ}}}\frac{\∂N_k}{\∂\mathrm{\ξ}}\·x_{\mathrm{Si}}^k$ (i＝1，2，3)(7)

$v_{{\mathrm{\η}}_i}=\underset{k=1}{\overset{4}{\mathrm{\Σ}}}\frac{\∂N_k}{\∂\mathrm{\η}}\·x_{\mathrm{Si}}^k$

By the vector v obtaining _ξand v _η, by formula (8) compute vector v _s, v _t, v _n

v _n=(v _ξ×v _η)||v _ξ×v _η|| ^-1

v _s=v _ξ||v _ξ|| ^-1（8）

v _t＝v _n×v _s

By the vector v obtaining _s, v _t, v _n, by formula (9), calculate transition matrix

Θ=(v _s，v _t,v _n)（9）

In formula (9), Θ is transition matrix; In formula (7), it is vector v _ξi component, i=1,2,3 corresponding x respectively, y, z coordinate; it is vector v _ηi component, i=1,2,3 corresponding x respectively, y, z coordinate; be k secondary node to the coordinate in i direction, by formula (3), obtain;

D. calculate the geometric matrix of each combined interface unit; By formula (10) computational geometry matrix

B=Θ ^TH(ξ，η)Φ（10）

In formula (10), Θ ^tit is the transposed matrix of transition matrix Θ; H (ξ, η) is the shape function matrix in formula (5); Φ=(I _{12 * 12}| I _{12 * 12}), I _{12 * 12}it is unit diagonal matrix;

Step 4, calculates the elastic matrix of each combined interface unit; This process comprises the following steps:

A. calculate the secondary displacement of joint of each combined interface unit; According to formula (1), extract the displacement of each combined interface unit major node obtaining, by formula (11)

u _S=T _δu _M（11）

Obtain the secondary displacement of joint of each combined interface unit:

$u_S={(u_S^{1-},v_S^{1-},w_S^{1-},{\mathrm{\θ}}_{\mathrm{Sx}}^{1-},{\mathrm{\θ}}_{\mathrm{Sy}}^{1-}...,u_S^{4+},v_S^{4+},w_S^{4+},{\mathrm{\θ}}_{\mathrm{Sx}}^{4+},{\mathrm{\θ}}_{\mathrm{Sy}}^{4+})}_{40\×1}^T---\left(12\right)$

In formula (12), u _sthe displacement of the secondary node of combined interface unit, represent that respectively i secondary node is to along x, y, the translation displacement of z direction (i=1,2,3,4); represent that respectively i secondary node is to around x, the rotation displacement of y axle; Subscript negative sign represents the displacement component of lower secondary node, and plus sige represents the displacement component of secondary node; T _δbe displacement transition matrix between the major node of combined interface unit and secondary node, be specially

$T_{\mathrm{\δ}}={\left[\begin{array}{cccc}{t}_{{\mathrm{\δ}1}^{-}}& & & 0\\ & t_{{\mathrm{\δ}2}^{-}}& & \\ & & ...& \\ 0& & & t_{{\mathrm{\δ}4}^{+}}\end{array}\right]}_{40\×40}---\left(13\right)$

In formula (13) be displacement transition matrix between the right major node of i rigidity unit and secondary node (i=1 ..., 4), the subscript of i is that negative sign represents lower rigidity unit, plus sige represents rigidity unit; concrete form be

$t_{{\mathrm{\δi}}^{\±}}=\left[\begin{array}{ccccc}1& 0& 0& 0& z_S^{i\±}-z_M^{i\±}\\ 0& 1& 0& -(z_S^{i\±}-z_M^{i\±})& 0\\ 0& 0& 1& y_S^{i\±}-y_M^{i\±}& -(x_S^{i\±}-x_M^{i\±})\\ 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 1\end{array}\right]---\left(14\right)$

B. calculate local relative displacement; Described local relative displacement refers to the relative displacement under local coordinate system; By formula (15), calculate local relative displacement

In formula (15), δ is local relative displacement, and B is the geometric matrix in formula (10), it is the transition matrix of using while removing rotational freedom; The described transition matrix of using while removing rotational freedom is specially

${\mathrm{\γ}}_{3\×5}=\left[\begin{array}{ccccc}1& 0& 0& 0& 0\\ 0& 1& 0& 0& 0\\ 0& 0& 1& 0& 0\end{array}\right]---\left(16\right)$

C. calculate total relative displacement of mixed type interlayer damage; Described mixed type interlayer damage refers to and comprises I type simultaneously, the interlayer type of impairment of II type and the damage of III type; Described total relative displacement refers to mixed type relative storey displacement amount; By formula (17), calculate total relative displacement of mixed type interlayer damage

${\mathrm{\&delta;}}_m=\sqrt{{\mathrm{\&delta;}}_1^2+{\mathrm{\&delta;}}_2^2+{<{\mathrm{\&delta;}}_3>}^2}=\sqrt{{\mathrm{\&delta;}}_{\mathrm{shear}}^2+{<{\mathrm{\&delta;}}_3>}^2}---\left(17\right)$

In formula (17), δ _mit is total relative displacement of mixed type interlayer damage; δ ₁, δ ₂, δ ₃be local relative displacement δ respectively along the x under local coordinate system, y, three components of z direction;

D. displacement calculating mixing ratio; Described displacement mixing ratio refers to the ratio of interlayer shear displacement and interlayer opening displacement, specifically

$\mathrm{\&beta;}=\frac{{\mathrm{\&delta;}}_{\mathrm{shear}}}{{\mathrm{\&delta;}}_3}---\left(18\right)$

In formula (18), β is displacement mixing ratio;

E. calculate mixed type interlayer and damage initial displacement and ultimate failure displacement; Described mixed type interlayer damages initial displacement and refers to that the damage of mixed type interlayer is initial when the total relative displacement of mixed type arrives this value; The displacement of described mixed type interlayer ultimate failure refers to that the damage of mixed type interlayer is developed when the total relative displacement of mixed type arrives this value; Specifically

${\mathrm{\&delta;}}_m^0=\left\{\begin{array}{cc}{\mathrm{\&delta;}}_3^0{\mathrm{\&delta;}}_1^0\sqrt{\frac{1+{\mathrm{\&beta;}}^2}{{\left({\mathrm{\&delta;}}_1^0\right)}^2+{\left(\mathrm{\&beta;}{\mathrm{\&delta;}}_3^0\right)}^2}},& {\mathrm{\&delta;}}_3>0\\ {\mathrm{\&delta;}}_{\mathrm{shear}}^0& {\mathrm{\&delta;}}_3\&le;0\end{array}\right.---\left(19\right)$

${\mathrm{\&delta;}}_m^f=\left\{\begin{array}{cc}\frac{2}{K{\mathrm{\&delta;}}_m^0}[G_{\mathrm{IC}}+(G_{\mathrm{IIC}}-G_{\mathrm{IC}}){\left(\frac{{\mathrm{\&beta;}}^2}{1+{\mathrm{\&beta;}}^2}\right)}^{\mathrm{\&eta;}}],& {\mathrm{\&delta;}}_3>0\\ {\mathrm{\&delta;}}_1^f,& {\mathrm{\&delta;}}_3\&le;0\end{array}\right.---\left(20\right)$

In formula (19), that mixed type interlayer damages initial displacement; single I, II, III type interlayer damages the initial displacement of corresponding damage; N, S, T is respectively corresponding boundary strength, and S=T; k is interface rigidity, is taken as K=10 ⁶; Described single I, II, the damage of III type interlayer refers to and only comprises I, II, the damage of III type is a kind of type of impairment wherein; In formula (20), it is the displacement of mixed type interlayer ultimate failure; G _iC, G _iIC, G _iIICrespectively I, II, the critical value of III type interlayer damage energy release rate; single I, II, II I type interlayer damages corresponding ultimate failure displacement, ${\mathrm{\&delta;}}_1^f=2G_{\mathrm{IIC}}/S,$ ${\mathrm{\&delta;}}_2^f=2G_{\mathrm{IIIC}}/T;$ B-K parameter η is by testing acquisition;

F. calculate damage variable; Described damage variable refers to for characterizing the amount of mixed type interlayer faulted condition; Described faulted condition refers to the state of injury severity score; By formula (21), calculate mixed type interlayer damage variable

$d=\frac{{\mathrm{\&delta;}}_m^f({\mathrm{\&delta;}}_m^{\max }-{\mathrm{\&delta;}}_m^0)}{{\mathrm{\&delta;}}_m^{\max }({\mathrm{\&delta;}}_m^f-{\mathrm{\&delta;}}_m^0)},0\&le;d\&le;1---\left(21\right)$

In formula (21), d is damage variable, δ _mthe maximal value reaching,

G. calculating elastic matrix; Elastic matrix is divided into 5 kinds of forms according to the difference of loading, is specially:

When time, unmarred elastic stage:

When and δ ₃>=0 o'clock, the damage softening stage: $D=\left[\begin{array}{ccc}(1-d)K& & \\ & (1-d)K& \\ & & (1-d)K\end{array}\right],$

When and δ ₃during <0, the damage softening stage: $D=\left[\begin{array}{ccc}(1-d)K& & \\ & (1-d)K& \\ & & (1-d)K\end{array}\right],$

When and δ ₃>=0 o'clock, damage the stage completely:

When and δ ₃during <0, damage the stage completely:

Step 5, calculates each combined interface unit joint forces; Described combined interface unit joint forces refers to the joint forces of major node, by formula (22), calculates:

In formula (22), F is the joint forces of each combined interface unit major node, and τ is stress, by formula (23), calculates

τ＝Dδ（23）

In formula (23), τ is stress, and D is elastic matrix, and δ is the local relative displacement in formula (15);

Step 6, calculates each combined interface element stiffness matrix; By formula (24), calculate the element stiffness matrix of each combined interface unit

In formula (24), K is the element stiffness matrix of each combined interface unit; T _δdisplacement transition matrix between the middle rigidity unit's major node of formula (11) and secondary node, it is the transposed matrix of T δ; the transition matrix of using while removing rotational freedom in formula (16), be transposed matrix; B is the geometric matrix in formula (10), B ^tit is the transposed matrix of B; D is elastic matrix;

Step 7, assembles the global stiffness matrix of whole finite element model.Described assembling global stiffness matrix is before separating balance equation, by finite element software, the element stiffness matrix of whole finite element model is assembled into global stiffness matrix.

Step 8, solves balance equation; Adopt Numerical Iteration Method, by finite element software, solve described balance equation; In solving, after completing an iteration, whether the displacement that judges all nodes in finite element model restrains, if obtain the convergence solution of all nodal displacements in finite element model, this iteration step is complete, otherwise repeating step 2～step 8, until obtain the convergence solution of all nodal displacements in finite element model.

The object of the present invention is to provide a kind of for calculating the method for analog composite material layer plywood interlayer damage.

The invention provides a kind of applicable combined interface unit being used in conjunction with plate shell unit, can consider translation displacement and the impact of rotation displacement on interlayer damage of plate shell unit node simultaneously.By the stiffness equations of derivation combined interface unit, the degree of freedom of the secondary node in combined interface unit is fallen in polycondensation, and make combined interface unit is an integral body in finite element model, only has eight major nodes.Compare with traditional modeling method of being combined with cohesion boundary element in current entity unit, it is less that finite element model of the present invention calculates scale, can improve counting yield.In model, use plate shell unit, make computational accuracy under same mesh density higher than Model of Solid Elements.In addition, because combined interface unit has thickness, so the cohesion unit of comparing zero thickness during pre-treatment modeling is very easily realized.

Shown in Fig. 4, be the double cantilever beam computation model that embodiments of the invention are described, this double cantilever beam computation model respectively in two beam arms face position set up four Node Quadrilateral Finite Element plate shell unit S4R; Except in the region of prefabricated initial crack, with the combined interface unit of limited thickness, connect upper and lower plate shell unit, record the load displacement course of loading end.In the DCB of embodiment computation model, as shown in Figure 5, in Fig. 5, the represented curve of solid line is test findings for the loading end load-displacement curves calculating and the comparing result of trial value, and the represented curve of dotted line is result of calculation; Y axle represents load value, and x axle represents load deflection.Adopt combined interface of the present invention unit and finite element model to calculate analog result better, calculating curve overlaps with trial curve substantially at linearity range, descending branch after peak value is consistent with the trend of trial curve, and the simulation of interlayer damage initial sum damage evolutionary process is all realistic.The load peaks 62.8N calculating, trial value 63.1N, error-0.4%.

Below in conjunction with drawings and Examples, the present invention is further described.

Accompanying drawing explanation

Fig. 1 is laminate interlayer breakdown diagnosis method flow diagram

Fig. 2 is finite element model schematic diagram of the present invention;

Fig. 3 is finite element model and combined interface cell cross-section schematic diagram;

Fig. 4 is double cantilever beam computation model;

Fig. 5 is the load-displacement curves result contrast of the loading end of double cantilever beam experiment and computation simulation;

Fig. 6 A is double cantilever beam finite element model grid, and this size of mesh opening is 4mm * 4mm;

Fig. 6 B is grid chart after the distortion of double cantilever beam finite element model, and this size of mesh opening is 4mm * 4mm.In figure:

1. upper plate shell unit 2. lower plate shell unit 3. 5.Xia sublayers, 4.Shang sublayer, combined interface unit

6. go up on major node 7. the cohesion boundary element of secondary node 9. zero thickness in rigidity unit 8.

10. 12. times major node 13. double cantilever beam test findings of 11. times rigidity of time secondary node unit

14. double cantilever beams calculate the plate shell unit of 16. times beam arms of plate shell unit of beam arm in analog result 15.

17. are used for connecting the combined interface unit of beam arm and lower beam arm

Embodiment

The present embodiment is the method for COMPOSITE MATERIALS laminate interlayer damage in a kind of double cantilever beam test.

As shown in Figure 2, by combined interface unit 3 simulation layer plywood interlayer glue-lines, in the face of adjacent sublayers, with upper plate shell unit 1 and lower plate shell unit 2, simulate respectively, the finite element model that combined interface unit 3 and upper plate shell unit 1 and lower plate shell unit 2 are formed is used for the interlayer damage process of analog computation laminate, reduce the scale of traditional calculations model, reduce calculated amount.

Fig. 3 is the cross sectional representation of finite element model shown in Fig. 2.As shown in Figure 3,4 use upper plate shell unit 1 simulations of upper sublayer, 5 use lower plate shell unit 2 simulations of lower sublayer, combined interface unit 3 simulations for interlayer glue-line.Combined interface unit 3 is that zero cohesion boundary element combines by eight rigidity units and a thickness.

A part of take in combined interface unit illustrates its structure as example: the node 6 that upper rigidity unit 7 is connected with upper plate shell unit 1 is called upper major node, the node 8 that upper rigidity unit 7 is connected with the cohesion boundary element 9 of zero thickness is called upper secondary node, the external applied load that upper major node 6 is subject to is delivered on the cohesion boundary element 9 of zero thickness by indeformable upper rigidity unit 7, and the process of interlayer damage is that the continuous decay by cohesion boundary element 9 rigidity of zero thickness realizes.Similarly, the external applied load that lower major node 12 is subject to is delivered on the cohesion boundary element 9 of zero thickness by indeformable lower rigidity unit 11.Upper major node 6 is called a major node pair with lower major node 12; It is right that upper rigidity unit 7 and lower rigidity unit 11 are called a rigidity unit; Upper secondary node 8 and lower secondary node 10 are called a secondary node pair.The cohesion boundary element of all rigidity units and zero thickness forms a whole, namely combined interface unit.As shown in Figure 2, combined interface unit is the boundary element with eight major nodes that has thickness, and each node comprises five degree of freedom, i.e. three translational degree of freedom u, v, w and two rotational freedom θ _x, θ _y.

As shown in Figure 4, DCB test specimen is the unidirectional fibre enhancement layer plywood of material T300/977-2, and machine direction is along plate length direction.The long L=150mm of DCB test specimen, wide w=20mm, the thick h=1.98mm of single beam arm, middle position near the prefabricated Initial crack length Lc=55mm of loading end.The clamped constraint of finite element model right-hand member, left end applies to upper beam arm the load that direction makes progress, and lower beam arm applies the downward load of direction.The material properties of lamina is as shown in table 1.

Table 1 T300/977-2 material properties

The detailed process of the present embodiment comprises the following steps:

Step 1, sets up finite element model.The described finite element model of setting up is the pre-treatment modeling function by finite element software, in neutral surface position separately, each sublayer of composite laminated plate, sets up geometric surface; Described sublayer refers to a lamina of laminate; Each sublayer of described laminate neutral surface separately refers to the central plane of each lamina thickness direction.In the present embodiment, the finite element software adopting is ABAQUS.The process of setting up finite element model comprises:

A. grid division and defined attribute.By finite element software, adopt conventional method that the neutral surface of each sublayer of composite laminated plate is divided into quadrilateral mesh.There is each neutral surface of grid to be defined as plate shell unit division, and give described plate shell unit material properties and material principal direction;

B. between the plate shell unit of neighbouring sublayer, generate the combined interface unit that has thickness, the finite element model that obtains establishing.Described combined interface unit is to be formed with cohesion unit combination by rigidity unit.

Step 2, extracts displacement of joint and the coordinate of each combined interface unit.The node of described combined interface unit refers to the major node of combined interface unit; Described major node refers to the node that combined interface unit is connected with plate shell unit.The process of extracting comprises:

A. extract displacement of joint.From the finite element model having established, extract the displacement of each combined interface unit major node.This major node displacement is expressed as by vector form:

$u_M={(u_M^{1-},v_M^{1-},w_M^{1-},{\mathrm{\&theta;}}_{\mathrm{Mx}}^{1-},{\mathrm{\&theta;}}_{\mathrm{My}}^{1-}...,u_M^{4+},v_M^{4+},w_M^{4+},{\mathrm{\&theta;}}_{\mathrm{Mx}}^{4+},{\mathrm{\&theta;}}_{\mathrm{My}}^{4+})}_{40\&times;1}^T---\left(25\right)$

In formula (25), u _mthe major node displacement of combined interface unit, represent that respectively i major node is to along x, y, the translation displacement of z direction (i=1,2,3,4); represent that respectively i major node is to around x, the rotation displacement of y axle, subscript negative sign represents the displacement component of lower major node, and plus sige represents the displacement component of major node, and node numbering is as shown in Figure 2.Described major node is to referring to x in same combined interface unit, and y coordinate is identical, two major nodes that z coordinate is different.

B. extract node coordinate.From the finite element model having established, extract the coordinate of each combined interface unit major node.This major node coordinate is expressed as by vector form:

$x_M={(x_M^{1-},y_M^{1-},z_M^{1-},...,x_M^{4+},{\mathrm{y\&theta;}}_M^{4+},{\mathrm{z\&theta;}}_M^{4+})}_{24\&times;1}^T---\left(26\right)$

In formula, x _mthe major node coordinate of combined interface unit, respectively i the x that major node is right, y, z coordinate (i=1,2,3,4), subscript negative sign represents lower major node coordinate components, plus sige represents major node coordinate components.

Step 3, calculates the geometric matrix of each combined interface unit.This process comprises the following steps:

A. calculate the right coordinate of the secondary node in each combined interface unit.Described secondary node is in the secondary node that refers to inside, combined interface unit and comprise, a pair of secondary node that position overlaps mutually; Described secondary node refers to the node that rigidity unit is connected with cohesion unit.According to extracting the coordinate of each combined interface unit major node obtaining in formula (26), by formula (27):

$x_S=\frac{1}{2}\left[I_{12\&times;12}\right|I_{12\&times;12}]x_M---\left(27\right)$

Obtain the right coordinate array x of the secondary node in combined interface unit _s:

$x_S={(x_S^1,y_S^1,z_S^1,...,x_S^4,y_S^4,z_S^4)}_{12\&times;1}^T---\left(28\right)$

In formula (27), x _sthe right coordinate array of secondary node of combined interface unit, I _{12 * 12}it is unit diagonal matrix.In formula (28), respectively i the right x of secondary node, y, z coordinate (i=1,2,3,4).

B. calculate the shape function matrix of each combined interface unit.By formula (29), calculate shape function matrix

$H(\mathrm{\&xi;},\mathrm{\&eta;})=\left[\begin{array}{ccccccc}{N}_1& 0& 0& & N_4& 0& 0\\ 0& N_1& 0& |...|& 0& N_4& 0\\ 0& 0& N_1& & 0& 0& N_4\end{array}\right]---\left(29\right)$

In formula (29), H (ξ, η) is shape function matrix; N _ibe shape function (i=1,2,3,4), be specially:

N ₁=(1-ξ)(1-η)/4，N ₂=(1+ξ)(1-η)/4（30）

N ₃=(1+ξ)(1+η)/4，N ₄=(1-ξ)(1+η)/4

In formula (30), ξ and η are the coordinates under natural system of coordinates, ξ, η ∈ [1,1].

C. calculate the transition matrix of each combined interface unit.By formula (31) compute vector v _ξand v _ηx, y, z component:

${v_{\mathrm{\&xi;}}}_i=\underset{k=1}{\overset{4}{\mathrm{\&Sigma;}}}\frac{\&PartialD;N_k}{\&PartialD;\mathrm{\&xi;}}\&CenterDot;x_{\mathrm{Si}}^k$ (i＝1，2，3)(31)

${v_{\mathrm{\&eta;}}}_i=\underset{k=1}{\overset{4}{\mathrm{\&Sigma;}}}\frac{\&PartialD;N_k}{\&PartialD;\mathrm{\&eta;}}\&CenterDot;x_{\mathrm{Si}}^k$

By the vector v obtaining _ξand v _η, by formula (32) compute vector v _s, v _t, v _n

v _n=(v _ξ×v _η)||v _ξ×v _η|| ^-1

v _s=v _ξ||v _ξ|| ^-1（32）

v _t＝v _n×v _s

By the vector v obtaining _s, v _t, v _n, by formula (33), calculate transition matrix

Θ=(v _s，v _t，v _n)（33）

In formula (33), Θ is transition matrix.In formula (31), it is vector v _ξi component, i=1,2,3 corresponding x respectively, y, z coordinate; it is vector v _ηi component, i=1,2,3 corresponding x respectively, y, z coordinate. be k secondary node to the coordinate in i direction, by formula (27), obtain.

D. calculate the geometric matrix of each combined interface unit.By formula (34) computational geometry matrix

B=Θ ^TH(ξ，η)Φ（34）

In formula (34), Θ ^tit is the transposed matrix of transition matrix Θ; H (ξ, η) is the shape function matrix in formula (29); Φ=(I _{12 * 12}| I _{12 * 12}), I _{12 * 12}it is unit diagonal matrix.

Step 4, calculates the elastic matrix of each combined interface unit.This process comprises the following steps:

A. calculate the secondary displacement of joint of each combined interface unit.According to formula (25), extract the displacement of each combined interface unit major node obtaining, by formula (35)

u _S=T _δu _M（35）

Obtain the secondary displacement of joint of each combined interface unit:

$u_S={(u_S^{1-},v_S^{1-},w_S^{1-},{\mathrm{\&theta;}}_{\mathrm{Sx}}^{1-},{\mathrm{\&theta;}}_{\mathrm{Sy}}^{1-}...,u_S^{4+},v_S^{4+},w_S^{4+},{\mathrm{\&theta;}}_{\mathrm{Sx}}^{4+},{\mathrm{\&theta;}}_{\mathrm{Sy}}^{4+})}_{40\&times;1}^T---\left(36\right)$

In formula (36), u _sthe displacement of the secondary node of combined interface unit, represent that respectively i secondary node is to along x, y, the translation displacement of z direction (i=1,2,3,4); represent that respectively i secondary node is to around x, the rotation displacement of y axle; Subscript negative sign represents the displacement component of lower secondary node, and plus sige represents the displacement component of secondary node.T _δbe displacement transition matrix between the major node of combined interface unit and secondary node, be specially

$T_{\mathrm{\&delta;}}={\left[\begin{array}{cccc}{t}_{{\mathrm{\&delta;}}^{-}}& & & 0\\ & t_{{\mathrm{\&delta;}2}^{-}}& & \\ & & ...& \\ 0& & & t_{{\mathrm{\&delta;}}^{+}}\end{array}\right]}_{40\&times;40}---\left(37\right)$

In formula (37) be displacement transition matrix between the right major node of i rigidity unit and secondary node (i=1 ..., 4), the subscript of i is that negative sign represents lower rigidity unit, plus sige represents rigidity unit. concrete form be

$t_{{\mathrm{\&delta;i}}^{\&PlusMinus;}}=\left[\begin{array}{ccccc}1& 0& 0& 0& z_S^{i\&PlusMinus;}-z_M^{i\&PlusMinus;}\\ 0& 1& 0& -(z_S^{i\&PlusMinus;}-z_M^{i\&PlusMinus;})& 0\\ 0& 0& 1& y_S^{i\&PlusMinus;}-y_M^{i\&PlusMinus;}& -(x_S^{i\&PlusMinus;}-x_M^{i\&PlusMinus;})\\ 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 1\end{array}\right]---\left(38\right)$

B. calculate local relative displacement.Described local relative displacement refers to the relative displacement under local coordinate system.By formula (39), calculate local relative displacement

In formula (39), δ is local relative displacement, and B is the geometric matrix in formula (34), it is the transition matrix of using while removing rotational freedom.The described transition matrix of using while removing rotational freedom is specially

${\mathrm{\&gamma;}}_{3\&times;5}=\left[\begin{array}{ccccc}1& 0& 0& 0& 0\\ 0& 1& 0& 0& 0\\ 0& 0& 1& 0& 0\end{array}\right]---\left(40\right)$

C. calculate total relative displacement of mixed type interlayer damage.Described mixed type interlayer damage refers to and comprises I type simultaneously, the interlayer type of impairment of II type and the damage of III type; Described total relative displacement refers to mixed type relative storey displacement amount.By formula (41), calculate total relative displacement of mixed type interlayer damage

${\mathrm{\&delta;}}_m=\sqrt{{\mathrm{\&delta;}}_1^2+{\mathrm{\&delta;}}_2^2+{<{\mathrm{\&delta;}}_3>}^2}=\sqrt{{\mathrm{\&delta;}}_{\mathrm{shear}}^2+{<{\mathrm{\&delta;}}_3>}^2}$

In formula (41), δ _mit is total relative displacement of mixed type interlayer damage; δ ₁, δ ₂, δ ₃be local relative displacement δ respectively along the x under local coordinate system, y, three components of z direction;

D. displacement calculating mixing ratio.Described displacement mixing ratio refers to the ratio of interlayer shear displacement and interlayer opening displacement, specifically

$\mathrm{\&beta;}=\frac{{\mathrm{\&delta;}}_{\mathrm{shear}}}{{\mathrm{\&delta;}}_3}---\left(42\right)$

In formula (42), β is displacement mixing ratio.

E. calculate mixed type interlayer and damage initial displacement and ultimate failure displacement.Described mixed type interlayer damages initial displacement and refers to that the damage of mixed type interlayer is initial when the total relative displacement of mixed type arrives this value; The displacement of described mixed type interlayer ultimate failure refers to that the damage of mixed type interlayer is developed when the total relative displacement of mixed type arrives this value.Specifically

${\mathrm{\&delta;}}_m^0=\left\{\begin{array}{cc}{\mathrm{\&delta;}}_3^0{\mathrm{\&delta;}}_1^0\sqrt{\frac{1+{\mathrm{\&beta;}}^2}{{\left({\mathrm{\&delta;}}_1^0\right)}^2+{\left(\mathrm{\&beta;}{\mathrm{\&delta;}}_3^0\right)}^2}},& {\mathrm{\&delta;}}_3>0\\ {\mathrm{\&delta;}}_{\mathrm{shear}}^0,& {\mathrm{\&delta;}}_3\&le;0\end{array}\right.---\left(43\right)$

${\mathrm{\&delta;}}_m^f=\left\{\begin{array}{cc}\frac{2}{K{\mathrm{\&delta;}}_m^0}[G_{\mathrm{IC}}+(G_{\mathrm{IIC}}-G_{\mathrm{IC}}){\left(\frac{{\mathrm{\&beta;}}^2}{1+{\mathrm{\&beta;}}^2}\right)}^{\mathrm{\&eta;}}]& ,{\mathrm{\&delta;}}_3>0\\ {\mathrm{\&delta;}}_1^f,& {\mathrm{\&delta;}}_3\&le;0\end{array}\right.---\left(44\right)$

In formula (43), that mixed type interlayer damages initial displacement; single I, II, III type interlayer damages the initial displacement of corresponding damage; N, S, T is respectively corresponding boundary strength, and S=T; k is interface rigidity, is taken as K=10 ⁶.Described single I, II, the damage of III type interlayer refers to and only comprises I, II, the damage of III type is a kind of type of impairment wherein.In formula (44), it is the displacement of mixed type interlayer ultimate failure; G _iC, G _iIC, G _iIICrespectively I, II, the critical value of III type interlayer damage energy release rate; single I, II, II I type interlayer damages corresponding ultimate failure displacement, b-K parameter η is by testing acquisition.Described B-K parameter η is according to MLBenzeggagh, M Kenane, at " Measurement of mix-mode delamination fracture toughness of unidirectional glass/epoxy composites with mixed-mode bending apparatus ", the test method of middle proposition obtains, the document is disclosed in Composite Science and Technology, 56:439-49.

F. calculate damage variable.Described damage variable refers to for characterizing the amount of mixed type interlayer faulted condition; Described faulted condition refers to the state of injury severity score.By formula (45), calculate mixed type interlayer damage variable

$d=\frac{{\mathrm{\&delta;}}_m^f({\mathrm{\&delta;}}_m^{\max }-{\mathrm{\&delta;}}_m^0)}{{\mathrm{\&delta;}}_m^{\max }({\mathrm{\&delta;}}_m^f-{\mathrm{\&delta;}}_m^0)},0\&le;d\&le;1---\left(45\right)$

In formula (45), d is damage variable, δ _mthe maximal value reaching,

G. calculating elastic matrix.Elastic matrix is divided into 5 kinds of forms according to the difference of loading, is specially:

1. work as time, unmarred elastic stage:

2. work as and δ ₃>=0 o'clock, the damage softening stage: $D=\left[\begin{array}{ccc}(1-d)K& & \\ & (1-d)K& \\ & & (1-d)K\end{array}\right];$

3. work as and δ ₃during <0, the damage softening stage: $D=\left[\begin{array}{ccc}(1-d)K& & \\ & (1-d)K& \\ & & K\end{array}\right];$

4. work as and δ ₃>=0 o'clock, damage the stage completely:

5. work as and δ ₃during <0, damage the stage completely:

Step 5, calculates each combined interface unit joint forces.Described combined interface unit joint forces refers to the joint forces of major node, by formula (46), calculates:

In formula (46), F is the joint forces of each combined interface unit major node, and τ is stress, by formula (47), calculates

τ＝Dδ（47）

In formula (47), τ is stress, and D is elastic matrix, and δ is the local relative displacement in formula (39).

Step 6, calculates each combined interface element stiffness matrix.By formula (48), calculate the element stiffness matrix of each combined interface unit

In formula (48), K is the element stiffness matrix of each combined interface unit; T _δdisplacement transition matrix between the middle rigidity unit's major node of formula (35) and secondary node, t _δtransposed matrix; the transition matrix of using while removing rotational freedom in formula (40), be transposed matrix; B is the geometric matrix in formula (34), B ^tit is the transposed matrix of B; D is elastic matrix.

Step 7, assembles the global stiffness matrix of whole finite element model.Described assembling global stiffness matrix is before separating balance equation, by finite element software ABAQUS, the element stiffness matrix of whole finite element model is assembled into global stiffness matrix.

Step 8, solves balance equation.Adopt Numerical Iteration Method, by finite element software ABAQUS solver, solve described balance equation.In solving, after completing an iteration, whether the displacement that judges all nodes in finite element model restrains, if obtain the convergence solution of all nodal displacements in finite element model, this iteration step is complete, otherwise repeating step 2～step 8, until obtain the convergence solution of all nodal displacements in finite element model.

As shown in Figure 4, the present embodiment is tested by the double cantilever beam of T300/977-2 unidirectional fibre enhancement layer plywood, analyzes the interlayer damage of this T300/977-2 unidirectional fibre enhancement layer plywood.The finite element model of the double cantilever beam specimen of T300/977-2 unidirectional fibre enhancement layer plywood is set up by business finite element software ABAQUS, and in two beam arms, four Node Quadrilateral Finite Element plate shell unit S4R are set up in face position respectively; Except in the region of prefabricated initial crack, with the combined interface unit of limited thickness, connect upper and lower plate shell unit, record the load displacement course of loading end.In Fig. 6, plate shell unit discrete in neutral surface 15 simulations for upper beam arm, plate shell unit discrete in neutral surface 16 simulations for lower beam arm; Between upper beam arm 15 and lower beam arm 16, except precrack region, all use combined interface unit 17 to connect.

In the DCB finite element model of the present embodiment, unit size 0.5mm * 0.5mm, as shown in Figure 5, in Fig. 5, curve 13 is test findings for the loading end load-displacement curves calculating and the comparing result of trial value, curve 14 is result of calculation; Y axle represents load value, and x axle represents load deflection.Adopt combined interface of the present invention unit and finite element model to calculate analog result better, calculating curve overlaps with trial curve substantially at linearity range, descending branch after peak value is consistent with the trend of trial curve, and the simulation of interlayer damage initial sum damage evolutionary process is all realistic.The load peaks 62.8N calculating, trial value 63.1N, error-0.4%.

Claims (1)

1. a method for COMPOSITE MATERIALS laminate interlayer damage, is characterized in that, concrete steps are:

Step 1, sets up finite element model;

Step 2, extracts displacement of joint and the coordinate of each combined interface unit; The process of extracting comprises:

A. extract displacement of joint; From the finite element model having established, extract the displacement of each combined interface unit major node; This major node displacement is expressed as by vector form:

$u_M={(u_M^{1-},v_M^{1-},w_M^{1-},{\mathrm{\&theta;}}_{\mathrm{Mx}}^{1-},{\mathrm{\&theta;}}_{\mathrm{My}}^{1-}...,u_M^{4+},v_M^{4+},v_M^{4+},{\mathrm{\&theta;}}_{\mathrm{Mx}}^{4+},{\mathrm{\&theta;}}_{\mathrm{My}}^{4+})}_{40\&times;1}^T---\left(1\right)$

In formula (1), u _mthe major node displacement of combined interface unit, represent that respectively i major node is to along x, y, the translation displacement i=1 of z direction, 2,3,4; represent that respectively i major node is to around x, the rotation displacement of y axle, subscript negative sign represents the displacement component of lower major node, plus sige represents the displacement component of major node; Described major node is to referring to x in same combined interface unit, and y coordinate is identical, two major nodes that z coordinate is different;

B. extract node coordinate; From the finite element model having established, extract the coordinate of each combined interface unit major node; This major node coordinate is expressed as by vector form:

$x_M={(x_M^{1-},y_M^{1-},z_M^{1-},...,x_M^{4+},y_M^{4+},z_M^{4+})}_{24\&times;1}^T---\left(2\right)$

In formula, x _mthe major node coordinate of combined interface unit, respectively i the x that major node is right, y, z coordinate i=1,2,3,4, subscript negative sign represents lower major node coordinate components, plus sige represents major node coordinate components;

The node of described combined interface unit refers to the major node of combined interface unit; Described major node refers to the node that combined interface unit is connected with plate shell unit;

Step 3, calculates the geometric matrix of each combined interface unit; This process comprises the following steps:

A. calculate the right coordinate of the secondary node in each combined interface unit; According to extracting the coordinate of each combined interface unit major node obtaining in formula (2), by formula (3):

$x_S=\frac{1}{2}\left[I_{12\&times;12}\right|I_{12\&times;12}]x_M---\left(3\right)$

Obtain the right coordinate array x of the secondary node in combined interface unit _s:

$x_S={(s_S^1,y_S^1,z_S^1,...,x_S^4,y_S^4,z_S^4)}_{12\&times;1}^T---\left(4\right)$

In formula (3), x _sthe right coordinate array of secondary node of combined interface unit, I _{12 * 12}it is unit diagonal matrix; In formula (4), respectively i the right x of secondary node, y, z coordinate i=1,2,3,4; Described secondary node is in the secondary node that refers to inside, combined interface unit and comprise, a pair of secondary node that position overlaps mutually; Described secondary node refers to the node that rigidity unit is connected with cohesion unit;

B. calculate the shape function matrix of each combined interface unit; By formula (5), calculate shape function matrix

$H(\mathrm{\&xi;},\mathrm{\&eta;})=\left[\begin{array}{ccccccc}{N}_1& 0& 0& & N_4& 0& 0\\ 0& N_1& 0& |...|& 0& N_4& 0\\ 0& 0& N_1& & 0& 0& N_4\end{array}\right]---\left(5\right)$

In formula (5), H (ξ, η) is shape function matrix; N _ishape function i=1,2,3,4, be specially:

N ₁＝(1-ξ)(1-η)/4, N ₂＝(1+ξ)(1-η)/4

(6)

N ₃＝(1+ξ)(1+η)/4,N ₄＝(1-ξ)(1+η)/4

In formula (6), ξ and η are the coordinates under natural system of coordinates, ξ, η ∈ [1,1];

C. calculate the transition matrix of each combined interface unit; By formula (7) compute vector v _ξand v _ηx, y, z component:

$\begin{array}{c}{v}_{{\mathrm{\&xi;}}_i}=\underset{k=1}{\overset{4}{\mathrm{\&Sigma;}}}\frac{\&PartialD;N_k}{\&PartialD;\mathrm{\&xi;}}\&CenterDot;x_{\mathrm{Si}}^k\\ v_{{\mathrm{\&eta;}}_i}=\underset{k=1}{\overset{4}{\mathrm{\&Sigma;}}}\frac{\&PartialD;N_k}{\&PartialD;\mathrm{\&eta;}}\&CenterDot;x_{\mathrm{Si}}^k\end{array},i=\mathrm{1,2,3}---\left(7\right)$

By the vector v obtaining _ξand v _η, by formula (8) compute vector v _s, v _t, v _n

v _n＝(v _ξ×v _η)||v _ξ×v _η|| ^-1

v _s＝v _ξ||v _ξ|| ^-1 (8)

v _t＝v _n×v _s

By the vector v obtaining _s, v _t, v _n, by formula (9), calculate transition matrix

Θ＝(v _s,v _t,v _n) (9)

In formula (9), Θ is transition matrix; In formula (7), it is vector v _ξi component, i=1,2,3 corresponding x respectively, y, z coordinate; it is vector v _ηi component, i=1,2,3 corresponding x respectively, y, z coordinate; be k secondary node to the coordinate in i direction, by formula (3), obtain;

D. calculate the geometric matrix of each combined interface unit; By formula (10) computational geometry matrix

B＝Θ ^TH(ξ,η)Φ (10)

In formula (10), Θ ^tit is the transposed matrix of transition matrix Θ; H (ξ, η) is the shape function matrix in formula (5); Φ=(I _{12 * 12}| I _{12 * 12}), I _{12 * 12}it is unit diagonal matrix;

Step 4, calculates the elastic matrix of each combined interface unit; This process comprises the following steps:

A. calculate the secondary displacement of joint of each combined interface unit; According to formula (1), extract the displacement of each combined interface unit major node obtaining, by formula (11)

u _S＝T _δu _M (11)

Obtain the secondary displacement of joint of each combined interface unit:

$u_S={(u_S^{1-},v_S^{1-},w_S^{1-},{\mathrm{\&theta;}}_{\mathrm{Sx}}^{1-},{\mathrm{\&theta;}}_{\mathrm{Sy}}^{1-}...,u_S^{4+},v_S^{4+},w_S^{4+},{\mathrm{\&theta;}}_{\mathrm{Sx}}^{4+},{\mathrm{\&theta;}}_{\mathrm{Sy}}^{4+})}_{40\&times;1}^T---\left(12\right)$

In formula (12), u _sthe displacement of the secondary node of combined interface unit, represent that respectively i secondary node is to along x, y, the translation displacement i=1 of z direction, 2,3,4; represent that respectively i secondary node is to around x, the rotation displacement of y axle; Subscript negative sign represents the displacement component of lower secondary node, and plus sige represents the displacement component of secondary node; T _δbe displacement transition matrix between the major node of combined interface unit and secondary node, be specially

$T_{\mathrm{\&delta;}}={\left[\begin{array}{cccc}{t}_{{\mathrm{\&delta;}1}^{-}}& & & 0\\ & t_{{\mathrm{\&delta;}2}^{-}}& & \\ & & ...& \\ 0& & & t_{{\mathrm{\&delta;}4}^{+}}\end{array}\right]}_{40\&times;40}---\left(13\right)$

In formula (13) be displacement transition matrix between the right major node of i rigidity unit and secondary node (i=1 ..., 4), the subscript of i is that negative sign represents lower rigidity unit, plus sige represents rigidity unit; concrete form be

$t_{{\mathrm{\&delta;i}}^{\&PlusMinus;}}=\left[\begin{array}{ccccc}1& 0& 0& 0& z_S^{i\&PlusMinus;}-z_M^{i\&PlusMinus;}\\ 0& 1& 0& -(z_S^{i\&PlusMinus;}-z_M^{i\&PlusMinus;})& 0\\ 0& 0& 1& y_S^{i\&PlusMinus;}-y_M^{i\&PlusMinus;}& -(x_S^{i\&PlusMinus;}-x_M^{i\&PlusMinus;})\\ 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 1\end{array}\right]---\left(14\right)$

B. calculate local relative displacement; Described local relative displacement refers to the relative displacement under local coordinate system; By formula (15), calculate local relative displacement

In formula (15), δ is local relative displacement, and B is the geometric matrix in formula (10), it is the transition matrix of using while removing rotational freedom; The described transition matrix of using while removing rotational freedom is specially

${\mathrm{\&gamma;}}_{3\&times;5}\left[\begin{array}{ccccc}1& 0& 0& 0& 0\\ 0& 1& 0& 0& 0\\ 0& 0& 1& 0& 0\end{array}\right]---\left(16\right)$

C. calculate total relative displacement of mixed type interlayer damage; Described mixed type interlayer damage refers to and comprises I type simultaneously, the interlayer type of impairment of II type and the damage of III type; Described total relative displacement refers to mixed type relative storey displacement amount; By formula (17), calculate total relative displacement of mixed type interlayer damage

${\mathrm{\&delta;}}_m=\sqrt{{\mathrm{\&delta;}}_1^2+{\mathrm{\&delta;}}_2^2+{\&lang;{\mathrm{\&delta;}}_3\&rang;}^2}=\sqrt{{\mathrm{\&delta;}}_{\mathrm{shear}}^2+{\&lang;{\mathrm{\&delta;}}_3\&rang;}^2}---\left(17\right)$

In formula (17), δ _mit is total relative displacement of mixed type interlayer damage; δ ₁, δ ₂, δ ₃be local relative displacement δ respectively along the x under local coordinate system, y, three components of z direction; $\&lang;{\mathrm{\&delta;}}_3\&rang;=\left\{\begin{array}{c}{\mathrm{\&delta;}}_3,{\mathrm{\&delta;}}_3\&GreaterEqual;0\\ 0,{\mathrm{\&delta;}}_3<0\end{array}\right.;$

D. displacement calculating mixing ratio; Described displacement mixing ratio refers to the ratio of interlayer shear displacement and interlayer opening displacement, specifically

$\mathrm{\&beta;}=\frac{{\mathrm{\&delta;}}_{\mathrm{shear}}}{{\mathrm{\&delta;}}_3}---\left(18\right)$

In formula (18), β is displacement mixing ratio;

E. calculate mixed type interlayer and damage initial displacement and ultimate failure displacement; Described mixed type interlayer damages initial displacement and refers to that the damage of mixed type interlayer is initial when the total relative displacement of mixed type arrives this value; The displacement of described mixed type interlayer ultimate failure refers to that the damage of mixed type interlayer is developed when the total relative displacement of mixed type arrives this value; Specifically

${\mathrm{\&delta;}}_m^0=\left\{\begin{array}{cc}{\mathrm{\&delta;}}_3^0{\mathrm{\&delta;}}_1^0\sqrt{\frac{1+{\mathrm{\&beta;}}^2}{{\left({\mathrm{\&delta;}}_1^0\right)}^2+{\left(\mathrm{\&beta;}{\mathrm{\&delta;}}_3^0\right)}^2}},& {\mathrm{\&delta;}}_3>0\\ {\mathrm{\&delta;}}_{\mathrm{shear}}^0,& {\mathrm{\&delta;}}_3\&le;0\end{array}\right.---\left(19\right)$

${\mathrm{\&delta;}}_m^f=\left\{\begin{array}{cc}\frac{2}{K{\mathrm{\&delta;}}_m^0}[G_{\mathrm{IC}}+(G_{\mathrm{IIC}}-G_{\mathrm{IC}}){\left(\frac{{\mathrm{\&beta;}}^2}{1+{\mathrm{\&beta;}}^2}\right)}^{\mathrm{\&eta;}},& {\mathrm{\&delta;}}_3>0\\ {\mathrm{\&delta;}}_1^f,& {\mathrm{\&delta;}}_3\&le;0\end{array}\right.---\left(20\right)$

In formula (19), that mixed type interlayer damages initial displacement; single I, II, III type interlayer damages the initial displacement of corresponding damage; N, S, T is respectively corresponding boundary strength, and S=T; k is interface rigidity, is taken as K=10 ⁶; Described single I, II, the damage of III type interlayer refers to and only comprises I, II, the damage of III type is a kind of type of impairment wherein; In formula (20), it is the displacement of mixed type interlayer ultimate failure; G _iC, G _iIC, G _iIICrespectively I, II, the critical value of III type interlayer damage energy release rate; single I, II, III type interlayer damages corresponding ultimate failure displacement, ${\mathrm{\&delta;}}_3^f=2G_{\mathrm{IC}}/N,{\mathrm{\&delta;}}_1^f=2G_{\mathrm{IIC}}/S,{\mathrm{\&delta;}}_2^f=2G_{\mathrm{IIIC}}/T;$ B-K parameter η is by testing acquisition;

F. calculate damage variable; Described damage variable refers to for characterizing the amount of mixed type interlayer faulted condition; Described faulted condition refers to the state of injury severity score; By formula (21), calculate mixed type interlayer damage variable

$d=\frac{{\mathrm{\&delta;}}_m^f({\mathrm{\&delta;}}_m^{\max }-{\mathrm{\&delta;}}_m^0)}{{\mathrm{\&delta;}}_m^{\max }({\mathrm{\&delta;}}_m^f-{\mathrm{\&delta;}}_m^0)},0\&le;d\&le;1---\left(21\right)$

In formula (21), d is damage variable, δ _mthe maximal value reaching,

G. calculating elastic matrix; Elastic matrix is divided into 5 kinds of forms according to the difference of loading, is specially:

When time, unmarred elastic stage: $D=\left[\begin{array}{ccc}K& & \\ & K& \\ & & K\end{array}\right];$

When and δ ₃>=0 o'clock, the damage softening stage: $D=\left[\begin{array}{ccc}(1-d)K& & \\ & (1-d)K& \\ & & (1-d)K\end{array}\right];$

When and δ ₃during <0, the damage softening stage: $D=\left[\begin{array}{ccc}(1-d)K& & \\ & (1-d)K& \\ & & K\end{array}\right];$

When and δ ₃>=0 o'clock, damage the stage completely: $D=\left[\begin{array}{ccc}0& & \\ & 0& \\ & & 0\end{array}\right];$

When and δ ₃during <0, damage the stage completely: $D=\left[\begin{array}{ccc}0& & \\ & 0& \\ & & K\end{array}\right];$

Step 5, calculates each combined interface unit joint forces; Described combined interface unit joint forces refers to the joint forces of major node, by formula (22), calculates:

In formula (22), F is the joint forces of each combined interface unit major node, and τ is stress, by formula (23), calculates

τ＝Dδ (23)

In formula (23), τ is stress, and D is elastic matrix, and δ is the local relative displacement in formula (15); Step 6, calculates each combined interface element stiffness matrix; By formula (24), calculate the element stiffness matrix of each combined interface unit

In formula (24), K is the element stiffness matrix of each combined interface unit; T _δdisplacement transition matrix between the middle rigidity unit's major node of formula (11) and secondary node, t _δtransposed matrix; the transition matrix of using while removing rotational freedom in formula (16), be transposed matrix; B is the geometric matrix in formula (10), B ^tit is the transposed matrix of B; D is elastic matrix;

Step 7, assembles the global stiffness matrix of whole finite element model; Described assembling global stiffness matrix is before separating balance equation, by finite element software, the element stiffness matrix of whole finite element model is assembled into global stiffness matrix;

Step 8, solves balance equation; Adopt Numerical Iteration Method, by finite element software, solve described balance equation; In solving, after completing an iteration, whether the displacement that judges all nodes in finite element model restrains, if obtain the convergence solution of all nodal displacements in finite element model, this iteration step is complete, otherwise repeating step 2～step 8, until obtain the convergence solution of all nodal displacements in finite element model.