Optimization of Wall Thickness and Lay-up Design for Composite Structure with Open-Shell Geometry

In recent years, many scientific and engineering groups have developed a variety of approaches and numerical methods for solving the problem of composite structural optimization. Most of these approaches and optimization tools are developed by the key specialists and scientific groups [1-3] and are presented in studies [4-6]. When optimizing composite structures with shell-like geometry as the objective function the total strain energy, the peak von Mises stress or the averaged strain are most often used, and they should be minimized at some design constraints (total weight, external or internal shape, etc.) [7]. The result of such optimization is also an increase of the structural vibrations eigenfrequencies.

The optimization of composite structures consists of the three basic phases [2, 3]. On the first, the free-size, or the topology optimization phase, a total laminate thickness and the optimal laminate composition can be changed 'continuously' throughout the structure [3]. On the next phase the thickness of each ply bundle varies in order to optimize all the objectives considering the manufacturing constraints. Third phase [2,3] represents a shuffling (redistribution) process to determine the optimal stacking sequence, providing the balance of laminate [5, 6]. This optimization procedure can be simplified for the composite parts, manufactured using quasi-isotropic laminates with the same Young's modulus Y and same Poisson's ratio $v$ in all in-plane directions, but not normal to the plane [5-7]. Such quasi-isotropic laminates can be produced when the total number of unidirectional layers is $n \geq 3$; all the layers have equal orthotropic elastic constants and equal thickness; the orientation of the $k^{\text {th }}$ layer of a n-layered laminate is taken in the form $\theta_k=\pi \cdot(k-1) / n$ [5]. Quasi-isotropic material of the optimized structure is considered as isotropic [8].

In cases of the real structure optimization a CAD model of the optimized structure is required. This model is further converted to the FE representation $[1-3,7,9-10]$. In order to determine the stress-strain state in the most stressed areas of the optimized part, a forward problem for the given boundary conditions and the operating load cases is formulated, and FE simulations are carried out, in some cases by using numerical CFD analysis [1, 3, 7].

Known optimization algorithms can be divided into two large groups. The first group includes methods for topology optimization [7, 11-13] (ESO, BESO, SIMP), which are inherent a great number of degrees of freedom (the design variables). The second group includes methods which reduce the dimension of the design variable space to the lowest possible number of parameters. In order to perform such a reduction, the varied geometry should be represented as a function that depends on a finite number of variable parameters, by a pointwise presentation of the material's distribution over the structure [14, 15], or by the gauge optimization models when the whole part is divided onto 10 and more distinct zones to be separately optimized [2, 3].

When the composite shell-like structure is optimized, it is required to determine the orientation of plies on a mould surface with bi-directional curvature by means of draping analysis [2, 6, 9]. A mapping approach to the problem of draping has been used in this work. This approach assumes a geometric mechanism to transformation of an initial unit square of fabric into the draped shape at such hypotheses [9]:

• Inextensibility of fibers;

• Continuity of displacements on the boundary between two adjacent layers;

• Perfect mapping of the composite plies onto the mould surface.

In spite of the existing CAD and highly targeted software [2, 9], simulation of the darts (cuts of roving or fabric) location is still a very complicated process. In the present work the problem of structural optimization of rotorcraft radar cowling with open-shell geometry, which undergoes the spatially distributed pressure from the airflow at the climbing / descent $(\pm 10 \mathrm{deg})$ modes of flight with a speed $300 \mathrm{~km} / \mathrm{h}$, was researched. The CAD model of the cowling has been converted into the FE representation of its external surface assumed to be a stiff shell according to the purpose of minimizing the total strain energy equivalent to the mean compliance of the structure. In order to find the most stressed mode of the flight an air blasting in a wind tunnel of the shell was simulated at the different orientation of the airflow. The spatial distribution of the air pressure acting on the shell surface was assumed as the solution of the CFD forward problem. On the next step the structural dynamics problem for this shell loaded by the previously found pressure field has been solved. Simple parameterization was proposed for this problem. According to this parameterization the varied radius and centre coordinates of an auxiliary sphere are used as the design variables.

An area where the shell walls should be stiffened by local thickening, defined as the smooth function, was bounded by the intersection of the shell and the auxiliary sphere. A direct search of an optimum design is performed inside a subset of 4D design space at constrained maximum wall thickness and total weight of the part, by means of the optimization algorithm, that calls the forward problem for the most stressed loading case. A thickness of the shell walls distribution, obtained as a result of the structural optimization, has been then used for the design of the composite prepreg draping. In order to obtain the draped shape of composite sheets lay-up the shell surface with two kinds of perpendicular geodesic lines was deployed on the plane.

Many aircraft and rotorcraft composite cowlings based on a glass-fiber epoxy represent a thin walled shell which is fastened on the stiffer part of the airframe. Such a cowling is subjected to a transient action of airflow that can invoke undesirable mechanical vibrations as well as deformations and noise. In order to determine an intensity and spatial distribution of air pressure on the shell surface, we have neglected the fluid-structure interaction, and assumed that the shell is a non-deformable rigid body, and its external surface has the shape of the studied shell (see Figure 1). All the dimensions were assumed to be small as compared to the dimensions of the airframe. Since a CAD model of the cowling was imported to Comsol Multiphysics software, all the side surfaces have been extruded to the distances much greater than real dimensions of the cowling.

The problem has been formulated by means of the Reynolds averaged Navier-Stokes (RANS) equations [16]. The $k-\omega$ turbulence model has been used due to its best behavior near the walls and more precise predictions of the free shear flows in comparison with the $k-\varepsilon$ model. The $k-\omega$ turbulence model with two additional equations relative to the turbulent kinetic energy $k$ and the dissipation per unit turbulent kinetic energy $\omega$, enclosed the RANS equations, which simulate the transport of a turbulent kinetic energy $k$ and $\omega$.

where the previous equations use the following constants and functions

where $\Omega_{i j}$ and $S_{i j}$ are the mean rate-of-rotation and rate-of-strain tensors, respectively, defined by

and the turbulent viscosity is

In order to determine the distribution of total pressure and velocity field along the cowling's wall, the transient statement of the turbulence flow problem has been formulated. The airspeed in inlet started from zero, slowly increased and reached its steady value of $\sim 83 \mathrm{~m} / \mathrm{s} (300 \mathrm{~km} / \mathrm{h})$, which is close to the helicopter's maximal speed. An outlet boundary condition was assumed to be $p=0$, when the conditions are represented in the form of wall functions

where $\mathbf{n}$ is the unit normal vector to a wall. These functions were used both on the walls of the wind tunnel and on the surface of the cowling that assumed to be absolutely stiff.

During a virtual test, the model has been placed into the semicircular tunnel's section at the different model orientation relative to airflow, which is aligned with the tunnel axis. The finite-element mesh of the air medium included about 90000 tetrahedron elements (see Figure 2).

The transient problem has been solved for the angles of model's orientation relative to the tunnel axis varying in the range $\pm 10 \mathrm{deg}$, with an increment of 2.5 deg. The spatial distributions of the pressure field, which acts on the model surface (see Figure 3), have been written in the text files for the subsequent digital processing.

The shell under study has been considered as isotropic during the determination of its stress/strain state. For the multilayered composite manufactured by laying-up of a fabric prepreg the following properties were accepted: the Young's modulus equals to $2 \cdot 10^{10} \mathrm{~Pa}$, Poisson ratio - 0.3 and density - $2,200 \mathrm{~kg} / \mathrm{m}^3$. The shell model of the cowling structure (see Figure 4) with homogenous wall thickness of 5 mm has been used in order to formulate the forward problem. The geometry of the cowling is converted from the CAD model. Before the optimization, its volume and weight were equal to $2.86 \cdot 10^{-4} \mathrm{~m}^3$ and 1.3 kg, respectively.

As the objective function the total elastic strain energy of the shell has been used. It was computed during the solving of the forward problem, when the surface of the shell has been loaded with the non-uniform quasi static pressure field calculated in previous study. Table 1 demonstrates the values of the energies. It is obvious that the airflow orientation value of -10 deg corresponds to the case of the maximal stresses.

Due to the formation of stillbirths structures and very low stability of such topology optimization methods as SIMP [12], ESO, BESO [13], they did not allow to optimize the shell's thickness distribution. It was decided first to parameterize the problem and thereby to reduce the dimensionality of the design space. The problem has been parameterized by a special artificial technique. According to this technique the thickness of the shell should be optimized in an area which is confined by the intersection of the fictive sphere with the shell surface (see Figure 5). In this area a smooth function $K\left(x_p, y_p, z_p, R\right)$ of 4 design variables describes the spatial distribution of the shell thickness.

where $h_0=5 \mathrm{~mm}$ is the initial thickness of the shell, $x_p, y_p, z_p$ and $R$ are the coordinates of center and radius of the auxiliary sphere, respectively.

The optimization process is realized by means of a cyclic call of the forward problem. During this process the domain in the space of design variables is uniformly covered with the points corresponding to different sets of parameters $x_p, y_p, z_p, R$. Then this domain is mapped onto the space of objective functionals. The optimization results have been saved into the text files, which contain the spatial thickness distributions for each set of design variables. All the designs that didn't satisfy the constraint on the cowling weight $W \leq 1.7 \mathrm{~kg}$ were not considered. For the angle of -10 deg the total number of simulated versions exceeds 500. In the Figure 6 the values of total strain energy combined with the values of the mass shell for the whole set of considered designs are presented. There is a one-to-one correspondence between every generated set of values of the design variables and the number of variants, so long as an each design variable takes on a consecutive value with an equal increment within its variation interval.

In Table 2, the best two solutions of the optimization problem are presented; they are combined with the values of total and maximum strain energies of a not optimized shell. These two solutions are very close to each other.

On this stage of investigation the coordinates of the optimized cowling's shell and the spatial distribution of the wall thickness have been numerically analyzed by means of MATLAB tools. At first the shell geometry and optimized wall thickness distribution have been represented by the wireframe parametric surfaces plotted from the set of discrete points.

A mapping approach has been used for the studied aircraft part with the bi-directional curvature. According to this approach the shell unfolding into a plane surface is required. The assumptions formulated in the introduction have been accepted for the draping analysis. A lay-up process starts with the two geodesic lines, which are perpendicular. The origin of a new coordinate frame was taken at the highest point of the dome. At this point the two afore-mentioned kinds of lines should be identified. All the new axes must be parallel to the global frame axes. In order to flatten the curves of the curved shell surface, and thereby to determine the dimensions of the pattern of prepreg, the two sets of scattered points have been approximated by the functions that interpolate the curved surface and wall thickness spatial distribution using triangulation algorithm. The used deployment algorithm includes the following steps. The first step consists in a cyclic section of the shell surface with the planes, passing through a vertical axis of new frame and turned around it with angular step $\Delta \phi=5$ deg. Further the length of the section line $z(x(r, \phi), y(r, \phi))$ is integrated according to the formula

where, the polar coordinates $r, \phi$ are determined in the area occupied by the shell. Because the shell is subjected to the manufacturing by open die molding and can be easily removed from the mould after the cure, this fact guarantees the dissimilarity of $z(r, \phi)$ along all the integration paths, which permits to carry out the transformation $(r, \phi) \rightarrow(l, \phi)$ and to express the function of wall thickness distribution through the new coordinates $h(r, \phi) \rightarrow \tilde{h}(l, \phi)$. The new function is defined on the plane (see Figure 7). Thus, the interpolation of the function $\tilde{h}(l, \phi)$ allows to construct the pattern of lay-up (see Figure 8).

In the Figure 8 the contour lines appear with a step of 0.5 mm thickness. The four layers of fabric oriented according to the scheme $\left[0^{\circ} ; 90^{\circ} ;+45^{\circ} ;-45^{\circ}\right]$ can be placed within the areas surrounded by each contour line. Such a plot allows to optimize a final stacking sequence by means of ply shuffling procedure without difficulty. In order to fix the patches of the fabric in the exact placement more than three knots must already be pinned onto the mould surface due to the complexity of a lay-up pattern. Figure 8 helps to make a conjecture about the correct positions of these knots. In this investigation, an excessively complex problem of the rational dart placement has not been studied due to the wrong and fruitless assumption about inextensible prepreg's fabric that it has to be used in known methods [2, 6, 9]. The present approach to the lay-up design has a particular application for the allround convex shells.

In this study an approach to the structural optimization of the aircraft convex-like shell structure is presented. This approach includes a preliminary CFD study of the structure at the operating airflow velocity and orientation relative to the shell to determine the most stressed case and pressure field, which acts to external surface of the shell. The forward structural mechanics problem for the shell loaded by the found airflow induced pressure has been numerically solved using the aforementioned results. For the finite element simulation, we assumed the shell to be quasi isotropic. The reinforcement of the shell's structure via the local thickening of its wall was proposed in order to optimize the global stiffness (and reduce a global elastic strain) of the shell. In order to define the set of design variables we consider the auxiliary sphere, whose intersection with the shell surface indicates the area where the shell's thickness should be increased. As the design variables we accept the radius of the sphere and the coordinates of its center. The growth of the wall thickness has been defined by the smooth function so as to eliminate the probable local stress concentration. During the optimization algorithm the finite element model with the varied design parameter is multiple called at the constrained total shell's weight. An optimal solution provides a two and half times less global stain energy, a reduction of $22 \%$ in the von Mizes stress peak, and $5 \%-15 \%$ growth of the vibration frequencies on the first six natural modes. Finally, the lay-up of prepreg has been constructed by means of the mapping method in respect to unfolding of the fiber-reinforced composite sheets placed on a curved tool surface.