Numerical simulation for expansion of preform and optimization of preform in thermoset composites

Abstract

The expansion of preform and the optimization of preform have become important steps in the molding process. At present, there are some questions in the expansion of thermoset composite material preform and precompression, for example, the inaccurate dimensions, cracks, and wrinkles. For the expansion of preform, the finite element inverse algorithm is used as the expansion algorithm, and then the initial solution is optimized by the arc length mapping method, the expansion of preform is realized by the iterative equation which is solved by the ABAQUS solver. The effectiveness of the expansion of preform is verified through the comparison between the finite element inverse algorithm with DYNAFORM. The optimization of the precompression process is researched in order to solved the problems of cracks and wrinkles in the integral precompression method of preform. Firstly, the precompression sequence is adjusted by the precompression method, and then the precompression direction is optimized by the genetic algorithm. Through numerical simulation, the maximum thinning rate is reduced to 13%, and the maximum thickening rate is reduced to 6%, which improve the problems of cracks and wrinkles of preform, and the effectiveness of the optimization method is verified.

Introduction

In the molding process, the quality of composite products can be improved by preforming the prepreg and then molding. First of all, the prepreg needs to be expanded into a certain shape of blank before preforming, the phenomenon of crack, and wrinkle during preforming can be reduced by the reasonable blank shape, as well as the trimming workload will be reduced after molding. ¹ Then, in order to improve product quality and improve product production efficiency, the prepreg is pressed into a certain shape through preforming. The quality of preform will be affected by the factors such as the direction of precompression when preforming.

The blank expand method of preform could be divided into two types: the geometric method and the mechanical method. In the geometric method, the geometric relationship between the expanded shape and the preform was only considered. The common geometric methods are: analytical method, surface parameterization method, ABF method, discrete Gaussian curvature to parameterize the surface, Gaussian curvature, and geodesic principle combined, complex surface development based on mesh edges, etc.^2,3 The influence of the constitutive equation of the stress-strain relationship of the prepreg on the expansion of the blank was not considered in the geometric method, and the result does not match the actual situation exactly. In the mechanical method, the constitutive equation of the raw material is applied to the blank expansion. Kim et al. ⁴ proposed a rigid-plastic finite element method based on the variational principle to calculate the relationship between stress and strain at each point. Yang et al. ⁵ combined with Kim’s idea and proposed an elastic-plastic finite element method, transformed the stress-strain matrix, derived the elastoplastic stress-strain matrix, and promoted the application and development of mechanical preform expansion. Liang and Bin ⁶ proposed a constrained nonlinear optimization method for the expansion of preform, but the global optimal results could not be got. Azariadis et al. ⁷ proposed an improved algorithm for the global optimization expansion of any blank surface, and the optimized expansion of preform was obtained by genetic algorithm. However, most of the current researches focus on the blank expansion of metal materials or thermoplastic prepreg, and the research on the blank expansion of thermoset composite prepreg is still very few.

The conventional method for solving the expansion of the preform is the finite element incremental method, but it has many preconditions and its use is restricted. The finite element inverse algorithm only considers the initial configuration and the final configuration, which increases the calculation speed. The inverse finite element algorithm requires an initial solution and iterative solution. At present, the vertical projection method and the spherical projection method are often used to determine the initial solution. However, the projection area will be zero when the projection point is not selected properly during the calculation in these two methods which determine the initial solution, and the effective solution can be obtained.

The preform is got through the precompression of multi-layer prepreg after blank expansion. The process of precompression can be divided into single-process precompression, compound mold precompression and continuous mold precompression. The traditional precompression design mainly relies on the designer’s experience, various precompression parameters are modified repeatedly, and it is time-consuming, high-cost, and difficult to adapt to the requirements of modern industry. ⁸ At present, many researches have been carried out on the forming law of preforming during precompression, and numerical simulation has been used in the field of molding widely.

First of all, the thickness and wrinkles of the preform were affected by the pressure and temperature of preform. Ivanov and Lomov ⁹ calculated the law of prepreg thickness change with pressure from the perspective of mechanics, and explained the thickness deformation mechanism of prepreg during compression. Ho et al. ¹⁰ reduced forming defects through improving the mold and considering the mechanical properties of laminated composite materials with different thicknesses. Somashekar et al. ¹¹ researched the characteristics of the compression and rebound of glass fiber preforms, and found that the rebound characteristics of preform were related to the maximum pressure, the rebound thickness was unchanged when the pressure exceeds the critical value. Sun et al. ¹² found that the generation of wrinkles in laminates was related to the friction mechanism closely, and the friction mechanism could be changed by increasing the preforming pressure.

Furthermore, the relationship between the layup sequence of the prepreg and the mechanical properties of the preform was studied. Wang et al. ¹³ simulated and analyzed the tensile properties of several preforms with different layup orders by using finite element theory, and found that the influence of the layup order on the rigidity of the preform. Lee et al. ¹⁴ evaluated the effect of different stacking sequences of prepregs on mechanical properties, then the validation of bending stiffness prediction model was performed by comparing the load per stress roller stroke through three-point bending analysis and test.

In addition, the cause of the defect was clarified by analyzing the interlayer change mechanism of the precompression prepreg through the numerical simulation of preforming. Wang et al. ¹⁵ analyzed the in-plane shear behavior of thermoset prepregs during preforming, and proposed that the relationship between stress and strain would change under different preforming speeds. Wang et al.^16,17 studied the interfacial behavior between prepreg laminates, and established a finite element model of prepreg preforming which considered the effects of preforming pressure, then the preforming experiments of the prepreg laminates were implemented and simulated by using the proposed FE mode, comparisons between experimental results and simulation data verified the numerical model. Boisse et al. ¹⁸ proposed the second gradient method to describe the generation of wrinkles based on the research of slippage phenomenon in the preforming process of prepreg, and simulated it with meso finite element models used for macroscopic forming. Zhang et al. ¹⁹ established a interaction model through the research of influence of the precompression direction on the interaction between different prepreg laminates during prepreg preforming, and proposed a comprehensive computational material engineering method, then the preforming process was simulated from the micro and macro perspectives.

Different preforming processes and optimization of process parameters have been studied, and the preforming winkles, mechanical properties and other indicators have been improved. Sorrentino and Bellini ²⁰ researched the hot drape forming process for preforming, which improved the wrinkles and irregular thickness. Kim et al. ²¹ evaluated the effect of different molding conditions on the mechanical properties of composite materials, and the optimized molding process parameters were obtained after experimental comparison. Di et al. ²² proposed that the maximum initial contact area was the main criterion for establishing the precompression direction based on the grid data of the cover, and established the algorithm model with the direction vector as a variable, then developed the corresponding module which could determine the precompression direction quickly.

However, there are many analyses on the preforming mechanism in the preforming process of prepregs, and there are few applied researches on the process optimization methods of the preforming process, especially the multi-objective optimization of prepreg preforming quality has not been studied.

In this paper, firstly, a finite element analysis model for composite material prepreg expansion is established. The mechanical expansion method based on the constitutive equation of the prepreg is used, combined with the inverse finite element algorithm to solve the equation, and then the arc length mapping method is used to optimize the initial solution which is obtained by the vertical projection method, and the numerical simulation of the composite material prepreg expansion is carried out. The effectiveness of the finite element inverse algorithm for composite material prepreg deployment is verified through theoretical analysis, numerical simulation, and comparison. Then, after the expansion of preform, multiple precompression regions of the preform will be selected, the precompression direction is used as the process parameter, the overall thinning rate of the multiple regions is used as the optimization target, and the genetic algorithm is used to optimize the precompression direction, the effectiveness of the optimization results are verified by the simulation.

Materials and methods

Materials

In this paper, glass fiber/epoxy resin prepreg is used to research, some structural parameters of the prepreg is shown in Table 1, including some parameters are provided by Taian Joy Composite Material Technology Co., Ltd., Taian, China.

Table 1.

Parameters of prepreg.

Weight ( $\mathrm{g}/{\mathrm{m}}^2$ )	Number of fibers in the fiber bundle	Resin content ( $\mathrm{vol}\mathrm{\%}$ )	Fiber density ( $\mathrm{g}/{\mathrm{m}}^3$ )	Number of fiber bundles ( $/\mathrm{cm}$ )
520 ± 20	3000	54 ± 4	2.55	5
Gap length ( $\mathrm{mm}$ )	Fiber bundle width ( $\mathrm{mm}$ )	Prepreg fiber bundle width (mm)	Thickness ( $\mathrm{mm}$ )	Fiber diameter ( $\mu \mathrm{m}$ )
0.32	1.7	0.11	0.22	9

Mathematical model of expansion

The basic principle of the finite element inverse algorithm

Lee and Huh ²³ proposed the finite element inverse algorithm. Assuming that the preform process of the blank is loaded proportionally, the initial shape and final configuration of the blank are considered only, and the intermediate deformation process is not considered. Figure 1 is a schematic diagram of the finite element inverse algorithm.

Figure 1.

Schematic diagram of finite element inverse algorithm.

In Figure 1, the blank A⁰ becomes the final configuration A under the action of the mold force, and the final configuration A is discretized triangularly. The position B⁰ of the node B in the blank A⁰ is obtained by the finite element method, and B is the final configuration at any point on the above, the initial blank shape of the preform and the thickness distribution after precompression can be obtained by comparing B and B⁰.

The configuration shape of the precompression product and the initial thickness of the blank are known conditions, and the geometric shape of the preform deployment and the thickness distribution of the final configuration are obtained. Through the known conditions and using the plastic work equation to obtain the unknown conditions, the desired geometry of the preform expansion can be obtained. The known and unknown quantities in the finite element inverse algorithm are shown in Table 2.

Table 2.

Known quantity and unknown quantity of finite element inverse algorithm.

Configuration	Known quantity	Unknown quantity
Initial configuration A0 (preform expansion)	Initial thickness	Shape and size
	Initial stress and strain	Initial node coordinates
Final configuration A (precompression product)	Final configuration shape and coordinates	Final configuration thickness
	Resistance and direction	Final stress and strain
	Node vertical displacement	The force of the mold

The equation of the finite element inverse algorithm

Blank motion description and strain equation

When t = 0, the configuration where the object is not deformed or does not move is the initial configuration. When t > 0, it is called the current configuration. In the theoretical research, the Lagrange coordinate system (Matter coordinate system) can be introduced to determine the position of the initial configuration. In order to specifically express a point in the initial configuration, it is represented by X( ${\mathrm{X}}_1$ , ${\mathrm{X}}_2$ , ${\mathrm{X}}_3$ ). And the Euler coordinate system is used to determine the position of the current configuration, any point of which is represented by x( ${\mathrm{x}}_1$ , ${\mathrm{x}}_2$ , ${\mathrm{x}}_3$ ).

As shown in Figure 2, the Lagrange coordinate system is combined with the Euler coordinate system, and the current configuration is obtained from the initial configuration through a certain time t.

Figure 2.

Two configurations at t = 0 and t > 0.

Combining the nodes in the graph with time t, the expressions of positions ${\mathrm{x}}_{\mathrm{i}}$ , ${\mathrm{X}}_{\mathrm{i}}$ and time t are shown as:

$$\mathrm{x}=\mathrm{x}(\mathrm{X},\mathrm{t})$$

(1)

The displacement of the node which is required to move from the initial configuration to the current configuration, $\mathrm{u}$ can be refined as:

$$\mathrm{u}=\mathrm{x}-\mathrm{X}$$

(2)

In the initial configuration, there are two nodes $\mathrm{X}$ and $\mathrm{X}+\mathrm{dX}$ are supposed that these are infinitely close, and the distance between them is the line element $\mathrm{dX}$ , which becomes $\mathrm{dx}$ in the current configuration, from equation (1), then the following equation can be obtained as:

$$\mathrm{x}+\mathrm{dx}=\mathrm{x}(\mathrm{X}+\mathrm{dX},\mathrm{t})$$

(3)

The distance between the two nodes in the current configuration, $\mathrm{dx}$ , can be calculated as:

$$\mathrm{dx}=\mathrm{x}(\mathrm{X}+\mathrm{dX},\mathrm{t})-\mathrm{x}(\mathrm{X},\mathrm{t})$$

(4)

The Taylor series expansion at $\mathrm{X}$ in the above equation is omitted, and the higher-order writing matrix form is omitted, $\mathrm{dx}$ , can be calculated as:

$$\mathrm{dx}=\left[\begin{array}{c}\mathrm{d}{\mathrm{x}}_1\hfill \\ \multicolumn{1}{c}{\mathrm{d}{\mathrm{x}}_2}\\ \multicolumn{1}{c}{\mathrm{d}{\mathrm{x}}_3}\end{array}\right]=\left[\begin{array}{ccc}\frac{\partial {\mathrm{x}}_1}{\partial {\mathrm{X}}_1}\hfill & \frac{\partial {\mathrm{x}}_1}{\partial {\mathrm{X}}_2}\hfill & \frac{\partial {\mathrm{x}}_1}{\partial {\mathrm{X}}_3}\hfill \\ \multicolumn{1}{c}{\frac{\partial {\mathrm{x}}_2}{\partial {\mathrm{X}}_1}}& \frac{\partial {\mathrm{x}}_2}{\partial {\mathrm{X}}_2}\hfill & \frac{\partial {\mathrm{x}}_3}{\partial {\mathrm{X}}_3}\hfill \\ \multicolumn{1}{c}{\frac{\partial {\mathrm{x}}_3}{\partial {\mathrm{X}}_1}}& \frac{\partial {\mathrm{x}}_3}{\partial {\mathrm{X}}_2}\hfill & \frac{\partial {\mathrm{x}}_3}{\partial {\mathrm{X}}_3}\hfill \end{array}\right]\left[\begin{array}{c}\mathrm{d}{\mathrm{X}}_1\hfill \\ \multicolumn{1}{c}{\mathrm{d}{\mathrm{X}}_2}\\ \multicolumn{1}{c}{\mathrm{d}{\mathrm{X}}_3}\end{array}\right]=\mathrm{FdX}$$

(5)

The matrix $\mathrm{F}$ is a linear transformation matrix, which is called the deformation gradient tensor, it reflects the deformation and motion of the line element. The deformation gradient contains the complete information about the deformation of each triangulated element after discretization.

The Lagrange coordinate system is taken as an independent variable, from equation (5), the lengths of line elements, $\mathrm{dX}$ and $\mathrm{dx}$ , are obtained as:

$$\{\begin{array}{c}\mathrm{d}{\mathrm{L}}^2={\left[\mathrm{dX}\right]}^{\mathrm{T}}\left[\mathrm{dX}\right]\hfill \\ \mathrm{d}{\mathrm{l}}^2={\left[\mathrm{dX}\right]}^{\mathrm{T}}{\mathrm{F}}^{\mathrm{T}}\mathrm{F}\left[\mathrm{dX}\right]\hfill \end{array}$$

(6)

The right Cauchy-Green deformation tensor is defined as $\mathrm{C}={\mathrm{F}}^{\mathrm{T}}\mathrm{F}$ , C is a symmetric positive definite tensor, since the square of the line element is taken as a positive value. C is only related to the deformation of the object, not to the rotation of the rigid body.

The Euler coordinate system is taken as the independent variable, from equation (5), the lengths of the line elements, $\mathrm{dx}$ and $\mathrm{dX}$ , are obtained as:

$$\{\begin{array}{c}\mathrm{d}{\mathrm{l}}^2={\left[\mathrm{dx}\right]}^{\mathrm{T}}\left[\mathrm{dx}\right]\hfill \\ \mathrm{d}{\mathrm{L}}^2={\left[\mathrm{dx}\right]}^{\mathrm{T}}{\mathrm{F}}^{-\mathrm{T}}{\mathrm{F}}^{-1}\left[\mathrm{dx}\right]\hfill \end{array}$$

(7)

The left Cauchy-Green deformation tensor is defined as $\mathrm{B}=\mathrm{F}{\mathrm{F}}^{\mathrm{T}}$ , similar to C, $\mathrm{B}$ is also a symmetric positive definite tensor.

Assuming that any amount of deformation in the current configuration is expressed as elongation $\lambda =\mathrm{dl}/\mathrm{dL}$ , from equation (7), ${\lambda }^{-2}$ is obtained as:

$${\lambda }^{-2}=\frac{{\left[\mathrm{dx}\right]}^{\mathrm{T}}}{\mathrm{dl}}{\mathrm{B}}^{-1}\frac{\left[\mathrm{dx}\right]}{\mathrm{dl}}={\mathrm{N}}^{\mathrm{T}}{\mathrm{B}}^{-1}\mathrm{N}$$

(8)

From equation (8), $\mathrm{B}$ is used to represent the deformation in the three main deformation directions. Let $\left[{\mathrm{N}}_1\right]$ , $\left[{\mathrm{N}}_2\right]$ , and $\left[{\mathrm{N}}_3\right]$ be the unit vectors in the three principal deformation directions, and ${\lambda }_1$ , ${\lambda }_2$ , and ${\lambda }_3$ are the principal elongations in the three principal deformation directions, a matrix of column vectors, $\mathrm{M}$ , is obtained as:

$$\left[\begin{array}{ccc}{\lambda }_1^2\hfill & \hfill & \hfill \\ \multicolumn{1}{c}{}& {\lambda }_2^2\hfill & \hfill \\ \multicolumn{1}{c}{}& \hfill & {\lambda }_3^2\hfill \end{array}\right]={\mathrm{M}}^{\mathrm{T}}{\mathrm{B}}^{-1}\mathrm{M}$$

(9)

${\lambda }_1$ , ${\lambda }_2$ , ${\lambda }_3$ , and $\left[{\mathrm{N}}_1\right]$ , $\left[{\mathrm{N}}_2\right]$ , $\left[{\mathrm{N}}_3\right]$ are the eigenvalues and eigenvectors of ${\mathrm{B}}^{-1}$ .

In the preforming process of the blank, strain is introduced to describe the deformation of the blank. The final configuration is discretized into multiple triangular cells by using triangular unit cells. The complex equations of the entire region can be estimated by the deformation analysis of the triangular cells. As shown in Figure 3, each unit belonging to the final configuration is mapped to the initial configuration unit. The strain of the entire element can be obtained by calculating the amount ${\mathrm{G}}_{\mathrm{i}}$ of the initial configuration element from the knowing edge vector $g_i$ in the final configuration element.

Figure 3.

Element mapping from final state to initial state.

The equation for obtaining unit strain is known as:

$${\mathrm{\epsilon }}_{\mathrm{ij}}=\mathrm{M}\left[\begin{array}{ccc}\ln {\lambda }_1\hfill & \hfill & \hfill \\ \multicolumn{1}{c}{}& \ln {\lambda }_1\hfill & \hfill \\ \multicolumn{1}{c}{}& \hfill & \ln {\lambda }_1\hfill \end{array}\right]\mathrm{M}$$

(10)

Let $\mathrm{\theta }$ be the angle between the main elongation $\lambda$ and the $\mathrm{X}$ axis of the final configuration coordinate system, from equations (9) and (10), the strain which in the two states from the initial configuration of the triangular element to the final configuration, $\epsilon$ , is calculated as:

$$\epsilon =\left[\begin{array}{c}{\mathrm{\epsilon }}_{\mathrm{x}}\hfill \\ \multicolumn{1}{c}{{\mathrm{\epsilon }}_{\mathrm{y}}}\\ \multicolumn{1}{c}{{\mathrm{\epsilon }}_{\mathrm{xy}}}\end{array}\right]=\left[\begin{array}{c}\ln {\lambda }_1\mathrm{co}{\mathrm{s}}^2\theta +\ln {\lambda }_2\mathrm{si}{\mathrm{n}}^2\theta \hfill \\ \multicolumn{1}{c}{\ln {\lambda }_1\mathrm{si}{\mathrm{n}}^2\theta +\ln {\lambda }_2\mathrm{co}{\mathrm{s}}^2\theta }\\ \multicolumn{1}{c}{\ln \frac{{\lambda }_1}{{\lambda }_2}\cos \mathrm{\theta }\sin \mathrm{\theta }}\end{array}\right]$$

(11)

Constitutive relationship

The constitutive relationship represents the relation between stress and strain. In this paper, the epoxy glass fiber prepreg is used for preforming. It can be expressed as:

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{\sigma }=-\mathrm{pI}+{\mathrm{T}}_{\mathrm{q}}\mathrm{m}+{\mathrm{T}}_{\mathrm{r}}\mathrm{n}+{\mathrm{\beta }}_1\mathrm{D}+{\mathrm{\beta }}_2{\mathrm{D}}^2}\\ +{\mathrm{\beta }}_3(\mathrm{Dm}+\mathrm{mD}+\mathrm{Dn}+\mathrm{nD})\hfill \\ +{\mathrm{\beta }}_4(\mathrm{Dp}+\mathrm{pD}+\mathrm{D}{\mathrm{p}}^{\mathrm{T}}+{\mathrm{p}}^{\mathrm{T}}\mathrm{D})\hfill \\ +{\mathrm{\beta }}_5({\mathrm{D}}^2\mathrm{m}+\mathrm{m}{\mathrm{D}}^2+{\mathrm{D}}^2\mathrm{n}+\mathrm{n}{\mathrm{D}}^2)\hfill \\ +{\mathrm{\beta }}_6({\mathrm{D}}^2\mathrm{p}+\mathrm{p}{\mathrm{D}}^2+{\mathrm{D}}^2{\mathrm{p}}^{\mathrm{T}}+{\mathrm{p}}^{\mathrm{T}}{\mathrm{D}}^2)\hfill \end{array}$$

(12)

Where $\mathrm{I}$ is the unit tension, p is any hydrostatic pressure, ${\mathrm{T}}_{\mathrm{q}}$ and ${\mathrm{T}}_{\mathrm{r}}$ are the pressure in the warp and weft directions, m is the warp tensor, n is the weft tensor, and $\mathrm{D}$ is the deformation rate tensor, The deformation rate tensor can be obtained from the right Cauchy-Green deformation rate tensor corresponding to the right Cauchy-Green tensor. $\mathrm{It}$ can be expressed as:

$$\frac{1}{2}\stackrel{·}{\mathrm{C}}={\mathrm{F}}^{\mathrm{T}}\mathrm{DF}$$

(13)

In this constitutive equation (12), there are six parameters, ${\mathrm{\beta }}_1$ and ${\mathrm{\beta }}_2$ are the shear viscosities of the resin material, ${\mathrm{\beta }}_3$ and ${\mathrm{\beta }}_5$ are the viscosities related to the warp and weft fiber directions, called longitudinal viscosity and transverse viscosity, and ${\mathrm{\beta }}_4$ and ${\mathrm{\beta }}_6$ are the viscosity related to the coupling direction of the weft fiber bundles is only effective when shear deformation occurs.

In this paper, epoxy resin, ${\mathrm{\beta }}_1$ and ${\mathrm{\beta }}_2$ can be determined by rheological experiments. Wang and Du ²⁴ tested the resin material used in the prepreg fabric, and obtained two parameters with a value of 406.84 Pa•s. Assuming that ${\mathrm{\beta }}_3$ and ${\mathrm{\beta }}_5$ are equal, Harrision et al. ²⁵ obtained the calculation equationo. ${\mathrm{\beta }}_3$ is shown as:

$${\mathrm{\beta }}_3={\mathrm{\beta }}_{\mathrm{m}1}\left[\frac{{\mathrm{d}}_{\mathrm{f}}({\mathrm{d}}_{\mathrm{f}}+\mathrm{g})}{\mathrm{g}({\mathrm{d}}_{\mathrm{f}}+{\mathrm{g}}_0)}\right]+{\mathrm{\beta }}_{\mathrm{m}2}[1-\frac{{\mathrm{d}}_{\mathrm{f}}\cos \mathrm{\gamma }}{{\mathrm{d}}_{\mathrm{f}}+{\mathrm{g}}_0}]$$

(14)

${\mathrm{\beta }}_{\mathrm{m}1}$ and ${\mathrm{\beta }}_{\mathrm{m}2}$ are the resin viscosity between the fiber bundle before and after deformation, $\mathrm{\gamma }$ is the shear angle, ${\mathrm{d}}_{\mathrm{f}}$ is the diameter of the fiber bundle, ${\mathrm{g}}_0$ is the initial width of the resin without deformation, and $\mathrm{g}$ is the width of the resin after deformation.

$$\mathrm{g}=({\mathrm{d}}_{\mathrm{f}}+{\mathrm{g}}_0)\cos \mathrm{\theta }-{\mathrm{d}}_{\mathrm{f}}$$

(15)

For coupling viscosities ${\mathrm{\beta }}_4$ and ${\mathrm{\beta }}_6$ , the calculation equation is shown as:

$${\mathrm{\beta }}_4={\mathrm{\beta }}_6=2{\mathrm{\beta }}_3\cos (\frac{\mathrm{\pi }}{4}-\frac{\mathrm{\gamma }}{2})$$

(16)

Minimum plastic work and iterative equation

Chung and Richmond ²⁶ combined the deformation theory with the extreme value work path and proposed the ideal deformation. The ideal deformation can be applied to the process of the preform and the finite element inverse algorithm program will be developed. Considering the deformation process as a uniform speed process, the incompressible condition of plastic deformation is also considered, the plastic work of uneven deformation is expressed as:

$$\mathrm{W}=\int (\int {\mathrm{\sigma }}^{*}\mathrm{Ddt})\mathrm{d}{\mathrm{V}}_0=\int \left({\overline{\mathrm{\sigma }}}^{*}\mathrm{d}\overline{\mathrm{\epsilon }}\right)\mathrm{d}{\mathrm{V}}_0$$

(17)

Here, $\mathrm{\sigma }$ is the stress, $\mathrm{D}$ is the strain rate, $\overline{\mathrm{\sigma }}$ is the equivalent stress, $\overline{\epsilon }$ is the equivalent strain, and ${\mathrm{V}}_0$ is the initial volume. When the deformation path of the material unit is determined, the plastic work can only be determined by the displacement, which can be regarded as:

$$\mathrm{W}=\mathrm{W}(\mathrm{U})$$

(18)

And the relationship between the coordinates of the initial configuration and the final configuration is:

$$\mathrm{U}=\mathrm{x}-\mathrm{X}$$

(19)

Here, $\mathrm{x}$ is the known final configuration coordinate, $\mathrm{X}$ is the sought initial configuration coordinate, and $\mathrm{U}$ is the displacement. Therefore, the plastic work can be written as a function whose initial coordinate $\mathrm{X}$ is an independent variable, and from the extreme value of the overall plastic work, the static equilibrium equation can be obtained as:

$$\mathrm{R}\left(\mathrm{X}\right)=\frac{\partial \mathrm{W}\left(\mathrm{X}\right)}{\partial \mathrm{X}}=0$$

(20)

Here, $\mathrm{R}\left(\mathrm{X}\right)$ is the residual stress, the equation (20) is solved by Newton-Raphson method, and the iterative equation system is obtained as:

$${\left[\frac{\partial \mathrm{R}\left(\mathrm{X}\right)}{\partial \mathrm{X}}\right]}_i\mathrm{dX}=\mathrm{R}{\left(\mathrm{X}\right)}_{\mathrm{i}},{\mathrm{X}}_{\mathrm{i}+1}={\mathrm{X}}_{\mathrm{i}+1}+\mathrm{\beta }\mathrm{dX}$$

(21)

In the equation (21), $\mathrm{\beta }$ is the deceleration factor. In this paper, $\mathrm{\beta }$  = 1, the solution can be solved by bringing equations (11), (12), and (17) into the equation (21).

Initial solution and optimization

The nonlinear equation (21) is difficult to solve directly, so it is necessary to give an initial estimate ${\mathrm{X}}_0$ during the solution. The initial solution has a great influence on the iterative operation. If the estimated value is reasonable, it can not only speed up the solution, but also determine whether the final result is correct. In this paper, the vertical projection method is used to solve, and the arc length mapping method is used to optimize the initial solution.

The initial solution

Lee and Huh ²³ proposed to use linear mapping to map the final configuration to the initial plane, and then obtained the initial solution in the finite element inverse algorithm. However, the vertical projection method is used in the linear mapping method, which will affect the correctness of the initial solution and the solution of the inverse algorithm in special cases. Figure 4 is a schematic diagram of the vertical projection method.

Figure 4.

Schematic diagram of vertical projection.

The principle of the vertical projection method is used to project the discrete nodes in the final configuration vertically onto the initial plane along the precompression direction, so the initial solution can be obtained. There are many problems in practical applications. When a part of the element is $90°$ from the initial plane and parallel to the projection direction, the projected area of this part is zero. As shown in Figure 4, the projection of two nodes A and B which on the vertical wall, will cause the overlap phenomenon on the initial plane, and distortion of the element, then the calculation has a large error and cannot be performed.

The optimization of initial solution

In order to solve the above problems, the arc length mapping method is used to optimize the initial solution in this paper. The arc length mapping method can avoid the situation where the projected area of the element is zero, and improve the accuracy of the result. And the initial solution is calculated by the arc length mapping method is closer to the result obtained by the inverse algorithm, which can accelerate the convergence speed of the inverse algorithm. As shown in Figure 5.

• (1) O is the expansion center point, A is the selected starting point, B is the selected current node, the section perpendicular to the X–Y plane formed by the two points is P, and the arc length between the section and the final configuration is L.

• (2) The projection of the section on the X–Y plane is supposed as $\overline{V}$ , the angle between $\mathrm{V}$ with the Y axis is θ, and the arc length L is converted to the same length, then the position of point C as:

$$\overrightarrow{\mathrm{OC}}=\mathrm{L}\frac{\overline{V}}{\Vert \overline{V}\Vert }$$

(22)

• (3) Transform point B into all the nodes of the final configuration, and traverse all the nodes to get the initial calculate value ${\mathrm{X}}_0$ .

Figure 5.

Schematic diagram of arc length mapping method.

The solution can be solved by the Newton-Raphson method iteratively after obtaining the initial solution ${\mathrm{X}}_0$ . A criterion is needed to terminate the iteration when the iterative formula is solved. In this paper, residual stress is used as the criterion, because of the formed configuration will produce a damaged state when the residual stress reaches a certain range, which is not allowed. The scope given in this article is $\Vert \mathrm{\Delta }\mathrm{R}\Vert \le \epsilon ={10}^{-6}$ .

The optimization of precompression direction

For the expansion prepreg, the integral precompression process is used firstly, then the step-by-step precompression method is used to change the precompression sequence, and the precompression direction is optimized by genetic algorithm to improve the preform quality.

Evaluation function

The angle $\mathrm{\alpha }$ between the precompression direction and the X axis or Y axis is used as a variable to establish the objective function, when the precompression direction model is established. According to Figure 6, the objective function can be established.

Figure 6.

Section analysis schematic diagram.

The initial contact area as:

$${\mathrm{f}}_1\left(\mathrm{\alpha }\right)=\sum _{\mathrm{i}=1}^{\mathrm{n}}({\mathrm{Y}}_{\mathrm{ie}}-{\mathrm{Y}}_{\mathrm{is}})$$

(23)

The uniformity of the distribution of initial contact points as:

$$f_2\left(\mathrm{\alpha }\right)=\sum _{\mathrm{i}=1}^{\mathrm{n}}\frac{{\mathrm{Y}}_{\mathrm{i}}}{\mathrm{n}}-{\mathrm{Y}}_0$$

(24)

The degree of dispersion of the initial contact points as:

$${\mathrm{f}}_3\left(\mathrm{\alpha }\right)=\max \left\{{\mathrm{Y}}_{\mathrm{i}}\right\}-\min \left\{{\mathrm{Y}}_{\mathrm{i}}\right\}$$

(25)

The uniformity of precompression resistance as:

$${\mathrm{f}}_4\left(\mathrm{\alpha }\right)=|{\mathrm{Z}}_{\mathrm{e}}-{\mathrm{Z}}_{\mathrm{s}}|$$

(26)

Where $\mathrm{n}$ represents the number of segments in the B line crosses the section line, ${\mathrm{Y}}_{\mathrm{is}}$ , ${\mathrm{Y}}_{\mathrm{i}}$ , and ${\mathrm{Y}}_{\mathrm{ie}}$ represent the Y coordinate values of the starting point, midpoint, and end point of the i-th intersection of the B line and the section line, ${\mathrm{Y}}_0$ represents the coordinate value of the geometric center line, ${\mathrm{Z}}_{\mathrm{e}}$ and ${\mathrm{Z}}_{\mathrm{s}}$ are the Z coordinates of the start and end points of the geometric figure.

The evaluation function method can be used to transform multi-objective optimization problems into single-objective optimization. The evaluation function U(α) can be expressed as:

$$\mathrm{U}(\alpha )=\sum _{\mathrm{i}=1}^4{\mathrm{\omega }}_{\mathrm{i}}{\mathrm{U}}_{\mathrm{i}}$$

(27)

${\mathrm{\omega }}_{\mathrm{i}}$ is the weight coefficient, $\sum ^{{\mathrm{\omega }}_{\mathrm{i}}}=1$ , it can be set as the reciprocal of the optimal value of the corresponding sub-objective function.

Genetic algorithm optimization

There is a nonlinear relationship between the thinning rate of each region and the direction of precompression. In this paper, the genetic algorithm is used to optimize the direction of precompression, and to find the global optimum of the thinning rate of each region.

During the operation of the genetic algorithm, the feasible solutions are operated to achieve the purpose of optimization. The process is shown as in Figure 7.

Figure 7.

Genetic algorithm workflow.

Coding

The z axis is taked as the reference direction, the decision variables are the rotation angle ${\mathrm{\alpha }}_{\mathrm{x}}$ that around the x-axis and the rotation angle ${\mathrm{\alpha }}_{\mathrm{y}}$ that around the y-axis; the binary 8-bit coding string is used to represent the single decision variable, since the slope ranges from – $90°$ to $90°$ , a total of 16 bits.

Fitness

The value range of the function is always non-negative, and that the maximum value of the function can be solved is the optimization goal. Therefore, this paper takes the individual fitness directly as the corresponding evaluation function value.

Choose

The proportional selection operator is selected in operation.

Crossover

The single-point crossover operator is used in crossover operation. In each generation group, the number of individuals performing the crossover operation is determined by the crossover probability. Here, this parameter is defined as a constant.

Variation

The basic bit mutation operator is used in mutation operation. The probability of mutation here is defined as a constant.

The generation of the next generation group is generated after the three operations of selection, crossover and mutation for the individuals of the previous generation group. In addition, the situation of each generation of individuals is output, including the genotype and phenotype of each generation of individuals, the value of the evaluation function, the average value of the evaluation function of each generation, and the current optimal value. Then, the optimization of the stamping direction is completed finally.

Results and discussion

This section may be divided by subheadings. It should provide a concise and precise description of the experimental results, their interpretation as well as the experimental conclusions that can be drawn.

Expansion of battery box preform

Expansion of the finite element inverse algorithm

The battery box of the electric car is selected as the object. The battery box is 1700 mm long, 1000 mm wide, and the highest part has a height of 300 mm. Discrete the triangular unit of the battery box CAD model, as shown in Figure 8, which contains 76,633 nodes and 23,834 units, Figure 9 is an expansion of the battery box model.

Figure 8.

Discrete model of triangular element for battery box.

Figure 9.

(a) Expansion diagram of the finite element inverse algorithm and (b) expansion diagram of DYNAFORM.

The initial solution is plugged into the solver for calculation, and that the blank shape of the product is get only takes 15 s. It can be seen that the optimized initial solution algorithm based on the arc length mapping method is suitable for finding the final solution of the finite element inverse algorithm. The expansion of the preform is shown in Figure 9(a). In order to verify the accuracy of the inverse algorithm, Figure 9(b) is the expansion diagram of DYNAFORM software.

Analysis of expansion calculation of battery box perform

For clear comparison, the size difference between the two expansion methods are compared, and the larger part of the difference for calculation is intercepted, as shown in Figure 10. Through calculation, it is found that compared with the calculation by DYNAFORM software, the expansion of the blank by the finite element inverse algorithm is in a compact state, with a maximum reduction of 10.2 mm, and the overall reduction of prepreg is 5.45%, which reduces the cutting workload of composite materials in the later period.

Figure 10.

Shape comparison between finite element inverse algorithm and DYNAFORM.

Precompression process design

Integral precompression process

DYNAFORM is used to simulate the integral precompression process. The blank is a glass fiber epoxy prepreg with a thickness of 2 mm, a precompression of 2000 N, and a precompression direction of −Z. Figure 11 shows the thickness distribution after precompression.

Figure 11.

Thickness distribution diagram of integral precompression preforming method.

It can be seen that the thickness of the three points A, B, and C is between 0.56 and 0.95 mm, with a thinning rate of 52.5% −72%, the phenomena of cracks are occurred. The thickness of the three points D, E, and F is between 2.3 and 2.4 mm, which indicates that the wrinkles are generated.

Step-By-Step precompression process

Division of precompression area

In view of the above problems, a Step-By-Step precompression process is proposed. Multiple precompression at different positions is performed on the same station and the same mold. The sequence of the precompression is shown in Figure 12.

Figure 12.

Step-By-Step precompression sequence diagram.

There are 11 cylinders for precompression, and the Roman numerals are used to represent the precompression sequence:

• I. Integral positioning and precompression.

• II. Two-wing positioning precompression.

• III.IV. The most severe part of cracks.

• V. VI. The most severe part of wrinkles.

Optimization of solution

According to the genetic algorithm, the solution of precompression direction optimization can be got. The precompression sequence I and II do not need to change the precompression direction, it is still −Z direction, the precompression sequence of the precompression sequence III and IV is (± $45°$ , $0°$ ), and the precompression direction of V is ( $0°$ , $45°$ ), the precompression direction of VI is (± $90°$ , − $45°$ ).

Verification of simulation

Through the Step-by-step method, the precompression direction of some areas are changed, the remaining conditions remain unchanged, and the thickness distribution map of the simulated stepwise precompression is obtained, as shown in Figure 13.

Figure 13.

Thickness distribution of Step-By-Step preforming after optimizing the direction.

From Figure 13, we can find that the preforming effect is better. The thickness of the H part is 1.74–1.92 mm, and the thinning effect is reduced obviously. The thinning rate is 4% −13%, which is a normal thinning situation. The thickness of the G part is between 2.05–2.12 mm, and the wrinkling is also slowed down significantly.

Conclusions

In order to solve the problem of the inaccurate dimensions of preform, the principle of blank expansion is analyzed in this paper, and the finite element inverse algorithm is used to expand the blank of preform according to the constitutive relationship of the thermoset composite prepreg, then a typical composite material is selected. The battery box is taken as an example, and it is realized by ABAQUS programming, and then compared with DYNAFORM, the expansion shape of the inverse finite element algorithm is more accurate, and the amount of prepreg can be reduced by 5.45%. Aiming at the problems of cracks and wrinkles in the preform during the integral precompression, the preform of the battery box is divided into 11 areas in this paper, and the Step-By-Step precompression method is used for preforming. Then, the evaluation function is established according to the precompression direction of different areas, and the precompression direction is optimized by the genetic algorithm, the optimal angle of the precompression direction is obtained in different regions. In the numerical simulation, the maximum thinning rate is reduced to 13%, and the maximum thickening rate is reduced to 6%, the cracks and wrinkles of the preform have been improved significantly.