## 1. Introduction

Composite materials have many structural advantages compared with general metals such as iron and aluminum. They are used in many areas because they are lightweight, have high specific stiffness and specific strength, have small deformations due to temperature, and are easy to develop into materials with various properties according to the application. In particular, many composite materials are applied to the aerospace, medicine, national defense, and advanced industries.

To produce comprehensive composite structures comprised of various parts, the composite parts must be combined in various shapes. The methods for combining composite materials can be largely divided into adhesion and fastening. Currently, the combination method of fastening is used most frequently because producing structures by adhesion has limitations. For fastening, holes must be created in composite structures, and the stress that causes these fastening holes has been researched as an important topic^{[1-4]}. Stress concentration is a very important design variable in the design of holes. Existing studies have suggested various methods to minimize stress concentration^{[4,5]}. However, they were mostly limited to isotropic materials.

This study proposes a method of determining the optimal positions of fastening holes in a laminated composite structure. The hole positions were determined on a two-dimensional plane, and the x and y coordinates were determined under various conditions. The hole positions were optimized with an optimization algorithm by using x and y coordinates on the plane as design variables. For optimization, the response surface method was used, which has been frequently applied recently due to its effectiveness^{[6]}. The stress that occurs around the holes can be minimized by optimizing the position of each bolt hole in a structure supported by multiple bolts, which improves the reliability and safety factor of the design.

## 2. Theoretical Background of Numerical Analysis

### 2.1 Formulation of finite element analysis and numerical analysis

The composite structure analyzed in this study is wooden laminated veneer lumber (LVL), which is made by laminating multiple layers of thin wooden sheets with wood grains. Wood is an anisotropic material with different degrees of stiffness and strength according to the direction. Therefore, a 3D composite finite element equilibrium equation must be employed to analyze this structure with the computational finite element method. The core part is determining the 3D anisotropic stiffness matrix.

The composition of the material stiffness matrix and the inducement process of the equilibrium equation are as follows. The equation is derived with the principle of virtual work, which derives an equation from the equilibrium of internal and external forces. The principle of virtual work considering the equilibrium of forces at time *t* and *t* + *Δt* is as follows.

The left term in Eq. (1) represents the internal force, and the right term represents the external force. Here, ${}^{t+}{S}_{ij}\Delta tlsub$ 0is the second Piola– Kirchhoff stress tensor and ${}_{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0}^{t+\Delta t}{\in}_{ij}$ is the Green–Lagrange strain tensor. When Eq. (1) is divided into linear and nonlinear strain components and expressed as an incremental equation, Eq. (2) is obtained.

Eq. (2) consists of a stiffness matrix of linear and nonlinear components, external force, and internal force. When it is expressed as a matrix, the following equation is obtained:

where *K _{L}* is the linear stiffness matrix and

*K*is the nonlinear stiffness matrix. The linear stiffness matrix can be expressed as Eq. (4).

_{NL}

The nonlinear stiffness matrix can be expressed as Eq. (5).

The internal force vector can be expressed as Eq. (6).

The material stiffness matrix [*D*] is as follows:

Each component of this matrix is expressed as follows:

where *m* = cos*θ* and, *n* = sin*θ* . When the coordinate that determines the directivity of the material coincides with the definition of the global coordinate system, *Q* is defined as follows:

where [*Q*] is a flexibility matrix.

### 2.2 Optimization by response surface method

The optimization is an algorithm that finds the maximum or minimum point in the region of the design variables. Traditionally, the method of finding the maximum and minimum points by the gradient method has been used frequently. However, many methods that have complemented the shortcomings of this gradient method have been proposed. The representative methods are the probabilistic model-based genetic algorithm and the response surface method. The response surface method first calculates the reaction surface by multiple calculations in the design function area of the variables. The optimization precision is determined by how accurately this reaction surface is implemented. To verify the accuracy of the maximum and minimum point values calculated on the reaction surface, the objective function of the actual optimization must be calculated by the combination of final optimized design variables. There are many methods for calculating the reaction surface from multiple design functions, and this paper uses the Kriging method. This method defines the sum *y*(*x*) of the polynomial function *f*(*x*) and the normal distribution function *z* (*x*) .

The covariance matrix *z* (*x*) , an interpolation model that uses multiple sample points, can be expressed as follows:

where *R* is the correlation matrix, *r*(*x ^{i}*,

*x*) is the spatial correlation function, which is expressed as follows:

^{j}

where *θ _{k}* is a variable and

*M*is the number of design variables. The large advantage of the response surface method is the reduction of time and cost because once the reaction surface is completed, the maximum and minimum are found on this surface.

## 3. Optimization Analysis

### 3.1 Analysis model and definition

Optimization was performed for a wooden laminated composite beam to determine the optimal positions of bolt holes, which are fastening elements. The size of the beam, which is the analysis model, was 4,270mm long, 406.4mm wide, and 129.5mm thick. This laminated composite beam had the boundary condition that it was fastened under the condition in Fig. 1. The diameter of each bolt hole was 19.05mm, and the bolts were fastened in four bolt holes. In the actual finite element analysis, a 1/4 model (1/2 in the thickness direction and 1/2 in the length direction) was applied due to the symmetry of the model. The beam was composed by laminating thin sheets in 12 layers. The shape, dimensions, and layout are shown in Fig. 2. The thin wooden sheets of each layer were an anisotropic material with directivity.

The beam received a distribution load of 2.7E-3MPa. At the bottom, it was loaded at four points, as shown in Fig. 3. The load was 10,666.5N at the center and 5,333.4N at the edge. This distribution of forces was calculated by the equivalent concentrated load assuming a constant pressure.

Fig. 3 shows the total analysis model. The distribution load acting on the top was the stress corresponding to the self-weight. The most important consideration for the optimization analysis for a beam that is subject to bending is minimizing the stress concentration that occurs around the bolt hole.

Table 1 shows the properties of each layer. *X _{t}* denotes the tensile strength in the

*x*direction,

*X*denotes the compressive strength in the

_{c}*x*direction, and

*S*denotes the shear strength of the

_{xy}*xy*plane.

Composite materials have different degrees of stiffness and strength by direction. Thus, various theories and methods can be applied to the definition of failure. This study examines the failure of composite materials during optimization by applying the Tsai–Wu failure theory, a representative theory that defines the failure of composite materials.

The failure of composite materials was evaluated by defining the Tsai–Wu failure index . Failure of the composite material was defined as occurring when this value reached 1.

To optimize the hole positions, the coordinates that determine the hole positions were set as design variables. Eight design variables and six constraint conditions were given for optimization conditions, including the coordinate variables *x*_{1}, *y*_{1}, *x*_{2}, *y*_{2}, *x*_{3}, *y*_{3}, *x*_{4}, *y*_{4} which determined each hole position and the constraint conditions that determined the interval of each hole. Many studies have been conducted on the optimization of anisotropic materials, such as composites and lumber^{[7-10]}.

The purpose of optimization is to minimize the von Mises stress that occurs around the bolt holes and ultimately find the optimal hole positions. This optimization algorithm, which finds such critical points, provides a method for solving complex engineering problems. The equation for performing this optimization analysis is as follows.

Object function

$F\left({x}_{1},\hspace{0.17em}{y}_{1},\hspace{0.17em}\hspace{0.17em}\xb7\hspace{0.17em}\hspace{0.17em}\xb7\hspace{0.17em}\hspace{0.17em}{x}_{4},\hspace{0.17em}{y}_{4}\right)={\left({\sigma}_{von\hspace{0.17em}Mises}\right)}_{\text{max}}$

Design variables

${x}_{1},\hspace{0.17em}{y}_{1},\hspace{0.17em}{x}_{2},\hspace{0.17em}{y}_{2},\hspace{0.17em}{x}_{3},\hspace{0.17em}{y}_{3},\hspace{0.17em}{x}_{4},\hspace{0.17em}{y}_{4}$

Constraints

$\begin{array}{l}1).\hspace{0.17em}76.2\le {x}_{2}-{x}_{1}\le 136.0\\ 2).\hspace{0.17em}323.8-{x}_{2}\ge 57.15\\ 3).\hspace{0.17em}76.2\le {x}_{4}-{x}_{3}\le 136.0\\ 4).\hspace{0.17em}323.8-{x}_{4}\ge 57.15\\ 5).\hspace{0.17em}114.3\le {y}_{1}+{y}_{3}\le 276.22\\ 6).\hspace{0.17em}114.3\le {y}_{2}+{y}_{4}\le 276.22\end{array}$

The constraint conditions were applied to determine the areas for hole positions.

### 3.2 Analysis results and review

The total model for computational structural finite element analysis was modeled with 3D Hexa elements. The stress calculation accuracy was improved by reducing the element size at the bolt hole position. The initial bolt hole positions are shown in Fig. 4. The specific positions are represented as the accurate coordinates of the holes in Table 2. For the finite element model in Fig. 5, the hole position that represents the minimum stress (von Mises stress) in the bolt hole under the boundary condition given in Fig. 3 was determined by the optimization method of the response surface method. The optimization results are outlined in Table 2. The maximum stress at the initial bolt hole position was approximately 23.2 MPa, and the stress value when it was optimized was 20.0 MPa. Thus, the maximum stress decreased by approximately 16%.

To calculate the stress around the bolt holes accurately, a boundary condition for calculating the stress that occurred when the bolt came into contact with the inner surface of the hole was given. The optimized hole positions and stress are shown in Figs. 6 and 7.

Fig. 6 shows a comparison between the initial and final bolt hole positions to facilitate the determination of the optimization results for the hole positions. The failure index presented in Eq. (13) was examined to verify the failure condition of the composite material at the optimized bolt hole positions. The failure index under this load condition was approximately 0.8. In future research, the most sensitive objective function could be found during the design by comparing the optimization results for a multi-objective function. The optimization of the multi-objective function requires calculations that are more complex. Fig. 7 shows the stress distribution of the bolt holes that represent the maximum stresses of the initial and optimized states. The maximum stress appeared at the bottom right in the initial state and at the top right in the optimized state.

## 4. Conclusions

The optimization of bolt hole positions for a laminated composite beam was calculated. The model, method, boundary condition, and stress calculation method for hole position optimization are presented. The optimization method reduced the maximum stress by approximately 16%. Bolt fastening is the most basic fastening method. The load support capacity of the total structure can be enhanced greatly by calculating the hole positions that are more reliable and have a higher safety factor than existing general hole positions with the optimization method proposed in this study.