Defects and damages in materials which will lead to a discontinuity in material properties are inevitable. In general,materials initially may be recognized as homogeneous,but will become discontinuous under various external environments and loading histories due to the micro-cracks,micro-holes,and local damages within,leading to a localized degradation in the strength and stiffness^{[1, 2, 3, 4]}. To make full use of materials and improve economic efficiency,limited size damages are often allowable in structural designs. Therefore,it is very necessary to develop a numerical method for effectively and practically introducing damages with different extents and distributions in implementing a finite element simulation so as to study the behavior of the damaged structure efficiently and assess its load-carrying capacity with the damage evolution.

During a continuous loading process,damages may initiate,propagate,and evolve in a structure,leading to significant changes in the material properties. During a loading process,the structure progresses through different damage states,resulting in changes in the material properties. The residual strength of composite materials subjected to impact is one of the typical examples^{[5, 6]}. When implementing a numerical simulation,the continuous loading process is usually divided into several time steps for the purpose of achieving convergence and capturing changes in material properties and/or boundary conditions. Thus,for a multi-step approach, damages of the previous step and their influence on the current step must be taken into account. Therefore,an effective numerical approach for damage information transfer from one step to the next under a multi-step loading process is essential. Furthermore,in the finite element method (FEM),one usually selects a geometric part and assigns a material property. Obviously,this approach will not be suitable for the case when material properties continuously vary with spatial coordinates and time. Of course,material properties can be assigned to each element in the model,but the complexity of this assigning task makes it unfeasible to implement,especially for models with a large number of elements. As pointed in Ref. ^[7] that the assignment of material properties in FEM remains an unresolved issue.

Efforts have been made in recent years for considering spatial changes in material properties resulting from damages existing randomly in a structure^{[8, 9, 10]} . In the biomechanical field,a breakthrough has been achieved in obtaining non-uniform properties of human bone^[8] by using a data processing technique for CT images,and then pixels within the scanned image are extracted. By assuming a relationship between the pixels and the material properties,the material constants at specific locations are obtained. These constants are then assigned directly to the integration points of each element by means of an ABAQUS subroutine. The change of material properties was also described with the spatial coordinates through a smooth function in Refs. ^{[9, 10]}. However,this approach is difficult to implement in the commercial software because it has to modify the finite element formulation and compile related code to introduce the distribution of the material properties.

The main drawback of the above mentioned methods is its dependence on the assumptions for the relationship between material properties and the image pixels or the spatial coordinates, especially for a complex structure. It can be even more difficult to determine the material variation within the structure. Adding further complexity,for engineering structures,the material properties may degrade with time. The degradation of material properties may depend on the load history and working environments,and it is hard to assess them accurately in advance.

Although there is no direct way to solve these modeling challenges,one can get inspiration from the analytical procedure of multi-step loading problems which can be divided into two categories,i.e.,multi-step analysis method and integrated analysis method. The multi-step analysis method needs to assume the damage zone,modify the geometric model accordingly, and then calculate the response under loading. For example,when calculating the residual strength of a composite structure after impact,it is assumed that the damage area after impact is an elliptical inclusion^[11] or a circular area^{[12, 13, 14, 15]}. This method is based on the damages that happened in the previous step. Due to the assumptions made for the damage area,only the main damages can be included,and some useful damage data may be overlooked. The numerical model also has to be updated before proceeding to the next step calculation. The integrated analysis method,on the other hand,has been demonstrated to be capable of overcoming those shortcomings and being effective in data transfer^{[16, 17, 18, 19, 20]}. However,the skills in finite element programming are required for the integrated analysis method. Difficulty may exist in achieving a joint solution of implicit and explicit analyses. An additional challenge occurs when several types of damages exist simultaneously in a structure when using the integrated analysis method.

By combining the virtues of the multi-step analysis method and the integrated analysis method,an effective numerical approach is proposed in this paper for the analysis of damaged structures,especially for structures where single or multiple damages with different extents and locations may exist. Based on the theory of damage mechanics and by introducing field variables,the evolution of the damage state of the structure is represented by field-variable functions describing the degradation in material properties. Then,through customization of the commercial finite element software,the field variables are directly assigned to the integration points of elements,thereby avoiding modification of the finite element model. Finally, the proposed numerical approach is demonstrated by a delamination analysis of a laminated composite beam. The proposed approach can be used to simulate the behavior of a damaged structure and to assess the effects of different types of damages,such as matrix damage,fiber fracture,and delamination on the behavior of the structure. 2 Analysis procedure and damage evolution description

Commercial finite element software provides users with a large number of element types, material models,and effective solution methods capable of handling various problems. Due to the diversity of engineering problems,a commercial finite element analysis (FEA) software also provides users a standard interface to customize the tool according to their needs. The general process of the secondary development of finite element commercial software for general structures is shown in Fig. 1(a). To simulate a damaged structure,it is necessary to prepare the initial damage data before simulation and update the data during the solving process,as shown in Fig. 1(b).

Fig. 1 Procedure of secondary development of finite element commercial software

Figure options

It would be noted that material properties of a damaged structure will degrade with time when the damages develop progressively. Therefore,based on damage mechanics,material properties can be defined as a function of damage related variables,i.e.,

Click Browse formula

where X = (x₁,x₂,x₃,· · · ) represents the damage related variables,and Y stands for the material properties such as Young’s modulus and density.

There are two types of damage descriptions,i.e.,ladder-like damage and progressive damage. Material properties with ladder-like damage describe a discontinuous damage state,while progressive damage corresponds to a smooth curve,as shown in Fig. 2. The material constant for ladder-like damage can be obtained from the discrete pairs of data (E_i,X_i) (i=1,2,· · · ); while for progressive damage,the material constant can be discretized by choosing a series of points (E_i,X_i),see Fig. 2(b).

Fig. 2 Description of material property for two types of damages

Figure options

To analyze and evaluate the performance of a progressively damaged structure,three key problems have to be resolved: acquiring the information on the distribution of the structural damages; inputting the damage information into the model; and simulating a related process for the damaged structure.

Currently,there are three methods in acquiring the information on the distribution of the structural damages: (i) direct measurement,(ii) image data processing^{[8, 9, 10]},and (iii) finite element calculation. The damage information obtained by direct measurement is limited to the measuring means. On the other hand,the relationship between damage information and CT image pixels is not clear,and it is also difficult to describe the internal damages. Therefore, the third method is used in this study.

Variables which are used to describe a state are called field variables. A set of field variables is defined to describe the damage state. For example,three failure modes are considered for a laminated composite plate under the compression loading,i.e.,matrix failure,delamination, and fiber fracture. In the simulation process,field variables F_V1 ,F_V2 ,and F_V3 represent the matrix damage,delamination,and fiber damage,respectively. Their magnitude and variation with time will be determined first,and then be assigned to the integration points of an element. These variables may range in the value from 0 to 1^[21]. The value 0 represents the state of no damage,and 1 represents the state of failure. The material properties of a damaged structure are related to the field variables. Figure 3 schematically depicts the material properties for the two states,where E₁ and E₂ are elastic modules,and T₁ and T₂ are material strengths.

Fig. 3 Material properties for two states of damages

Figure options

Field variables will be updated with the damage evolutions. For example,in a composite plate after being impacted,field variables will be updated by Eq. (2) so as to obtain the new damage distribution,and determine the residual strength of the plate,i.e.,

Click Browse formula

where D_m,D_f,and D_del are the damage variables and denote matrix damage,fiber damage,and delamination,respectively. The field variables can be determined based on the Hashin failure criterion or Tsai-Wu failure criterion. Once the field variables are determined,the damage information between the steps is transferred.

The key point in the proposed approach is the transfer of damage data from one step to the other during the numerical simulation for a damaged structure. The approach mainly contains two steps and the flow chart of this approach is shown in Fig. 4.

Fig. 4 Flow chart of proposed approach

Figure options

Step 1 To generate the data file of the initial damage state. The damage distribution is obtained from the structure response under the initial load by using the FEA software ABAQUS/Standard. Then,Python language is used to establish a data file for the damage state.

Step 2 To define related field variables and to describe degraded material properties. ABAQUS/Standard reads the damage data file when starting an analysis through a user subroutine “UEXTERNALDB”,and implements the job according to user’s demands,such as compression or tension after impact. The user subroutine “UEXTERNALDB” is used to manage user-defined external databases and calculate the model-independent history information in ABAQUS/Standard^[22].

For the structures existing damage,only Step 2 is necessary,and the material properties in the damaged area can be assigned as being completely damaged.

The proposed approach is applicable for simulating a damaged structure with single or multi damages with various extent and locations on the structure,or a structure with initial imperfection,both resulting in progressive degradation of material properties. During the analysis process,the damage information can be preserved completely so as to avoid making assumptions for the damage state and changes for the finite element model. The most advantageous feature of this approach is that it can easily introduce partial or all damage information. Combining with the related criteria,the approach can also be used to assess the effects of different types of damages. 3 Verification and applications 3.1 Verification of proposed approach

When subjected to impact loading,even with low impact energy,laminated composites can be easily damaged,resulting in degradation of material properties. Delamination is one of the frequently occurred damages. To validate the proposed method,a double cantilever beam (DCB)^[23] (as shown in Fig. 5) with pre-delamination under Mode I is analyzed.

Fig. 5 Double cantilever beam: (a) Model I: initial crack is 55 mm in length; (b) Model II: initial crack is modeled by setting material property with full degradation

Figure options

The composite beam consists of T300/977-2 carbon fiber reinforced epoxy lamina laid in sequence (0^◦)₂₄. The material parameters are^[22] E₁₁ = 150.0 GPa,E₂₂ = E₃₃ = 11.0 GPa, G₁₂ = G₁₃ = 6.0 GPa,G₂₃ = 3.7 GPa,V₁₂ = 0.25,V₂₃ = 0.45,GIC = 0.268 kJ/m²,and T = 45 MPa,where GIC is the critical fracture toughness of Mode I. The length and width of the beam is 150 mm and 20 mm,respectively. The thickness is 3.96 mm. The initial crack length is 55 mm lying in the middle layer of the beam. A pair of pull forces is imposed at the free end. Cohesive elements are set in the middle layer of the beam. The constitutive relationship of the cohesive elements is shown in Fig. 6,where the relationship between the traction and the separation is described. k₀ is the initial stiffness,k₁ and k₂ are residual stiffness for different damage states.

Fig. 6 Constitutive relationship of cohesive element

Figure options

The initial delamination at the free end can be represented by two models. One is the same as that given in Ref. ^[23] by setting an initial crack with 55 mm in the length at the free end. The other is by setting a layer of cohesive elements in the pre-delamination area of the central layer of the DCB. Degradation of the cohesive elements can be expressed through a damage variable d,which takes values from 0 to 1. The relationship between the degraded stiffness k₁ and the initial stiffness k₀ is given as

Click Browse formula

The damage variable d in ABAQUS is represented by the field variable F_V1 . It takes 0 when no damage and 1 when the material is completely damaged.

As an example,the material properties with no damage are K = K₀,G = G₀,and F_V1 = 0, while the material properties when the material is completely damaged are K = 0,G = 0,and F_V1 = 1.

Figure 7 gives the stress distributions of the laminated composite beam of two numerical models under Mode I loading. The maximum stresses in Figs. 7(a) and 7(b) are 211.7 MPa and 212.7 MPa,respectively,and the difference is about 0.4%. Figure 8 shows the forcedisplacement curves and the experimental results^[23]. It is seen clearly that the two curves from the two models are consistent and they are in good agreement with the experimental results, showing that the structure with pre-crack can be modeled either by presetting the damaged region or by specifying the material state of the damaged areas. The proposed method is simple, effective,and easy to investigate the effects of the damaged areas and their locations on the performance of the structure.

Fig. 7 Stress distribution of laminated composite beam of two numerical models under Mode I loading

Figure options

Fig. 8 Comparison of force-displacement curves of delaminated beam (curves of Model I and Model II are almost coincident)

Figure options

The most outstanding feature of the proposed approach is that it can easily introduce partial or full damage information,combined with mechanical criteria,to assess the consequence resulted from partial or full damage. In this section,the effect of the impact damage on overall performance of the structure will be investigated by using the proposed method to demonstrate the advantages of the method. 3.2.1 Simulation of laminated composite under impact

Consider a laminated composite plate T700/3234^[24] with stacking sequence [0/90^◦/90^◦/0]. Its size is 100 mm× 76.2 mm× 0.8 mm. The plate is impacted by a cylindrical projectile with a spherical head and the projectile weights 6.17 kg,as shown in Fig. 9. The impact energy is 24 J. The Hashin criterion^[25] is used to assess the damages of the fiber and matrix,and the Yeh and Kim criterion^[26] is used to judge the delamination. After degradation of the material, Young’s modulus is assumed to be E_i′=0.3E_i.

Fig. 9 Laminate composite plate under impact

Figure options

During the analysis on the impact damage of the laminated composite,fiber fracture,matrix cracks,and delamination are described by three variables,i.e.,F_V1 ,F_V2 ,and F_V3,respectively. The value 0 or 1 of these variables will represent a failure state or a state without damage, respectively. Implementing the numerical analysis through the proposed approach,the maximum impact forces are obtained and compared in Table 1. By comparing with its experimental data given in Ref. ^[24],it can be seen that the prediction of the maximum impact force by this approach agrees better.

Table 1. Comparison of maximum impact force

Table options

Damages such as fiber fracture,matrix cracks,and delamination can be described by three variables,i.e.,F_V1 ,F_V2 ,and F_V3. The value 0 or 1 represents the state without damage or with a failure state,respectively. The failure state and corresponding field variables are listed in Table 2.

Table 2. Field variables and corresponding material properties

Table options

The first five modes of the laminated plate without damage and the laminated plate containing impact damage are obtained through ABAQUS,as shown in Fig. 10. The column (a) is the modes of laminated plates without damage,while columns (b) and (c) show the modes for the damaged laminates with damaged areas with different impact energy. Little difference can be found for the first three modes,but obvious differences can be seen for the fourth and fifth modes,indicating that the impact damage can cause changes in higher order modes of the structure.

Fig. 10 First five modes of plate with and without damage

Figure options

Table 3 lists the first five natural frequencies of laminates with and without damage. With the increase of the impact energy,the damage zone becomes larger,and the natural frequency of the structure decreases. The phenomena agrees well with those described in Refs. ^{[27, 28]}. The natural frequencies of the laminated plate can be obtained either directly from the related finite element model or through setting the field variables,F_Vi ,correspondingly. The natural frequencies obtained from both methods are consistent,showing that the method by setting the field variables is correct and effective.

Table 3. First five-order natural frequencies of laminates with and without damage

Table options

In order to see clearly the influence of damage on natural frequencies,the natural frequencies under different damage conditions normalized by the corresponding natural frequencies of the plates without damage,which is defined as the nominal frequency,are given in Table 4. In most cases,low-order frequencies are important. Therefore,the effects of the impact damage on the first five frequencies are considered. It may be observed that the first-order frequency decreases by 10% and the second-order frequency reduces by 12% for the impact energy of 9.4 J. When the impact energy increases to 24 J,the first-order frequency decreases by 11%,and the second-order decreases by 14%. Thus,the variation of the natural frequencies resulted from the impact damage would be paid attention in the structural design.

Table 4. Nominal frequencies

Table options

In addition,the proposed approach is also convenient to investigate the effects of the different damages. The effects of the three types of damages by defining the percentage reduction of the natural frequencies are shown in Table 5 under the impact energy of 24 J. As can be observed in Table 5 and Fig. 11,the matrix damage occurs over a larger area,thus leading to a greater decrease in nominal frequencies. The assessment of the effects of structural damages on the natural frequency of laminated structures has a great significance in the process of material selection and design of the laminates.

Table 5. Percentage of reduction of natural frequencies

Table options

Fig. 11 Comparison of damage area or extent of third layer for laminated plate

Figure options

By defining a set of field variables to describe damages,an effective numerical approach for damaged structure under various loading conditions is proposed. The effectiveness of the approach is demonstrated by simulating the delamination development under Mode I loading for a pre-delaminated cantilever beam and analyzing the impact damage of a laminated composite. The main advantages of the proposed approach can be described as follows.

(i) It is applicable for simulating the performance of a structure with single or multiple damages with various extents. It is also applicable in analyzing a damaged structure with different types of damages.

(ii) In this approach,the damage in a structure can be described conveniently through the directly assigning values of the field variables to the integration points of an element and can be updated during the solving process,avoiding making assumptions for the damage and changing the finite element model during the numerical simulation.