Bearing Failure Optimization of Composite Double-lap Bolted Joints Based on a Three-step Strategy Marked by Feasible Region Reduction and Model Decoupling

To minimize the mass and increase the bearing failure load of composite doublelap bolted joints, a three-step optimization strategy including feasible region reduction, optimization model decoupling and optimization was presented. In feasible region reduction, the dimensions of the feasible design region were reduced by selecting dominant design variables from numerous multilevel parameters by sensitivity analyses, and the feasible regions of variables were reduced by influence mechanism analyses. In model decoupling, the optimization model with a large number of variables was divided into various sub-models with fewer variables by variance analysis. In the third step, the optimization sub-models were solved one by one using a genetic algorithm, and the modified characteristic curve method was adopted as the failure prediction method. Based on the proposed optimization method, optimization of a double-lap single-bolt joint was performed using the ANSYS code. The results show that the bearing failure load increased by 13.5% and that the mass decreased by 8.7% compared with those of the initial design of the joint, which validated the effectiveness of the three-step optimization strategy.


Introduction
Advanced composites play important roles in aircraft structures because of their high strength/stiffness-to-weight ratios and the weight-saving requirements in the aerospace industry [Kumar and Srinivas (2018) ;Xu, Yang, Zeng et al. (2016) ;Zhang, Li, Wang et al. (2015); Sane, Padole and Uddanwadiker (2018)]. However, joints are commonly regarded as weak parts of composite structures because of the significant concentration of stress around bolt holes [Zhao, Yang, Cao et al. (2018); Zhao, Qin, Zhang et al. (2015); Liu, Lu, Zhao et al. (2018)]. In composite bolted joint design, increasing the load carrying capacity and minimizing joint mass have become crucial but challenging issues in the efficient use of composites. The load carrying capacity and joint mass are influenced by numerous multilevel parameters. Therefore, it is necessary to develop an efficient design approach to seek the best performance of composite bolted joints in engineering. To achieve a better structural performance and lower structural mass, various investigations have focused on the optimization of composite structures [Wu and Burgueño (2006); Lindgaard and Lund (2010); He and Aref (2003); Pelletier and Vel (2006); Manh and Lee (2014); Coelho, Guedes and Rodrigues (2015); Bakar, Kramer, Bordas et al. (2013); Kradinov, Madenci and Ambur (2007); Li, Huong, Crosky et al. (2009)]. Wu et al. [Wu and Burgueño (2006)] proposed an integrated approach for finding the optimal free-form three-dimensional shape and laminate stacking sequence design of fiber reinforced polymer shells. Lindgaard et al. [Lindgaard and Lund (2010)] proposed an approach for optimizing the nonlinear buckling fiber angle of laminated composite shell structures. He et al. [He and Aref (2003)] used a genetic algorithm to find the optimum design parameters of a fiber reinforced composite bridge deck, which included the number of stiffeners, the thickness and the orientations of outer skin layers. Pelletier et al. [Pelletier and Vel (2006)] presented a multi-objective optimization method for composite laminates based on an integer-coded genetic algorithm. The fiber orientations and volume fractions of the lamina were chosen as the primary optimization variables. Manh et al. [Manh and Lee (2014)] proposed an optimization model for aligning fibers in imperfect laminates based on a genetic algorithm and NURBS-based finite element isogeometric analysis. In addition to macroscopic parameters, such as the geometry configuration and stacking sequence, the mechanical performance of composite structures is also influenced by microscopic parameters. Coelho et al. [Coelho, Guedes and Rodrigues (2015)] used a multiscale topology optimization model to optimize bimaterial composite laminates in which the fiber cross-section was designed to obtain the optimal microstructure. Bakar et al. [Bakar, Kramer, Bordas et al. (2013)] presented a genetic algorithm to optimize the elastic properties of woven fabric composites. The design variables included the gap length, yarn thickness, material constituents and effect of the shape factor. Compared with the above simple laminates, only a few investigations have focused on the optimization of composite bolted joints. Based on a genetic algorithm and stress analyses, Kradinov et al. [Kradinov, Madenci and Ambur (2007)] achieved an optimum design of bolted composite joints. They considered the laminate thickness, lay-up, bolt locations, bolt flexibility and bolt sizes as design variables. Li et al. [Li, Huong, Crosky et al. (2009)] improved the bearing performance of composite bolted joints with reinforced z-pins. The volume contents and sizes of fibrous z-pins around the bolt holes were adopted as design variables. Clearly, recent investigations have generally focused on some of the design parameters of joints at the single scale, which has apparently restricted further exploration of the potential of composites. The performance of composite bolted joints is influenced by numerous parameters [Saleeb, Wilt, Al-Zoubi et al. (2003); Friedrich, Wu, Thostenson et al. (2011);Zhang, Zhou, Chen et al. (2016)], including dimension parameters, assembling parameters, unidirectional lamina, and the fiber and matrix, and an ideal optimization model involving all the relevant parameters for composite bolted joints is needed. In addition, an optimization strategy with high calculation efficiency is essential in engineering. Therefore, in this paper, an ideal optimization model for composite bolted joints that includes all the mechanical and geometric parameters is presented. Then, a three-step optimization strategy, including feasible region reduction, optimization model decoupling, and optimization based on the modified characteristic curve method, is presented. Finally, optimization of a single-bolt, double-lap joint is conducted to validate the effectiveness of the proposed ideal optimization model and three-step optimization strategy.

Ideal optimization model for composite bolted joints
A typical composite single-bolt double-lap joint is shown in Fig. 1. This figure shows a general multi-level framework of composite structures: from the microscopic fiber/matrix to the mesoscopic unidirectional lamina, macroscopic laminates and structure. An optimization model of composite bolted joints involving all the relevant parameters is an ideal model, which is the precondition for obtaining the optimal design of joints. The three levels of parameters that affect the failure performance and mass of composite bolted joints were teased out and are listed in Tab. 1. The macroscopic geometry variables include dimension parameters and assembling parameters, which can directly impact the ultimate failure load, failure mode and structural mass. The majority of the mesoscopic parameters are unidirectional lamina properties, which affect the structural mass and failure load. The microscopic parameters are fiber/matrix constituent properties, which have a small effect on the structural mass but a more significant effect on the failure load. In Tab. 1, the subscripts 1, 2 and 3 are defined with the material principal coordinate systems, as shown in Fig. 1. X and Y represent the longitudinal and transverse strength properties of the unidirectional lamina, respectively, and S refers to the shear strength of the materials. The subscripts T and C represent the tensile and compressive strengths, respectively. At the microscale, the fiber and matrix are regarded as isotropic materials. For the composite joint, D, e, Sw and t are decisive variables for the mass, and the mass can be mathematically expressed as follows: where ρ is the average density of the composite material. Obviously, m increases with the increase of the macroscopic dimension factors e, Sw and t of the laminate. The objective of optimization in the present work is to simultaneously maximize the bearing failure load and minimize the joint mass. Adopting the method of division, the ideal optimization model of composite single-bolt double-lap joints can be expressed as: S u b je c t to : b e a rin g fa ilu re m o d e where f is the mass-load ratio, m is the mass, and F is the failure load. The shear-out and tension failure modes, among other modes, are catastrophic failure modes and are undesirable in composite bolted joints design, while the bearing failure mode is the expected failure mode, and therefore, the optimization constraint is the bearing failure mode.

Three-step optimization of composite bolted joints
The three-step optimization strategy is presented. Then, the details of the strategy, including feasible region reduction, model decoupling, and one by one optimization with the modified characteristic curve method, are introduced.

Three-step optimization strategy
The flow chart of the three-step optimization strategy is shown in Fig. 2. For the ideal optimization model of composite bolted joints, the optimization is divided into three steps: feasible region reduction, model decoupling and one by one optimization. First, in feasible region reduction, the modified characteristic curve method is applied as the failure prediction method and the finite element model is established. Then, the dominant design variables are selected by sensitivity analysis, which reduces the dimensions of the feasible design region; influence mechanism analysis is conducted to reduce the feasible region of the dominant variables to further reduce the feasible design region. Finally, the simplified optimization model for composite bolted joints is proposed. By model decoupling, the simplified optimization model, which is multivariable, is further divided into several sub-models with relatively few variables by experimental design, failure analysis and variance analysis. In the optimization, based on the initial design variables and modified characteristic curve method, the sub-models can be solved one by one with the genetic algorithm to obtain the global optimal joint.

Ideal optimization model
Failure load prediction method Sensitivity analysis

Simplified optimization model
Optimal joints Step Step 3 The characteristics of the three-step optimization strategy include (1) sensitivity analyses coupled with influence mechanism analyses are conducted to reduce the feasible design region. Reduction of the feasible design region and computational complexity is achieved in two ways; (2) the original optimization model is decoupled into several sub-models to decrease the amount of calculations, which is an innovative idea. This idea works because the parameters of the joints are multi-level and easy to decouple; (3) the modified characteristic curve method, which is accurate and computationally efficient, is applied. The details of the optimization process of the ideal optimization model for composite bolted joints will be introduced in the following sections.

Modified characteristic curve method for failure load prediction
The large computation cost restricts the application of the progressive damage method (PDM) [Qin, Zhao, Xu et al. (2019)] for the optimization of composite structures, which is the most commonly used failure load prediction method for complex composite structures. Meanwhile, the characteristic curve method is suitable for predicting the failure load of composite bolted joints and attracts attention in engineering because of its simplicity and low computation cost as well as its good ability to accurately predict the failure load and failure mode. In the current work, a modified characteristic curve method constructed by Zhang et al. [Zhang, Liu, Zhao et al. (2014)] was adopted to predict the failure loads and failure modes of composite bolted joints for subsequent sensitivity analysis, variables decoupling and optimization solving. Joint failure was detected when the stresses in any ply satisfies the modified Yamada-Sun failure criterion at any point on the characteristic curve. The modified characteristic curve, as shown in Fig. 3, is expressed as follows: and r is the distance from the center of the hole to any point on the modified characteristic curve; R is the diameter of the hole; Rt, Rc and Rs are the tensile, compressive and shear-out characteristic lengths, respectively; and θ is anticlockwise measured from the symmetry axis, which ranges in the interval of [-90°, 90°]. Tensile failure occurs when the failure point is located at 75°≤  ≤90°; bearing failure emerges at the region of 0°≤  ≤15°, and shear-out failure occurs at 30°≤  ≤60°. The composite material system and joint configuration used by Zhang et al. [Zhang, Liu, Zhao et al. (2014)] is the same as the composite bolt joint in the current work. Therefore, this method and the corresponding characteristic coefficients are applied to the failure prediction in the present work.

Reducing the feasible design region to simplify the optimization model 3.3.1 Selection of dominant variables based on sensitivity analysis
Multi-level parameters have different influences on the ultimate failure load of composite bolted joints. To reasonably choose the dominant design variables, sensitivity analysis [Camanho, Bowron and Matthews (1998) ;Li, Gu and Zhao (2017)] was applied to establish a ranking of the various parameters according to the influences of the parameters on the failure load. Because only four parameters (D, e, Sw and t) are related to joint mass, the sensitivity of other parameters to the joint mass m will give results of zero, and the sensitivity of other parameters to the load transfer efficiency f is only related to their sensitivity to the failure load. Therefore, the failure load F instead of the load transfer efficiency f is the objective of the sensitivity analysis, and D, e, Sw and t are considered to be dominant design variables. In the sensitivity analysis, other dominant design variables are selected from other multi-level parameters xi (i=1, 2, · · · , n), where n is the number of parameters. The sensitivity analysis can be expressed by the following formula, which represents the change of the failure load per unit of relative variation of the parameters.
where xi refers to the variation of parameter xi. Each xi is determined by the standard deviation σ of the factor xi, which is referenced to the maximum deviation of each parameter specified in the ASTM D5961 standard [ASTM D 5961 (2013)]. Considering the different units of the parameters, a relative non-dimensional increment form /  ii xx was used in the sensitivity analysis to make the object evaluation comparable and fair.

Feasible region reduction of dominant variables based on the influence mechanisms
The dimensions of the feasible design region decrease with the decrease in the number of optimization variables, and thus, the amount of optimization calculations decreases. Another way to reduce the feasible design region is to reduce the feasible region of the dominant variables, which also reduces the computational cost. This paper introduces influence mechanisms of variables that affect the failure behavior to further reduce the feasible design region. This analysis of influence mechanisms involves three aspects: fiber orientation, geometry parameters and tightening torque. Each analysis yields feasible regions of some dominant variables, and finally, a simplified optimization model with various dominant variables for composite bolted joints is obtained.

Model decoupling based on variance analysis
There is still a high amount of computation in the optimization of the composite bolted joint because the number of dominant variables is large. As is known, the computation cost of an optimization model with n variables is much higher than that of the sum of n optimization models with one variable. To reduce the computational complexity, variance analysis is performed to study the interaction of variables. If the interaction term is not significant, a multivariate optimization can be decomposed into several optimizations with fewer variables or a single variable. Specifically, when the number of variables is known as n, the interaction term of any two variables is regarded as a factor and the quadratic term of the variables is taken into consideration; then, the total number of factors is 2n+n(n-1)/2. The experiments of the 2n+n(n-1)/2 factors are arranged with an orthogonal table, and the level number of each factor is chosen according to the amount of calculations, such as 3. After the failure load of each test is obtained by using the failure prediction method presented in Section 3.2.1, variance analysis is carried out to determine whether the influence of the interaction term is significant. If the n-1 interaction terms between a variable and other variables have little influence, the variable can be separately optimized, assuming that there are j variables in this case; if the interaction terms between mi variables and the other n-mi variables are not crucial, an optimization model with mi variables can be separated, assuming the number of optimization models in this case is k and i=1...k. Apparently, n=j+m1+…+mk. Finally, the original optimization model is decoupled into j+k sub-models, and each sub-model has fewer variables than the original one.

Conducting optimization with the one by one method
The sub-models will be solved one by one after the optimization model is decoupled, and then, the optimization results converge to a global optimal result. When a sub-model is executed, the variables in the other optimizations remain constant, which substantially reduces the computation time and simplifies the optimization procedure. The genetic algorithm [Manh and Lee (2014); Dong, Wei, Tian et al. (2017); Dong, Wei and Liu (2018)] is adopted for the optimization of sub-models. The optimization schedule is illustrated in Fig. 4. Step 1. Based on the initial values of the design variables (or the update values of the design variables), a genetic algorithm is used to obtain the candidate value of the variable being optimized, and the finite element model of a composite bolted joint is built to obtain an accurate stress field. Then, a modified characteristic curve method is used to predict the failure load and failure mode of the composite bolted joint.
Step 2. If the failure mode is bearing, but the sub-model has not yet been optimized according to the genetic algorithm, a new candidate value will be obtained with the genetic algorithm and Step 1 is conducted again; if the failure mode is not bearing, the candidate value is abandoned, a new candidate value will be obtained with the genetic algorithm, and Step 1 is conducted again; if this sub-model has been optimized by the genetic algorithm, the variables are updated and the next sub-model with other variables will be optimized until all the sub-models are optimal.

Optimization result
The three-step optimization strategy is applied to an ideal optimization model of a composite bolted joint, including feasible region reduction, model decoupling and one by one optimizing for the optimization of the sub-models.

Initial design of the single-bolt double-lap composite joint
The initial design of the single-bolt double-lap composite joint according to the ASTM D5961 standard and the investigations described in Zhang et al. [Zhang, Zhou, Chen et al. (2016); Zhang, Liu, Zhao et al. (2014] is listed in Tab. 2, and the laminates in the joint are made of T800 carbon/epoxy.  [Zhao, Shan, Liu et al. (2017)] was applied as shown in Fig. 5. All the components were modeled using the ANSYS ® SOLID185 element. The finite element model provided stress around the hole for the modified characteristic curve method, and the latter could predict failure load and failure mode with the stress. The computational mass m of the joint is 78.4 g, and the failure load F predicted by the modified characteristic curve method is 19.08 kN, which is in good agreement with the experimental result [Camanho, Bowron and Matthews (1998)] of 20.37 kN with a relative error of -6.3%. b) Selection of dominant variables based on the sensitivity analysis For the single-bolt, double-lap composite bolted joint with the initial parameters presented in Tab. 2, Tab. 3 lists the actual increment xi and the results of the sensitivity analysis Si (Eq. (5)). Fig. 6 also illustrates the sensitivity analysis results. It follows that the Si for the variables θ8, θ5, θ4, v23, t, µ, v12, θ1, θ7, θ9 and θ10 are less than 0.005 kN/1. They are trivial variables for the joint failure load, while other variables are of great significance, especially the three variables Vf, Ef and XC, which have the largest Si and are all larger than 1 kN.

Figure 6: Sensitivity analysis results
It can be proven that the majority of influential variables are material parameters, such as Vf, Ef, XC, E11, YC, Em, Gm, vm, tply, vf, Gf, S12, G12, S23, G23, E22, XT, YT, and η. The optimization design of the composite material parameters is performed by selecting excellent material systems rather than investigating a new material system. Optimizing the material parameters described above may lead to non-existent materials; therefore, this optimization is not regarded as the optimization content of this paper. The dominant design variables ultimately include the fiber orientation θ2, θ3 and θ6; geometry parameters e, Sw and D; and tightening torque N.

Reducing the feasible region of dominant variables
The influence mechanisms of the three kinds of dominant variables were studied to reduce the feasible region and thus reduce the amount of optimization calculations. a) Fiber orientation The various fiber orientations and the significant stress concentration around bolt hole cause complex stress status and further lead to different failure modes. The commonly used fiber orientations in the engineering industry are 0°, 90° and ±45°. The relationships between different fiber orientations and the failure modes of composite bolted joints are discussed.
If θ is 0°, as shown in Fig. 7, the fiber direction is parallel to the loading direction. Because of the poor matrix and fiber-matrix interface properties, the initial damage at point D is easy to propagate along the shear-out plane and leads to a sudden shear-out failure. Akatasa et al. [Akatasa and Dirikolu (2004)] and Park [Park (2001)] proved that excessive 0° layers in a composite bolted joint caused unexpected failure load reduction and shear-out failure mode. If θ is 90°, as shown in Fig. 8, the stress at point B is in the same direction as the fibers, so that the fibers inhibit damage initiation at point B. Because the matrix properties and fiber-matrix interface properties are too small to sustain corresponding stresses in the direction perpendicular to the fiber, interface debonding or matrix cracking occurs at both points C and D, which leads to tension failure. If θ is 45° or -45°, the angles between fiber direction and stress direction at points B, C and D are all 45°, as shown in Fig. 9. Both compressive and tensile stress can be decomposed into the longitudinal and transverse directions of fibers, which means the fibers bear a significant degree of load at points B, C and D. Thereby, the possibility of the occurrence and propagation of damage can be reduced compared with that at the 0° and 90° layers. The selection of the stacking sequence is important in composite bolted joint design. The reasonable proportion of the 0° layers contributes to increasing the failure load of composite bolted joints [Quinn and Matthews (1987)]. The ±45° layers can reduce the stress concentration factor and avoid tension and shear-out failure mode. The 90° layers are able to prevent shear-out failure and maintain the transverse loading capacity. The other suggestion is that the 0° and 90° layers should not be placed at the surface, which has the largest bearing stress at point C. In addition, adjacent layers should have different fiber orientations to obtain a better interlaminar performance [Quinn and Matthews (1987); Collings (1982)]. Therefore, the constraints on the fiber orientations are θi≠θi+1 and 10%≤ n45/-45/0/90 ≤50%. b) Geometry parameters The failure behavior and joint mass are influenced by the geometric dimensions. Fig. 10 illustrates the typical stress status of the central laminate, which shows the effects of the geometric dimensions on the failure behavior of a composite joint. Those stress statuses are related to four common laminate failure modes in a composite bolted joint, namely, bearing, tensile, shear-out and cleavage [Maimi, Camanho, Mayugo et al. (2007); Camanho and Matthews (1997)]. Obviously, Section A experiences a distributed force, and a high stress concentration is exhibited at points D and D'. For the small edge distance Sw, a high stress distribution appears at section A, which directly results in a sudden tension failure, as shown in Fig.  10(b). A sufficient edge distance Sw ensures the stress at section A is low enough to avoid catastrophic tension failure. In addition, a short end distance e deteriorates the stress states at points B and D, which leads to cleavage or shear-out failure, as shown in Fig.  10(c). Therefore, the end distance e should be long enough to prevent the laminate from experiencing cleavage and shear-out failure. It is worth noting that tension, shear-out and cleavage failure always lead to catastrophic structure failure [Chang and Chang (1987); Pisano, Fuschi and Domenico (2013)]. However, the bearing damage slowly propagates near point C (Fig. 10(d)) with the increase of the external load, which provides a failure alert [Chang and Chang (1987); Guo and Nairn (2017)]. Hence, bearing failure is the expected failure mode for a composite bolted joint. Appropriate e, Sw and D not only contribute to increasing the failure load but also to avoiding catastrophic failure. Recent investigations have suggested that the ratios of e/D and Sw/D should be greater than 2.5 to ensure the bearing failure mode [Zhang, Liu, Zhao et al. (2015); Hart-Smith (1976)]. In addition, the upper limit of e and Sw should be restricted to decrease the composite joint mass. In the current work, e/D and Sw/D, denoted by  e and  S , respectively, were studied instead of e and Sw, and the optimized ranges of  e and  S were set as 2.5≤ e ≤4 and 2.5≤ S ≤4.
Under an external load, the bending deformation of a bolt leads to the reduction of the contact areas between the bolt and laminates. Then, a significant stress concentration occurs at the reduced contact areas, which deteriorates the damage status around the bolt hole and decreases the load transfer efficiency. Moreover, the thicker the central laminate, the smaller the bending deformation of the bolt, which results in a smaller stress concentration at the contact area, as shown in Fig. 11. To obtain the highest bearing load transfer efficiency, the recommended thickness of the laminate is 1≤D/t≤2 [Collings (1977)].  The failure load of a composite bolted joint is affected by the tightening torque N. The tightening torque causes interlaminar compression pressure around the bolt hole, which can inhibit the interlamination damage and bearing damage around the bolt hole, as shown in Fig. 12 [Sun, Chang and Qing (2002)]. The tightening torque also increases the friction force in a composite bolted joint, which contributes to the failure load improvement and prevents the adjacent layers from shear slipping [Phillips (1989); Khashaba, Sallam and Shorbagy (2006)]. Furthermore, the tightening torque provides an axial preload to the bolt, which can reduce the bending deformation of the bolt and relieve the bearing stress concentration at the outer layer of the laminate, as shown in Fig.  12(d). Thus, the tightening torque N should be included in the optimization, and the constraint is 0≤N≤10.
" "is stress caused by load

Simplified optimization model for composite bolted joints
Because of the bolt bending deformation, the bolt-hole bearing stress in the outer layers is much higher than that in the inner layers. Therefore, the failure behavior of the joint is more affected by the outer layers, and the fiber orientations of the outer layers are decisive parameters of the failure behavior in the optimization design process. Thus, this paper only takes the fiber orientations of the outer five layers as optimization parameters. The final design variables in the present work are D,  e ,  S , N, and θi (i=1, 2, 3, 4, 5). The constraint conditions of the optimization variables are based on the discussion in Section 4.2.2. Consequently, the simplified optimization model for composite bolted joints can be expressed as: where n is the total number of plies, n±45/0/90 is the number of plies in different fiber orientations, and θ is the fiber orientation.

Variable decoupling
The fiber orientation θi can be decoupled from the other four variables D,  e ,  S , and N because they are parameters in different levels. These four variables can be transformed in the following way: A=d/6, B=  e , C=  S , and D=N/5. With the addition of six interaction terms, AB, AC, AD, BC, BD, and CD, and four quadratic terms, A 2 , B 2 , C 2 , and D 2 , of the four variables, there are fourteen factors in total. Five levels within the constraints are set for each factor for the orthogonal experimental design. Then, the failure prediction of models under different orthogonal experiment conditions is carried out by using the method presented in Section 3.2, and variance analysis is conducted on the failure prediction results. The P value of the significance test for each factor is shown in Tab. 4.

Optimization based on the modified characteristic curve method
Optimization of the five optimization sub-models for the design of the composite bolted joint is performed according to the procedure presented in Section 3.5. The optimization history is shown in Fig. 13, where the objective function values are plotted for each optimization iteration. There is only one step in optimizing the optimization sub-model (1) because the initial design value of 4.76 mm, as shown in Tab. 2, is the best choice from the available standard bolt diameters. The other four sub-models require more steps for analysis, among which the optimization sub-model (5) has the longest iteration because it includes five variables, although each variable only has four candidate values. The minimum objective function value of 3.307 in the bearing failure mode is noted. The optimization design results are shown in Tab. 5 and are compared with those of the initial design. It can be seen from Tab. 5 that the bearing failure load of the optimized joint is increased by 13.5% and the mass is reduced by 8.7%. Therefore, the proposed optimization method for composite bolted joints is meaningful and effective.

Conclusions
This paper presents a three-step optimization strategy based on feasible region reduction, model decoupling and the modified characteristic curve method for composite single-bolt double-lap joints. The three-step optimization strategy mainly has the following three characteristics: (1) the feasible region is reduced from two aspects to reduce the computational complexity, (2) the optimization model is decoupled into several submodels that are easier to solve to reduce the computation cost, and (3) the modified characteristic curve method is applied for failure prediction in the optimization. An ideal optimization model for composite bolted joints is established and simplified into the optimization model with dominant design variables; then, the simplified optimization model is decoupled into 5 sub-models and solved. Comparing the result of optimization design with the result of the initial design joint, which is based on the recommendations in the ASTM D5961 standard, the final failure load of the optimized composite joint is increased by 13.5% and the mass is reduced by 8.7%. This comparison illustrates that the proposed optimization strategy is an effective method for minimizing the mass and increasing the bearing failure load of composite bolted joints.