Fastening solutions on composite structures: model and verifications of contact with friction, progressive composite damage, and ductile metal damage

In this study, the fundamental steps required for the mechanical analysis of bolted laminated composite structures were revealed. Advanced methods were developed in Fortran and Python to implement nonlinearity to the composite material model using MARC/MENTAT finite element software. Friction and damage parameters for HTA/6376 CFRP material are verified using experimental data sources in the literature. The critical Cockroft–Latham ductile damage parameter of AISI 304L sheet material with its anisotropic properties, which is required for the design of plastically deformable metal components in the composite joint, is computed. Good agreement was obtained between the friction and progressive damage of CFRP composite, ductile damage of AISI 304L, and the external experimental references.


Introduction
The friction mechanism defines the stress distribution in bolted composite connections. So, it is important to choose a proper method for numerical analyses. The two distinct friction methods, continuous and discontinuous or the 'Stick-Slip' model exist in MARC [1]. Investigated the friction modeling on MSC MARC FEA software package. The experimental results [1] revealed that the discontinuous model fits the load-deflection curves well for net-fit and large bolt-hole clearance values. However, in a continuous model, fitting the experimental results could be possible by adjusting the 'relative sliding velocity' parameter in MARC. The value of this parameter should be less than the maximum value at which the instability begins. According to the experimental results [1], the coefficient of friction (COF) between the two laminates is the most effective parameter, since a large amount of friction force is expected to be produced here. The results showed that COF 0.7 represents better than COF 0.45 for joint displacement below 0.4 mm. However, COF 0.45 matches closer to the experimental curve above 0.4 mm of joint displacement. Therefore, the COF was recommended to be chosen between the two as 0.40-0.45 for numerical studies [1]. The result of [1] highlights the immense superiority of the use of the 'Stick-Slip' model in the finite element analysis (FEA) simulation of composite bolted joints. It is known from the literature that increasing the clamp load as high as possible improves joint stability. So, COF under the head of the bolt should be minimum to maximize the preload gained by fastening torque. Preload keeps the assembly parts together. A practical view of composite bolted connections in [2] shows that the conservation of preload is a critical design criterion for controlling the friction force between assembly parts. The 'relaxation' can cause the loss of preload and hence should be eliminated or minimized. The relaxation initiates by the end of the tightening session and mainly depends on the matrix material of the composite. The root cause of this can be the heat input produced during the assembly phase by friction because the temperature can cause drastic changes in the matrix behavior.
CFRP composite plates are susceptible to the loads, perpendicular to their surface. Dry tightening of bolts on the CFRP plate may cause micro damage around the hole, especially if the tightening is repeated. Therefore, it is better to avoid recursive tightening procedures for CFRP laminates. Sealants, sleeves, or washers can be used for protecting the bearing region against wear. Also, corrosion must be considered in CFRP-metal bolted joints.
Contact method is an important factor in FEA. The node-to-segment method is widely used in numerical analysis. However, a new method, segment-to-segment was first introduced in [2] promising more precise contact results. MARC/Mentat FEA software has an option in contact controls to use this method if needed.
Damage is the most critical parameter in the joining process. It is advantageous to know if any damage occurs inside the bolt and composite plates. An experimental study in [3] reveals the microscopic damage mechanism and shows the intra-ply damage with useful SEM images. Interference fit is known to increase mechanical performance of the joint. However, higher interference fit ratios can result in early damage around the bolt-hole region. The study in [4] reveals the damage mechanisms induced by bare bolt insertion. Methods in [5,6] can be applied to the post-process of FEA for composite damage analysis. A numerical study in [7] predicts failure load for two different tightening torque conditions. A numerical and experimental study in [8] clarifies the damage mechanism of woven fiber-reinforced composites by using ABAQUS commercial FEA software. A detailed experimental study in [9] reveals four important steps of damage, initiation, progression, non-linear softening, and catastrophic failure. A comprehensive study [10] covering the design and failure analysis of composite bolted joints for aerospace composites can be considered a good resource for fundamental issues.
In this study, a new methodology was developed for the composite verification process using Python and Fortran inside the MARC FE software. The introduced method was for the verification of the friction behavior of a bolted HTA/6376 CFRP composite connection, and for the determination of the critical damage indices of the composite. The Cockroft-Latham damage parameter for AISI 304L sheet metal is computed which does not exist in the literature. These parameters are the fundamental requirements for a reliable assembly analysis. The results were compared with references in the literature and found to be in good agreement. The proposed model in this study is advantageous where automated parametric design iterations of a model consisting of AISI 304L sheet metal forming process onto the HTA/6376 CFRP composite material, and are needed to be done, seamlessly.

Material and method
The analyses were carried out using Python. The flow chart of the main program is given in figure 1. The code structure is designed for multiple execution cycles to carry out parametric multi-numbered analyses. The main Python program listens to the open ports and initiates the subsequent analysis once the submitted analysis is completed. The input parameters of the FE model are stored in a Python function and executed sequentially. The most important function of the code is the mesher algorithm, which is written for the present study in Python. The FE mesh of the composite parts is refined and optimized according to the density function. The element and node numbers are stored before the execution of MARC job(s), to be used to generate the report file.
The nonlinearity in the shear behavior of the composite plate is successfully implemented in the composite material model and combined with the damage model of the composite material using a new set of codes that are newly written for the present study.
The FE studies are divided into three portions, contact and friction modeling, composite damage modeling, and metal damage modeling. HTA/6376 CFRP material was preferred for the first two steps of the verification analyses since it is currently used in the aerospace industry. AISI 304L stainless steel material is also a common material for metal forming applications, and that is chosen for the sheet metal forming application over the composite parts.

Contact and friction modeling
An experimental and numerical study for three-dimensional FEA of single-lap bolted joints in [11] is chosen for the first stage of the FEA verification procedure. The analysis models are generated following [11] and [1]. The models in [12] were implemented in Python to generate various configurations of mesh geometries. MARC internal mesher is not used for this case. Each of the 2D quad elements is generated, and expanded to 3D-hex, individually via Python code following figure 1. An overview of the FEA model is given in figure 2. According to [11] one of the composite plates is fixed at the end, and the other is pulled up to 5 kN of force in the long-edge direction of the composite plates. Strain and load values are read by using eight strain gauges attached to the composite plates in figure 3.   The joint geometry is designed with the ratios, w/d = 6, e/d = 3, d/t = 1.54, where w, e, d, and t represent the width of the lamina, distance from the center of the bolt hole to the edge, the diameter of the bolt hole, and the thickness of the plate, respectively. Quasi-isotropic HTA/6376 lay-up is used with the stacking sequence of [45/0/ − 45/90] 5s . Each ply has a nominal 0.13 mm thickness. So, the laminate thickness is around 5.2 mm when cured. The diameter of the bolt hole, the thickness of the washer, and the inner and outer diameter of the washer are designed as 8 mm, 1.6 mm, 15 mm, and 8.4 mm, respectively. Two clearance properties C1 = 0 µm, and C4 = 240 µm are analyzed for the verification of this study. In the experimental studies of [11], the bolts are torqued to 0.5 Nm (finger tight), which corresponds to 360 N of bolt pre-load or 7.2 MPa of bolt pre-stress. Bolt and nut are merged in the analysis. In the material properties of the washer, the thermal expansion coefficient is defined only in the thickness direction. External heat is applied as an initial condition to the nodes of the washer on top laminate. So, the washer elements will expand in the direction of thickness and produce a thermal prestress condition. The plasticity option is also activated for the washer elements to compensate for the excessive stresses during contact computations.
The purpose of this section is to verify a reliable friction model for a real-world scenario. In this regard, 'Continuous' and 'Stick-Slip' (Coulomb) friction models are compared in [1] and it is concluded that the 'Stick-Slip' friction model describes the friction condition better than the continuous model. In conjunction with the results in [1], a lower friction coefficient for the laminate-laminate interface gives closer results to the experiment. Therefore, kinetic friction coefficients for bolt-laminate, washer-laminate, and laminate-laminate are chosen as 0.1, 0.3, and 0.42, respectively.
The control parameters α, β, and e represent the friction coefficient multiplier (µ static /µ kinetic ), slip-to-stick transition region distance, and the convergence tolerance on the friction force, respectively. These parameters are set to α = 1, β = 10 −4 and e = 0.1. The single-sided method is used to describe the contact. Overall contact tolerance is set to 0.01, and the bias factor is defined as 0.9. Outside and inside contact tolerance values are set to (1 − 0.9) 0.01 = 0.001, and (1 + 0.9) 0.01 = 0.019, respectively. It is known that setting the contact bias parameters correctly, reduces the increment splitting because of penetration of contacting body node to the contacted body. MARC has the advantage of using continuous surface representation for 2D curves and 3D surfaces. It is known that discrete elements can cause uneven penetrations and can produce heterogeneous stress/strain distribution resulting in longer solution time and coarse outputs. Hence, smoothed boundary description option is activated for all contact bodies by setting the C 0 Continuity parameter. Linear or quadratic segments are replaced by coons surfaces with this option.
Element 149, a 3D, an eight-node composite brick element having three degrees of freedom (DOF) in each node was chosen for composite 3D bodies. It has four integration points on each layer and can unify up to 510 layers in one element. 20 layers are used for this case and each layer consists of 4 laminas. 17200 elements, and 11440 nodes in each laminate, and a total of 34400 elements, and 22880 nodes are used in the analysis. The minimum element edge length is 0.15 mm. Contacted elements of a laminate by bolt, washer, and the corresponding laminate were grouped and set as an individual contact body to speed up or ease contact detection, according to [11]. A simplified illustration of Element 149 [13] is given in figure 4. Element 7, an 8-node (3 DOF per node) hexahedral element type with 'full integration mode' was used for the outer volumes of the deformable bolt, and the nut bodies. 2560 elements, and 2880 nodes for the bolt, and 5760 elements, and 4800 nodes for washers are used. Element 136, a 6-node (3 DOF per node) pentahedral element with full integration mode, is used for the inner (core) volume of the bolt body. About 1120 elements and 1202 nodes are used for the bolt's inner core. Then the separately defined contact bodies which belong to a single physical body are tied together with the Glue option in contact properties.
The material properties of composite laminates and steel washers are given in table 1. Poisson's ratio υ 31 is derived from υ 13 for MARC input. The purpose of orthotropic material for steel washers is to set the thermal expansion coefficient only in the thickness direction, z. X5CrNiMo18_10_h (as rolled) material in the MARC material database for the washers. The main aim of this selection is to compensate for contact-induced peak stresses by small plastic deformations and is to eliminate singularity. This material has an isotropic hardening rule. The piecewise linear strain rate method is chosen for stress/strain approximation. The titanium bolt is modeled as an elastic-plastic isotropic material. Modulus of elasticity and Poisson's ratio are set to 110 GPa and 0.29, respectively.
For all load cases, the 'Non-Positive Definite' option is activated to increase accuracy, reduce rigid body motions, and improve the convergence ratio for user-defined time steps. In the second, longitudinal pulling load step, the face load is applied as a uniform and exponentialtime dependent pressure function to improve solution time and numerical stability.
a Temperature-dependent experimental data at T = 100 • C from MARC internal material database.

Composite damage modeling
Damage modeling is an important step for the main analysis in this study. However, precise accuracy is not needed. Damage results will be used as an indicative factor in further analyses to decide which design solution is viable and more reliable for HTA/6376 material. There are two studies selected from the literature for the verification analyses [14], and [15]. The first one presents a realistic solution for the longitudinal tensile case. An applicable model for implementing material nonlinearity is presented by using FE software. Since HTA/6376 material is used in the experimental and numerical studies of [14], the nonlinear portion of the material model is derived to be used in [15]. The source of the studies relies on [16] and tries to improve the relation between shear stress and strain. Equation (1) can give reasonable results for small strains. However, shear stress deviates from reality for larger strain values (figure 12) where, γ 12 , τ 12 , G 0 12 , and α are the in-plane shear strain, in-plane shear stress, interlaminar shear modulus, and the experimental coefficient, respectively.
ASTM standard D3518/D3518M-13 [17] explains the experimental setup and the calculation method of engineering shear strain. In conjunction with [17] and [14], 28.8 mm × 121.2 mm HTA/6376 lay-up with the stacking sequence of [±45] 4s was modeled by using 'Element 149' as given in figure 5. Totally, 4 layers (16 plies) and each layer contains 4 laminas with material properties given in table 2. 960 elements and 960 nodes are used for these analyses. The minimum element edge length is 1.4 mm.
The displacement constraint ( figure 5) is applied by a linear function from zero to 10 mm in a second. The fixed end constraints are applied as three individual boundary conditions. This method is similar to the well-known 3-2-1 method which aims to eliminate stress concentration nearby the fixed region and improve numerical accuracy. The three separately applied boundary conditions are given in figure 5.
According to the standard [17], the set of equation (2) is used to obtain G 12 , which is the input for MARC material property where, τ 12, i , P i , w, h, γ 12, i , ε x, i , ε y, i , τ 12 , G 12 , and γ 12 are shear stress at ith data point, force at i th data point, coupon width, coupon thickness, engineering shear strain at i th data point, longitudinal normal strain at ith data point, and lateral normal strain at i th data point, shear stress, shear modulus, and engineering shear strain, respectively. MARC/MENTAT has an option to utilize user-defined subroutine codes which the user can change or add features to the FEA. Python and Fortran are the two programming languages to meet this requirement. Even though modern Python has more options than Fortran, there still exist some components of MARC which is written in Fortran language. It is because Fortran is fast. Therefore, the Fortran language is used for this material modeling case. To implement nonlinearity, a new subroutine was compiled in Fortran. The flow chart is given in figure 6.
In the first step, G 12 is derived from the experimental data by using equation (2c), then the shear modulus vs shear strain curve is fitted to a 4th-degree polynomial function. This polynomial function is written as a software function in Fortran, then the program run as given in figure 6. In the flow chart, ubginc is the Fortran program that starts at the beginning of each increment. Then if the flag is set by the material properties, md_hooklw subroutine is processed after the ubginc. md_hooklw allows the user to manipulate material properties by changing the compliance or the stiffness matrix of orthotropic materials. On the other hand, an internal function elmvar allows the user to get variables for elements such as strains. These variables are stored in the memory for each layer and each integration point. Hence, by the implementation of the polynomial function, the shear modulus is changed in every increment to maintain the correlation between experimental data.
In the next step of damage modeling the main damage rules for the HTA/6376 material are applied from [15]. In this study, the numerical setup was designed as given in figure 7. A CFRP laminate and an AA7475-T76 aluminum plate are joined by a 6Al114VSTA titanium bolt and nut couple, then pulled in the direction shown in figure 7. The CFRP stacking sequence was [(0/ ± 45/90) 2 ] s . The coupon and the aluminum plate have the same width and height, 28.8 mm and 150 mm, respectively. The bolt hole diameter is 4.85 mm The thickness of the aluminum plate and the composite laminates are 4 mm and 4.16 mm, respectively. The distance from the center of the bolt hole to the edge is 14.4 mm.
The composite plate is numerically designed in 4 layers and each layer is made up of 4 laminas. The total number of elements in the model is 12480. 5056 elements and 3163 nodes (type 149) for composite laminate, 5056 elements (type 7) and 3163 nodes for aluminum plate, 1344 elements (type 7) and 725 nodes for bolt (outer), 640 elements (type 136) and 66 nodes  Light springs are 2D linked elements consisting of two nodes, which can stabilize the stiffness matrices without any force contribution. One node can be attached to an existing node, and the other can be attached to the ground or to another existing node. 3D contact bodies have rigid body modes before they come into contact. Therefore, a small stabilizer is attached to the critical locations. The value of these stabilizers is small as 10 N mm −1 . They are linked and added to the model as numerical stabilizers (figure 7). 10 N mm −1 of bolt torque resulting from 10.4 kN preload is divided by the surface area of the bolt shank and applied to the lower end surface of the bolt elements as a uniform pressure. The pressure load is applied in a linear form from zero to its max. value through the load case. During the preload case, the nut is fixed in the z direction at its lower end while the bolt is moving in the −z direction. The 'Follower force' option is activated to maintain the clamp load direction through the load case. Elements are set not to be released at the end of the load case. A new load case is added before the pulling stage to immediately glue the nut elements with the bolt to maintain clamping condition. The displacement is also applied as a linear function from zero to its max. value for numerical stabilization. The material properties for HTA 7/6376 are similar to table 2, except for G 12 and G 13 , which are nonlinear, and strain-dependent specifications as mentioned above. The endurance limits for the composite plate are given in table 3.
Damage criteria in [15] are summarized as follows: Matrix crushing failure: Matrix cracking failure: Fiber-matrix shearing failure: Fiber failure: If any of the criteria defined in equations (3)-(6) is satisfied, the material will degrade in the corresponding direction(s). Therefore, these criteria are defined in ufail subroutine code in accordance with table 4. Multiple equations of matrix crushing failure for compression and tension cases were combined in one failure index per case, f 1 and f 2 , respectively.
As mentioned in the explanation of md_hooklw, the compliance matrix can be changed if necessary. Since the shear modulus are nonlinear, and the resultant damage factors are not applied to all elasticity modulus (S 12 , S 13 , S 23 ), it is necessary to define a compliance matrix for all elements, layers, and integration points. Therefore, the new compliance matrix must be defined as follows: where, d 1 , d 2 , d 3 , d 4 , d 5 , and d 6 are the reduction factors (table 5) for material degradation.   The Fortran flow chart is given in figure 8 to explain how MARC will utilize the method. Some internal common blocks, such as concom, elemdata, dimen, array2, tbwbcs were included in the code to get the required variables, like the increment number, load case id, elemental data, variable dimensions, etc.

Metal damage modeling
Thin stainless-steel material is commonly used as an interface material for bolted composite structures. However, when it was used in a sheet metal form as a sleeve material, it is likely to be deformed by large strain values. Therefore, rupture may inevitably occur. For the feasibility of design studies, the damage must be considered. Othmen et al [18] investigated the ductile fracture of AISI 304L stainless steel sheet in stretching by FEA. In the essence of [18], the experimental cupping test results were compared with FEA to verify displacement-strain and displacement-force relations. Then the damage value in FEA was noted at the stroke when rupture begins in the experiment.
For this analysis, an identical 2D-axisymmetric type of an FEA model (figure 9) is generated in MARC/MENTAT. 'Element 10', a 4-node quadrilateral element with 2 DOF in the z and r directions per node, is used. The minimum element size is set to 0.1 mm, and 1753 elements with 2077 nodes are used in the analysis. The FE mesh of the sheet metal is refined from zero to 40 mm in diameter.
Even though there is a standard [19] ruling the test geometries, the parameters in [18] are used for modeling and analysis. The lower and upper dies are modeled identically, where they have 65 mm inner diameter with 5 mm fillet radius, and 120 mm outer diameter. The punch nose is modeled as a hemisphere having a radius of 60 mm. All three dies are modeled as rigid tools, and the sheet metal is deformable elastic-plastic as a workpiece. The friction coefficient between the punch and the workpiece is set to 0.2. The mechanical properties of AISI 304L material are set as E = 200 GPa, υ = 0.3 where, E and υ represent the modulus of elasticity, and the Poisson's ratio, respectively. The flow stress of AISI 304L material was defined in a table using the formulation given in equation (8) from ε p = 0 to the necking point ε p = 0.451, and the experimental tabular data from the necking point to the rupture in [18] where: σ 0 is the yield stress at zero plastic strain, ε p is the equivalent plastic strain, while Q and b represent the material parameters. The complete flow data [18] is given in figure 10, where the strain rate, . ε = 0.0012 s −1 . Anisotropic properties in [18] are entered into the MARC by using the Hill'48 anisotropic yield criterion [20]. Required values are given in table 6.
The simulation process is divided into two cases. In the first load case, the upper die is clamped to the specimen with 158 kN of axial force. Then the punch is introduced in the second load case by moving in the −z direction.
Apart from [18] Cockroft&Latham damage model is used and matched with the rupture strain. The damage model is given in equation (9) σ max σε dt ⩾ C where, σ max is the maximum principal stress,σ is the effective von Mises stress and . ε is the effective plastic strain rate. C is the material constant threshold for damage.

Verification of results
2.4.1. Verification 1: contact and friction modeling. The mesh model in [12] and the friction properties in [1] were combined with [11] using the new code written in Fortran. The FE model was designed using a moderate mesh density. The results C1-Exp, C1-Base, C1-Imprv, and  figure 3, and represent the experimental [11], the base (coarse) mesh [11], improved (fine mesh) [11], and the present model results for C1 (0 µm) clearance, respectively. The C4-MI and C4-FEA-Present in figure 3 indicate the results for C4 (240 µm) clearance of improved (fine mesh) [11] and the present model, respectively. According to the comparative plots in figure 11, strain values of gauges 6, and 7, which are parallel and perpendicular to the loading direction (figure 3), fit the best to the experiment in all numerical analyses. The other results are close enough to be used as a reference for further analyses in this study.

Verification 2:
composite damage modeling. The most significant factor in the tensile load-displacement behavior of a composite plate is the nonlinearity of the shear modulus. The proper parameters were obtained after fine-tuning the initial polynomial function by trial and error. The results given in figure 12, where the curves represented by α are derived by using equation (1), are detailed in the previous section.
The results closely fit the experimental data in figure 12, which means that the current model perfectly represents material nonlinearity. Therefore, the results obtained here were applied to progressive damage analysis studies.
For ply-by-ply damage representation, advanced methods were applied to the MARC using PLOTV Fortran subroutine and a new Python script. In this way, even though each of the elements consists of 4 composite layers, and even the MARC allows the user to plot one failure index at a time, it is now possible to plot every layer, individually. It should be noted that each of the failure indices plotted in figure 13 consists of multiple layers (more than 4) in different orientation angles. In other words, the simulation model has 4 layers for composite ply, however, the plotted results have 16 layers.
According to figure 13, the matrix material is the first and the most damaged component in the assembly. The analysis was stopped when numerical instability started due to intensive and accumulated material degradation. The fiber and fiber-matrix damage patterns 3, 4, and 5   in figure 13 are in good agreement with the results presented in [16]. It is obvious in figure 13 that for the first and the last layers of the composite plate, the bolt pretension and the secondary bending cause the fiber to fail nearby the bolt seating surface.
It should be considered that the secondary bending effect becomes more significant when the plate gets thinner. The displacement plot for the present configuration is given in figure 14.
However, the destructive effect of secondary bending can be reduced to a point by making the stress distribution as homogeneous as possible. At this point, additional parts like insert, washer, or sleeve can be used. This topic is still open to researchers for years hiding great opportunities for composite structures.
The tensile load-displacement curve of the experiment is another indicator to verify the analysis. Based on the results of load-displacement curves in figure 15, the present study represents the worst-case scenario for damage, which is much better for safety. In other words, the present study follows the lower load path of the experiment.

Verification 3: metal damage modeling.
The verification of ductile damage of AISI 304L stainless-steel material relies on the consistency of the load-displacement curve of the punch, major strain, and the maximum total strain of the sheet metal. However, punch force is susceptible to the clamping region where the jaws are holding the sheet metal specimen. As mentioned above, the COF was set to 0.2 for the punch-specimen interface. Nevertheless, to reveal the effect of COF within the clamping region, after the prestressing load case, analyses have been repeated for four COF values, 0.2, 0.8, 0.86, and 0.9, respectively. On the other hand, additional analysis covering a glued contact between the specimen and the jaws was computed. The results in figure 16(a) showed that the punch force gets closer to the experiment as the COF gets higher values. Similarly, the maximum value of the major strain, which has a significant effect on the formability of the sheet metals, converged to the maximum value of the experimental total strain value when the COF increased ( figure 16(b)). Moreover, the maximum of the major strain occurred close to the rupture stoke at displacement x = 25 mm. As for the pattern of experimental total strain history, it perfectly matches the glued condition if the results are divided by 2.00 in figure 16(c). In addition, the mean value of the stress  triaxiality ratio from 0.2 to the rupture strain in the Erichsen test [18] was 0.64, while it is found to be 0.65 in this study. In this respect, the results are consistent with [18].
Experimental 2D plots are given in figure 17 to compare the major strain, max principal strain, and damage values with the experiment. According to figure 17, the major and maximum principal stresses are concentrated in regions similar to those in the experiment.
The history of damage parameters until the rupture is given in figure 18. The Cockroft&Latham damage curve is obtained from the FEA in the present study, while the other curves were extracted from [18]. The critical value of the Cockroft-Latham damage index is found to be 0.52, where the major strain is 0.32 which is close to the experiment [18].
Damage parameters of composite and stainless steel materials obtained above will be monitored as an indicator of damage in future studies.

Conclusion
The present study describes the fundamental analyses of bolted composite joints using the most common materials in the aerospace industry. A new set of codes were developed in the present study using Python and Fortran for batch FEA modeler for autonomous parametric studies, advanced mesh generation for composites, implementation of non-linear orthotropic material model into the FEA, and progressive composite damage modeling. An advanced Fortran code was developed under the postprocessing segment for damage visualization. The two similar models in [14] and [15] for progressive damage were combined and the results of the new method are compared with the experiment in [15]. The critical Cockroft-Latham damage index for AISI 304L sheet metal material was obtained using FEA. Apart from [18], the strain pattern of the experiment was compared with the major and the maximum principal strain, and it was observed that the maximum principal strain pattern is close to the experimental results. The results of each analysis are verified by experimental and numerical data sources in the literature. Therefore, the material models and the simulation parameters can be used in the related design studies. In the continuing study, a new assembly method for bolted composite structures consists of two composite plates, a metal insert, bolt and the nut, involving plastic deformation of the metal insert (AISI 304L) and damage results of all materials in the assembly, will be analyzed and discussed based on the present study.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).