## 3 Finite Element Modeling  
The approach based on the finite element method is the most suitable tool and is capable of simulating the failure process of the composite structure under complex loads. The material property degradation method using CDM and the cohesive zone method are implemented in FE modeling to characterize the initiation, evolution, and interaction of the damage modes for the failure analysis of composite structures. A meso-level FE model is developed in the commercial FE package ANSYS to predict the progressive failure of intralaminar damage, interlaminar damage, and their interaction with the laminate composites. Composite layers are meshed with structural solid elements, and the interface between two layers is modeled with the cohesive zone elements in the ANSYS element library. The approach based on the coupling of continuum damage mechanics and the cohesive zone method is implemented by using ANSYS Parameter Design Language (APDL). The linear orthotropic material property is used before damage initiation, while the anisotropic elastic (ANEL) modulus matrix is used during damage evolution. Damage variables, crack direction parameters, the transform matrix, and the damaged stiffness matrix are all calculated using the element table and updated with the time step after the solution. Matrix multiplication and invert operations are performed using ANSYS matrix process functions. When simulating FPR composites failure, Hashin criteria and equivalent damage scale-based criteria are used for judging the damage initiation and evolution, respectively. Once a layer element is damaged, the damage variable and crack direction parameter are calculated at this element, and then, stiffness degradation is performed according to the damage evolution model based on the CDM theory. This stiffness degradation of layer elements leads to the stress redistribution and stress concentration on neighboring layer elements. The damage variable calculations are performed at all damaged elements under this stress concentration, and this process is repeated continuously under the same load level until there is no variation of the damage variable distribution of the element. Comparison of the damage variable distribution between two sequence processes is performed to guarantee the gradual increase of the damage variable at every element for every process. The concept of the newly-damaged element set is applied to the calculation of the crack direction parameter. The newly-damaged element set contains the element in which the damage variable at the current load step is greater than the one at the previous load step, and the damage variable at the previous load step is equal to zero. That is, in all elements of this set, the damage variables would have variations at the current time step, and the damage variables are equal to zero at the previous time step. To save the computational cost, the calculation of the crack direction parameter is only performed in the newly-damaged element set. On the other hand, the calculation of the damage level is also performed for the CZEs between two layers to predict the delamination evolution. Once the intralaminar elements adjacent to a CZE have failed, the degradation of critical energy release rates is performed according to Equation (15) for the CZEs. Then, the load increases, and the iterative process repeats until a final catastrophic failure occurs in the composite laminates. Convergence difficulties are often seen in numerical simulation for fracture problems of composites by using the material property degradation method of CDM or the traction-separation law of CZM. Although the energy stabilization method helps enhance the numerical convergence, it could not completely overcome the convergence difficulties. To alleviate these difficulties, a viscous damage variable is introduced in this work.