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

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{\text{ 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 ${\left[ {45/0/ - 45/90} \right]_{5s}}$. Each ply has a nominal $0.13{\text{ mm}}$ thickness. So, the laminate thickness is around $5.2{\text{ 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{\text{ mm}}$, $1.6{\text{ mm}}$, $15{\text{ mm}}$, and $8.4{\text{ mm}}$, respectively. Two clearance properties $C1 = 0{\text{ }}\mu {\text{m}}$, and $C4 = 240{\text{ }}\mu {\text{m}}$ are analyzed for the verification of this study. In the experimental studies of [11], the bolts are torqued to $0.5{\text{ Nm}}$ (finger tight), which corresponds to $360{\text{ N}}$ of bolt pre-load or $7.2{\text{ 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 $\alpha$, $\beta$, and $e$ represent the friction coefficient multiplier $\left( {{\mu _{{\text{static}}}}/{\mu _{{\text{kinetic}}}}} \right)$, slip-to-stick transition region distance, and the convergence tolerance on the friction force, respectively. These parameters are set to $\alpha = 1$, $\beta = {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 $\left( {1 - 0.9} \right)0.01 = 0.001$, and $\left( {1 + 0.9} \right)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}{\text{ 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{\text{ 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.

Zoom In Zoom Out Reset image size

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 ${\upsilon _{31}}$ is derived from ${\upsilon _{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{\text{ 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 exponential-time dependent pressure function to improve solution time and numerical stability.

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)

$$\begin{align}{\gamma _{12}} = \frac{{{\tau _{12}}}}{{G_{12}^0}} + \alpha \tau _{12}^3\end{align} \tag{ 1 }$$

where, ${\gamma _{12}}$, ${\tau _{12}}$, $G_{12}^0$, and $\alpha$ 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{\text{ mm }}\, \times 121.2{\text{ mm}}$ HTA/6376 lay-up with the stacking sequence of ${\left[ { \pm 45} \right]_{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{\text{ mm}}$.

Zoom In Zoom Out Reset image size

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

$$\begin{align}{\tau _{12,\,i}} = \frac{{{P_i}}}{{2\left( {w\,h} \right)}}\,\,\,\,\,\,\,(a)\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,{\gamma _{12,\,i}} = {\varepsilon _{x,\,i}} - {\varepsilon _{y,\,i}}\,\,\,\,\,\,\,\,(b)\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,{\tau _{12}} = {G_{12}}\,{\gamma _{12}}\,\,\,\,\,\,\,(c)\end{align} \tag{ 2 }$$

where, ${\tau _{12,{\text{ }}i}}$, ${P_{{\text{ }}i}}$, $w$, $h$, ${\gamma _{12,{\text{ }}i}}$, ${\varepsilon _{x,{\text{ }}i}}$, ${\varepsilon _{y,{\text{ }}i}}$, ${\tau _{12}}$, ${G_{12}}$, and ${\gamma _{12}}$ are shear stress at $i{\text{th}}$ 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.

Zoom In Zoom Out Reset image size

In the first step, ${G_{12}}$ is derived from the experimental data by using equation (2 c), 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 ${\left[ {{{\left( {0/ \pm 45/90} \right)}_2}} \right]_s}$. The coupon and the aluminum plate have the same width and height, $28.8{\text{ mm}}$ and $150{\text{ mm}}$, respectively. The bolt hole diameter is $4.85{\text{ mm}}$ The thickness of the aluminum plate and the composite laminates are $4{\text{ mm}}$ and $4.16{\text{ mm}}$, respectively. The distance from the center of the bolt hole to the edge is $14.4{\text{ mm}}$.

Zoom In Zoom Out Reset image size

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 for bolt (core), and $384$ elements (type 7) and 449 nodes for nut are used in the analyses. The minimum element edge length around the hole is $0.23{\text{ mm}}$. Friction parameters are defined as mentioned earlier. Fixed constraints in figure 7 are applied as explained in figure 5.

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⁻¹. They are linked and added to the model as numerical stabilizers (figure 7). 10 N mm⁻¹ 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:

$$\begin{align}{\sigma _2} < 0 \Rightarrow \,\,\,{d_{{\text{mc}}2}}& = \sqrt {{{\left( {\frac{{{\sigma _2}}}{{{Y_{\text{C}}}}}} \right)}^2} + {{\left( {\frac{{{\tau _{12}}}}{{{S_{12}}}}} \right)}^2} + {{\left( {\frac{{{\tau _{23}}}}{{{S_{23}}}}} \right)}^2}} \,\,(a)\,\,\,\,\,\nonumber\\ {\sigma _3} < 0 \Rightarrow \,\,\,{d_{{\text{mc}}3}}& = \sqrt {{{\left( {\frac{{{\sigma _3}}}{{{Z_{\text{C}}}}}} \right)}^2} + {{\left( {\frac{{{\tau _{13}}}}{{{S_{13}}}}} \right)}^2} + {{\left( {\frac{{{\tau _{23}}}}{{{S_{23}}}}} \right)}^2}} \,\,(b).\end{align} \tag{ 3 }$$

Matrix cracking failure:

$$\begin{align}{\sigma _2} > 0 \Rightarrow \,\,\,{d_{{\text{mt}}2}}& = \sqrt {{{\left( {\frac{{{\sigma _2}}}{{{Y_{\text{T}}}}}} \right)}^2} + {{\left( {\frac{{{\tau _{12}}}}{{{S_{12}}}}} \right)}^2} + {{\left( {\frac{{{\tau _{23}}}}{{{S_{23}}}}} \right)}^2}} \,\,(a)\,\,\,\,\,\nonumber\\ {\sigma _3} > 0 \Rightarrow \,\,\,{d_{{\text{mt3}}}}& = \sqrt {{{\left( {\frac{{{\sigma _3}}}{{{Z_{\text{T}}}}}} \right)}^2} + {{\left( {\frac{{{\tau _{13}}}}{{{S_{13}}}}} \right)}^2} + {{\left( {\frac{{{\tau _{23}}}}{{{S_{23}}}}} \right)}^2}} \,\,(b).\end{align} \tag{ 4 }$$

Fiber–matrix shearing failure:

$$\begin{align}{\sigma _1} < 0 \Rightarrow \,\,\,{d_s} = \sqrt {{{\left( {\frac{{{\sigma _1}}}{{{X_{\text{C}}}}}} \right)}^2} + {{\left( {\frac{{{\tau _{12}}}}{{{S_{12}}}}} \right)}^2} + {{\left( {\frac{{{\tau _{12}}}}{{{S_{12}}}}} \right)}^2}} .\end{align} \tag{ 5 }$$

Fiber failure:

$$\begin{align}{\sigma _1}\langle 0 \Rightarrow \,\,\,{d_{{\text{fc}}}} = \left| {\frac{{{\sigma _1}}}{{{X_{\text{C}}}}}} \right|\,\,\,\,(a)\,\,\,\,\,\,\,{\sigma _1}\rangle 0 \Rightarrow \,\,\,{d_{{\text{ft}}}} = \sqrt {{{\left( {\frac{{{\sigma _1}}}{{{X_{\text{T}}}}}} \right)}^2} + {{\left( {\frac{{{\tau _{12}}}}{{{S_{12}}}}} \right)}^2} + {{\left( {\frac{{{\tau _{13}}}}{{{S_{13}}}}} \right)}^2}} \,\,\,\,(b).\end{align} \tag{ 6 }$$

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 $\left( {{S_{12}},{\text{ }}{S_{13}},{\text{ }}{S_{23}}} \right)$, 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:

$$\begin{align}{\left[ S \right]_{ijkl}} = \left[ {\begin{array}{*{20}{c}} {\frac{1}{{{E_1}{d_1}}}}&{ - \frac{{{\upsilon _{12}}}}{{{E_1}}}}&{ - \frac{{{\upsilon _{13}}}}{{{E_1}}}}&0&0&0 \\[3pt] {}&{\frac{1}{{{E_2}{d_2}}}}&{ - \frac{{{\upsilon _{23}}}}{{{E_2}}}}&0&0&0 \\[3pt] {}&{}&{\frac{1}{{{E_3}{d_3}}}}&0&0&0 \\[3pt] {}&{}&{}&{\frac{1}{{{G_{12}}{d_4}}}}&0&0 \\[3pt] {}&{{\text{symmetric}}}&{}&{}&{\frac{1}{{{G_{23}}{d_5}}}}&0 \\[3pt] {}&{}&{}&{}&{}&{\frac{1}{{{G_{31}}{d_6}}}} \end{array}} \right]\end{align} \tag{ 7 }$$

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.

Zoom In Zoom Out Reset image size

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{\text{ 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{\text{ mm}}$ in diameter.

Zoom In Zoom Out Reset image size

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{\text{ 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 $\upsilon$ 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]

$$\begin{align}\sigma = {\sigma _0} + Q\left( {1 - {{\text{e}}^{ - b{\varepsilon _{\text{p}}}}}} \right)\end{align} \tag{ 8 }$$

where: σ₀ 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, $\mathop {\bar \varepsilon }\limits^. = 0.0012\;{{\text{s}}^{ - 1}}$.

Zoom In Zoom Out Reset image size

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{\text{ 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)

$$\begin{align}\int \frac{{{\sigma _{\max }}}}{\overline {\sigma }}{\dot {\bar \varepsilon}} {\text{d}}t \geqslant C \end{align} \tag{ 9 }$$

where, ${\sigma _{{\text{max}}}}$ is the maximum principal stress, $\bar \sigma$ is the effective von Mises stress and $\mathop {\bar \varepsilon }\limits^.$ is the effective plastic strain rate. C is the material constant threshold for damage.

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 C1-FEA-Present in the legend of figure 11(a) are the numerical micro strain value read from the locations depicted in figure 3, and represent the experimental [11], the base (coarse) mesh [11], improved (fine mesh) [11], and the present model results for C1 $\left( {0{\text{ }}\mu {\text{m}}} \right)$ clearance, respectively. The C4-MI and C4-FEA-Present in figure 3 indicate the results for C4 $\left( {240{\text{ }}\mu {\text{m}}} \right)$ 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.

Zoom In Zoom Out Reset image size

2.4.2. 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 $\alpha$ are derived by using equation (1), are detailed in the previous section.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

2.4.3. 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,{\text{ }}0.8,{\text{ }}0.86,{\text{ }}$ 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 ${\text{displacement }}x = 25{\text{ 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].

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

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].

Zoom In Zoom Out Reset image size

Damage parameters of composite and stainless steel materials obtained above will be monitored as an indicator of damage in future studies.