Energy absorption assessment of conical composite structures subjected to quasi-static loading through optimization based method

The geometry of a composite frusta-conical structure with 20° apex angle was optimized using topology optimization method. The ATOM topology solver was used to perform optimization process. The structure was subjected to an axial compressive quasi-static loading where the criterion of optimization is the largest absorption of energy per unit volume. The optimization resulted in concave shell geometry with a volume reduction of 15%. To compare the effect of geometry on energy absorption capability, two types of GFRP specimen with different thickness were fabricated using filament winding and hand lay-up methods. The axial crushing test results were validated with numerical analysis using the commercial finite element software ABAQUS. The collapse mode and failure mechanism were investigated and identified. The presented methodology helps to find the better geometry instead of testing different geometries.


Introduction
Laminate composite materials have been introduced in automotive vehicles in order to decrease weight but with added benefits of high specific strength, specific density and modulus. Recent use of laminate composites has increased due to high crashworthiness capabilities [1]. Evaluation of different types of composite materials and structural geometry is necessary and this has resulted in many experimental tests using different conditions of loading [2][3][4][5]. These studies were focused on the effective parameters involving crashworthiness indicators such as specific energy absorption, mean crushing load and peak load [6]. Utilizing different geometries [4], thickness, fiber orientation angle [7,8] as well as type of material, the effective parameters have been examined for the best crashworthiness efficiency.
Various methods of fabrication such as filament winding, hand lay-up and vacuum infusion were introduced to study the effect of fabrication methods on the energy absorption capabilities [9]. Different impact tests were conducted [9][10][11] to assess the role of material properties and production method on failure and fracture modes. Conical thin-walled composite structure is one of the most studied geometries in crashworthiness analysis. Conical shapes with different chord angles, thickness, number of laminates and aspect ratios were investigated [12,13].
Different analytical approaches have been carried out to examine crashworthiness. Experimental studies were done to evaluate the effect of apex angle in the range of 5°-25°. Previous studies [14] showed that increasing structure thickness induces better progressive collapse in glass/epoxy composites. Mamalis et al. [15] have performed a numerical study on axial crushing for thin walled fiberglass composites with different apex angles of 5°-20°. They examined the effect of aspect ratio on measured peak load. Esnaola et al. [3] conducted a study on the effect of different shapes but with the same aspect ratio of t/d (wall thickness and outer radius). He reported that tubes have progressive failure when the aspect ratio was 0.083 but the energy absorption of round tubes was better than square ones. The square geometries have catastrophic failure due to crack initiation in sharp edges, which in turn caused weakness in capability of the samples to absorb energy [3].
Mamali et al. [15] did an analytical study on semi apical 5°-20°conical shell samples and showed that the 5°semi apical samples have the highest impact resistance. Ansari et al. [13] proposed a new analytical approach to evaluate buckling deformation of conical composite shell based on minimum potential energy theory using Ritz solver. Hosseini et al. [11] conducted experimental tests on carbon epoxy conical structure to compare the energy absorption with glass epoxy structure. Mamali et al. [15] experimentally proved that conical structure with 0°layers was not suitable for energy absorption.
The numerical finite element method (FEM) is the key tool to predict behavior of structures with complex geometries under dynamic loading [5]. In order to observe failure mode under dynamic loading, many simulations were done using the finite element method, as was carried out by Luo et al. [16]. Kathiresan et al. [12] have done a holistic study on the effect of different vertex angles, types of glass fiber and showed that the specimen with woven fiber and 18°apex angle has the maximum peak load and energy absorption. Many studies have carried out in quasistatic loading through FE Modeling. Their studies focused on the crash worthiness indicators when a structure is under axial loading. Moreover the effect of filler material to enhance material properties has been investigated [17][18][19][20][21][22].
Topology optimization is used to initiate concept design and then the feasible concept can be modified to obtain the final product. This process reduces manufacturing cost and time in reaching the feasible design. Coehlo et al. [23] executed a hierarchical approach regarding design domain and material properties. Sohouli et al. and Zuo et al. performed a topology optimization for beam cross section laminates using minimum compliance criteria [24,25]. Some studies [26,27] have been performed discrete topology optimization method to specify the fiber orientation for multi-material composite laminate in order to get the best performance of the structure. The topology optimization was used for minimizing the strain energy in different materials, hence it can be used as a guideline to optimize the composite laminated shells.
In the present study, topology optimization was exploited to modify a conical geometry under quasi-static loading. The method yielded a design of conical composite structure which has the highest specific energy absorption capability.

Optimization method
The mathematical constraint for a general optimization problem requires the minimum and maximum value of an objective function that is subjected to equality or inequality constraint. h and g are design constraint functions expressed in terms of the vectorp ð Þ.p ð Þ is a vector of n independent variables of the optimization parameters. The minimum and maximum values of the constraints can be assumed as the peak values of the constraint functions. The constraint of the optimization process is defined as: To optimize laminate structures, the function op ð Þ and the values of the constraints h kp ð Þ and g lp ð Þ are computed through FE analysis. In this case the ABAQUS Topology Optimization Module (ATOM) is used as an optimization engine [28]. Figure 1 shows the flowchart of topology optimization deployed in ABAQUS.
From a previous study, one of the most researched geometry in crashworthiness is conical frusta which has straight conical chord [12]. To optimize this geometry, a conical frusta with 20°apex angle has been chosen and minimizing the volume is selected as the objective. The upper and lower cone diameters and the height are defined as the optimization constraints. The volume reduction is set from 0% to 30% with increment of 5%. After 50 design cycles and 15% volume reduction, a feasible model was achieved. The topology optimization was based on properties of chopped strand mat lamina [29]. The final optimized design was modeled and transferred to the CAD program in STL format and the structure with smoothly curved chord was redesigned. Figure 2 illustrates the optimized design concept and final design.

Specimen preparation
Based on the final design, two metal mandrels were fabricated using CNC turning machine. Figure 3 shows the two metal mandrels used to fabricate two types of geometry. Extension cylinders were fabricated to be fitted to the ends of the mandrels for ease of fabrication of the composite cones.
For each cone geometry, three samples were fabricated using filament winding and hand lay-up methods. The filament wound samples were produced with woven roving glass fiber and epoxy with 70% fiber weight. After curing, the specimens were separated from the mandrels. The filament winding process is shown in Figure 4. Extension lengths were used at both ends of the mandrels to overcome the problem of fiber orientation continuity in these regions.
In the hand lay-up method, the E-glass fiber mat was used. The chopped strand mat was curved around the rotating mandrel and polyester resin was manually brushed onto the CSM mat. Figure 5 illustrates the hand lay-up process and finished specimen. Specimens with two different thicknesses were made. Table 1 shows the dimension of different samples used in this study. The specimens are categorized based on the shape of conical mandrel (CC), concave mandrel (CV), method of fabrication, hand lay-up (H) and filament winding (F).
After finishing the fabrication process the specimens were controlled. The prepared samples based on ASTM D 3171-99-G, the resin of three samples was burned off and the volume fractures of the fabricated samples with ±2% were controlled. Moreover, a CMM machine was used to measure the thickness of different samples. In filament winded samples due to different tension in different sides of samples the thickness is varied.

Setup of experiment
A 600 kN Instron universal testing machine (UTM) was used to perform quasi static compression test. The loading rate was set at 2.5 mm/min. For repeatability verification, each test was repeated three times. Figure 6a shows the setup of the experiment performed on the universal testing machine. Moreover, Figure 6b demonstrates the progressive failure during the quasi static loading. Figure 7 shows the load-displacement graphs from the test data. The mean results from three test samples for each category were used in subsequent analysis.

Finite element modeling
For comparison of the behavior of different geometries under quasi-static loading, the hand lay-up samples were used. Comparison of crashworthiness capabilities between   two geometries were analyzed qualitatively. The geometries in consideration were the optimized geometry and regular conical specimen. The 3D model was created in ABAQUS/Explicit using material properties shown in Table 2. The four-node general-purpose shell element S4R was allocated to the deformable specimen and the rigid element was assigned to the top and bottom loading plates. A friction coefficient of 0.15 was set between the top and bottom loading plates and the specimen. To ensure that the mesh size is appropriate for this simulation, a mesh convergence study was carried out. In order to mesh convergence study, a unit axial load was applied on the top of specimen and by using Lagrangian method (H Method) the size of elements was decreased until reached a constant corresponding stress [30]. The satisfactory convergence was achieved when element size of 1.5 mm was utilized. The constant level of strain energy also can be a criteria for mesh quality verification. Figure 8 shows the mesh convergence study of a conical composite structure. For quasi-static loading, the explicit approach was used, whereby the internal energy should be less than 5% of the kinetic energy. Displacement control was applied where the specimens were compressed by 0 to 50 mm, as practiced during experimental tests. The rigid lower plate was constrained for all degrees of freedom and the top plate was free to move vertically. The finite element model of conical shell is shown in Figure 9.

Damage initiation and failure criteria
To demonstrate the damage initiation and progressive failure, Hashin's criteria was applied in ABAQUS where the onset of damage is related to material degradation. Hashin's failure prediction is evaluated by four different criteria, namely, fiber failure in tension, matrix failure in      tension, fiber failure in compression and matrix failure in compression [31]. When all these criteria are not reached, the structure is safe but when one or more of these equations are satisfied, damage is initiated. The element will be then degraded and finally the element will be removed from the calculation. Hashin's failure criteria model is shown in equations (4)- (7). 7 Results and discussion

Crashworthiness characteristic evaluation
The axial crushing test aims to obtain a controllable crushing pattern which can relate to crashworthiness capability. In this present study, crashworthiness indicators such as energy absorption (EA), peak force (P max ), mean crush force (P ave ), specific energy absorption (SEA) and crush force efficiency (CFE) have been addressed to evaluate crashworthiness capability of different cone geometries. The energy absorption is the criteria to demonstrate the stability limit of collapse when comparing between different shapes and geometries [32]. The magnitude of EA can be calculated by integration of the area under the load-displacement curve, and is shown in equation (8).
where P is the applied load and d and d are displacement and crushing deformation, respectively. In another form, where P represents the instantaneous load and d and H are displacement and crushed specimen length, respectively.
In every loading process, three stages of collapse are observed for all categories of specimens, namely, pre buckling, buckling, and post buckling. In the pre buckling stage, the laminates resist bending and after reaching the first peak load, the buckling stage is started followed by a  drop in load. Upon further loading, the maximum load P max can be obtained. Due to non-symmetric collapse in composite structures, it is possible that P max be different from the first peak load. Figure 11 shows the deformation and comparison of different geometry crushing responses. The mean crushing load P ave is calculated from the average force of the collapsing modes throughout the test. If the crushing load increases the crashworthiness ability can be increased as given by equation (9).
The mean crushing load P ave for a given crushed tube is calculated by: Specimens with higher P ave have more capacity to absorb energy from a controlled deformation. The structure is crushed but without catastrophic failure.
Energy absorption per volume, SEA is considered a vital indicator to get a better understanding about the crashworthiness of a structure, and is formulated as follows: where m is the mass of the specimen. Enhancement of SEA leads to increasing of crash resistance capacity. The other important indicator for assessing crashworthiness performance is CFE. High value of this indicator leads to more ideal energy absorption. CFE can be presented as follows [17]: Table 3 shows the crashworthiness indicators obtained from experiments and FE simulation. In all cases, it can be seen that the toughness of concave samples are higher than conical samples as indicated by the area under the load-displacement curve for each component. By changing the geometry from conical to concave, the energy absorption is enhanced by 32%. In Figure 11b, the mean experimental results are considered and other sample results are shown via the error bars. Figure 14a and b show the comparison of SEA values for different components. Regardless of thickness, it can be said that the optimized shape of every sample can result in up to 27% SEA enhancement. The FE method shows a maximum difference of 18% compared with test results and in most cases the difference is around 12%. This difference may be related to manufacturing issues and other damage criteria that were not considered in the simulation. Figure 15a and b show P max for different samples experimentally and numerically. Figure 15a shows the effect of thickness on energy absorption at the maximum peak load. Although the manufacturing process is different, increasing the thickness by 50% can raise the peak load by 80%. However, the dominant effects on all crashworthiness indicators are related to geometry of specimen.
The CFE illustrates the progressive collapse ratio. If the ratio approaches 100%, it can be interpreted that the progressive collapse occurs perfectly. Figure 16a and b show the comparison of CFE between numerical simulation and experiment tests. In contrast with metal specimen [31], composite structure shows near ideal progressive failure when P ave is close to P max . The effect of geometry can be seen as the dominant effect on CFE crashworthiness indicator. The maximum error between numerical and experimental analyses is 15%. The CFE indicator between CC-H-5.4 and CV-H-5.4 samples shows a difference of 7% but the result between CC-H-4.2 and CV-H-4.2 samples shows a difference of 14%. The difference between FE and experiment results can be attributed to localized stiffness which causes increase of the maximum peak load and subsequently affects the CFE value. The error bars show that other samples behave differently too. This phenomenon can be attributed to the manufacturing process or nonuniform volume fraction at specific cross section of the specimen which in turn causes the increase of the maximum peak load.

Collapse modes
The mode of collapse depends on factors related to the arrangement of the fibers, the properties of the matrix and fibers of the composite material as well as the specimen geometry [4]. For all present geometries, progressive crushing can be clearly observed in the form of "mushrooming". Collapse initiates at the smaller diameter end of the tube where delamination occurs in circumferential area (see Fig. 17). The inner layers of the laminate mushroom inwards (internal frond crushing zone a) and the outer layers mushroom outwards (external frond crush zone c) while some areas remain in debris wedge ends (debris wedge crush zone b) (see Fig. 18).

Mechanism of damage
The fracture process is governed by the properties of the fibers, resins, fiber-matrix interfaces and the arrangement and distribution of the fibers in the matrix [33]. All samples  in this experiment show controlled progressive failure modes. Progressive crushing in composite structures is more sensitive to delamination [15]. By increasing the compression load progressively, matrix cracking, crack initiation between layers and circumferential delamination start to appear (Fig. 18). This is accompanied by a drop in the load carrying as indicated on the load-displacement graph as a decrease in load after the first peak point. The cracks propagate between plies in the crush region of the tube causing the lamina to bundle. When a critical stress value is reached, the laminate proceeds to bend or buckle and creates different crushing zones.
For failure, both geometries show axial cracks which cause the formation of petals during progressive crushing. For concave and conical CSM samples, failure is accompanied by a combination of fiber splaying and axial tearing.
From Figure 18, region a in the concave structure is formed by layers bending in inflated curve shape but for conical geometry, the crushed region bends inside and down. Comparing the mechanism of damage for the two geometries, the concave geometry needs more load to bend the layers inside so this geometry has resulted in higher peak crush load, the corresponding mean crush load and eventually, higher total absorbed energy.
The major difference of energy absorption between conical and concave geometries lies in the direction of layers folding. The shape of internal frond and external frond zones for the conical geometry is a downward bend but for concave geometry, these crushing zones bend upwards (see Fig. 19).

Conclusion
A crashworthiness analysis was carried out on two different geometries, namely conical and concave. Using topology optimization process, the conical geometry was optimized, resulting in a concave geometry with 15% volume reduction. Both types of specimen were then fabricated using different manufacturing processes. The crashworthiness indicators such as maximum peak load, mean crush load, crash force efficiency, energy absorption and specific energy absorption were addressed. The results between experimental and numerical methods show the good agreement between them. The SEA was the most important indicator for crashworthiness evaluation and the optimized concave geometry exhibits a 32% increment in SEA. The major amount of the crushing energy was absorbed due to the increasing number of axial cracks and subsequent bending of crushed regions. The study reveals that the shape of structure has an important role in crashworthiness behavior due to damage distribution patterns in different regions.