Mathematical Modeling and Numerical Optimization of Composite Structures

Abstract

This chapter is devoted to modeling the properties of composite materials and structures. Mathematical relations describing the nonlinear elastic three-point bending of isotropic and reinforced beams with account of different strength and stiffness behavior in tension and compression are obtained. An algorithm for numerical solution of corresponding boundary-value problems is proposed and implemented. Results of numerical modeling were compared to acquired data for polymer matrix and structural carbon fiber reinforced plastics. A computational technology for analysis and optimization of composite pressure vessels was developed and presented.

1. Introduction

Carbon fiber reinforced plastics (CFRP) are the most promising modern composite materials. High-duty structures used in aviation and space industry, car manufacturing and building sector require new CFRPs as well as ways to improve their characteristics. Applying computer modeling techniques significantly reduces both the time and cost of investigations aimed at searching optimal parameters of CFRP structures [1]. Mathematical modeling provides an opportunity for comprehensive analysis of both CFRPs and CFRP structures. It has become an effective tool for solving important applied problems.

To build a mathematical model of composite materials, including those made of carbon fibers, one relies upon the experimental data acquired in mechanical testing. Wide application of digital testing machines has brought such experiments to higher level of quality. By measuring a large number of parameters with a high discretization frequency, modern testing machines allow for high amount of information on material deformation and failure to be obtained within a single experiment. Therefore, data processing has become an important step for mathematical modeling of CFRPs and CFRP structures. This is exceptionally important because of quite specific behavior of CFRPs and of their components: fibers and matrices.

One of the features of such materials is their different strength and stiffness behavior in tension and compression combined with nonlinearities of stress-strain curves. Multiple studies for epoxy matrices showed that their ultimate strains in tension were much lower than in compression: approximately 4 versus 20% and more [2]. Moreover, under tension and compression, the deformation behavior of epoxy matrices significantly differs. The corresponding stress-strain curves have different stiffness (secant modulus) at the same values of strain. The similar difference can be observed for CFRPs. In [3, 4, 5], it was shown that in tension tests of carbon fiber specimens with reinforcement angles less than ${20}^{\circ }$ the stiffness grows together with strains (stiffening), whereas for epoxy matrices softening is observed. This phenomenon was explained by the properties of carbon fibers.

Contrasting behavior in tension and compression, stiffening, softening and other nonlinearities are forcing researchers to build and use special mathematical models and computing algorithms. Mathematical models taking into account the abovementioned properties of materials were proposed and studied theoretically by Timoshenko [6] and Ambartsumyan [7, 8] in the mid-twentieth century. Later, Jones had experimentally, theoretically and numerically studied the nonlinear behavior of several fiber-reinforced composites. The main focus of the research was on the difference in stiffness and strength behavior under tension and compression [9]. After Ambartsumyan’s and Jones’ researches, a lot of studies were dedicated to this problem. Most of them were dealing with linear bi-modulus models of materials or 3D finite elements. In [10], Ambartsumyan with a coauthor suggested a theoretical approach to modeling of multimodulus nonlinear elastic beams under bending, but still without calculations.

Another trend is studying the behavior of sandwich panels or beams with a CFRP faces having differences in tension and compression along with the other mentioned nonlinearities [11, 12, 13]. These works concern the problem of flexure of CFRPs and similar materials. They consider bending of specimens as a reference test. The first paper [11] is devoted to experimental investigations and shows most of the nonlinearities we supposed such materials should have: stiffening in tension, softening in compression, different moduli even at the origin of coordinates. Other two works [12, 13] present more complex studies including full cycle of mathematical modeling spanning from the experimental investigations to numerical ones.

A comprehensive approach to modeling and simulation of nonlinear elastic deformation of polymer matrices and different CFRPs was presented in [14]. This chapter deals with different strength and stiffness behavior of the materials in tension and compression exemplified by a case of three-point bending. This approach implements a full cycle of model development and validation, which comprises the following stages: carrying out tests and acquiring experimental data, data prepossessing and building stress-strain curves, analytical approximation of acquired curves, mathematical modeling and numerical simulation of deformation processes, comparative analysis of results of numerical modeling to acquired data.

2. Structural models of composite materials

For most of the composite materials models, we can write the relations between average stresses ${\sigma }_{\mathit{\alpha \beta }}$, ${\tau }_{\alpha 3}$ and strains $e_{\mathit{\alpha \beta }}$, ${\gamma }_{\alpha 3}$ (generalized Hooke’s Law):

$$\begin{array}{l}{\sigma }_{\mathit{\alpha \alpha }}=a_{\mathit{\alpha \alpha }}{e}_{\mathit{\alpha \alpha }}+a_{\mathit{\alpha \beta }}{e}_{\mathit{\beta \beta }}+a_{\alpha 3}\cdot 2e_{\mathit{\alpha \beta }}-a_{\alpha \mathrm{\Theta }}\mathrm{\Theta },\\ {\sigma }_{\mathit{\alpha \beta }}=a_{\alpha 3}{e}_{\mathit{\alpha \alpha }}+a_{\beta 3}{e}_{\mathit{\beta \beta }}+a_{33}\cdot 2e_{\mathit{\alpha \beta }}-a_{3\mathrm{\Theta }}\mathrm{\Theta },\\ {\gamma }_{\alpha 3}=q_{\mathit{\alpha \alpha }}{\tau }_{\alpha 3}+q_{\mathit{\alpha \beta }}{\tau }_{\beta 3},\end{array}$$E1

where $\mathrm{\Theta }$ is the increase of temperature. Relations (Eq. (1)) are called the thermoelasticity relations, or, when no temperature influence is considered, they are simply elasticity relations.

The structural model of fiber reinforced composite described in [15, 16, 17, 18] has become a foundation for a large number of current researches. Now it is widely used while simulating the behavior of composite structures. The model is based on the following assumptions: the stress-strain state into isotropic elastic fibers and into entire volume of isotropic ideally elastic matrix is homogeneous; fibers and matrix are deformed jointly along the direction of reinforcement; stresses in fibers and in matrix corresponding to other directions are equal.

For computing the effective elastic modulus of unidirectional fiber-reinforced composite, the Reuss-Voigt average was used giving the following formulae

$$\begin{array}{l}{E}_1={\omega }_f{E}_f+{\omega }_m{E}_m,E_2=\frac{E_f{E}_m}{{\omega }_f{E}_m+{\omega }_m{E}_f},\\ {\nu }_{12}={\omega }_f{\nu }_f+{\omega }_m{\nu }_m,G=\frac{G_f{G}_m}{{\omega }_f{G}_m+{\omega }_m{G}_f},\end{array}$$E2

where all the terms having squared Poisson coefficients are neglected.

Herewith $E_1$, $E_2$ are effective moduli along and across the direction of reinforcement, $G$ is effective share modulus, ${\nu }_{12}$ is effective Poisson coefficient in the plain of layer; $E$, $\nu$, $\omega$ with “$f$” and “$m$” indices are elastic moduli, Poisson coefficients and volume fractions of matrix and fibers correspondingly, hereby ${\omega }_m+{\omega }_f=1$.

On the ground of symmetry of compliance tensor, one has

$${\nu }_{21}={\nu }_{12}{E}_2{E}_1^{-1}.$$

Formulae for effective coefficients of thermal expansion have the following form

$${\alpha }_1={\omega }_f{\alpha }_f+{\omega }_m{\alpha }_m,{\alpha }_2=\frac{{\omega }_f{\alpha }_f{E}_f+{\omega }_m{\alpha }_m{E}_m}{{\omega }_f{E}_f+{\omega }_m{E}_m}.$$E3

In description of the model, it is noted that among formulae for effective moduli, those obtained using Reuss averaging (in particular formulae for $G$) lead to the worst results. Estimations for $G$ obtained using variational method are also obtained, and it is shown that lower boundary

$$G=\frac{\left(1+{\omega }_f\right)G_f+{\omega }_m{G}_m}{\left(1+{\omega }_f\right)G_m+{\omega }_m{G}_f}{G}_m$$E4

gives more accurate approximation than (Eq. (2)) does. Hereinafter share moduli of matrix and fibers are

$$G_m=\frac{E_m}{2\left(1+{\nu }_m\right)},G_f=\frac{E_f}{2\left(1+{\nu }_f\right)}.$$

Components of effective stiffness tensor for unidirectionally reinforced layer in case of state of plane stress have the following form:

$$A_{\mathit{\text{αααα}}}=\frac{E_{\alpha }}{1-{\nu }_{12}{\nu }_{21}},A_{1122}=\frac{{\nu }_{21}{E}_1}{1-{\nu }_{12}{\nu }_{21}},A_{1212}=G.$$E5

Unwritten expressions can be obtained using symmetry rule or vanish. Hereinafter, we assume $\alpha ,\beta =1,2$ and $\alpha \ne \beta$.

The coefficients in the relations (Eq. (1)) for example are defined by the formulas given in [1, 17, 18].

3. Mathematical model and numerical analysis of reinforced beam deformation

Three-point bending flexural test has been one of the standard techniques to determine physical and mechanical characteristics of materials. Figure 1 shows a scheme of physical model of three-point bending of a beam with the rectangular cross section $b\times 2h$, and the span $l$ between the supports. The left edge of the beam is hinged, while the right one is supported freely. The force $P$ is applied to the center of the beam. The model neglects the shape of the supports and assumes the occurring load $P$ and support reactions $R_A$ and $R_B$ to be concentrated. In addition, the model neglects the possible heterogeneity of the deformations in the direction normal both to the longitudinal direction and to the load direction.

Figure 1.

Three-point bending of a rectangular-sectioned beam.

In this case, the beam’s upper part undergoes compression strain in the longitudinal direction, bottom part—tension strain. VSE-1212 polymer matrix and VKU-28 (T-800 carbon yarn plus VSE-1212 epoxy matrix) structured CFRP react differently to tension and compression. VKU-28 has been one of the most promising types of CFRPs that is going to be used in the latest generations of aircrafts. The effect of accounting for this factor on the computational results is essential. Further, these results are compared to acquired ones.

Due to very low deformation rates, the classical theory of beam bending can be regarded as satisfactory for description of the equilibrium state. To this end, it is convenient to consider the beam’s median surface as a reference one.

The beam’s stress-strain state is characterized by the following values determined on the reference surface: the shear force $Q\left(x\right)$, the bending moment $M\left(x\right)$, the longitudinal force $N\left(x\right)$, and by the longitudinal displacement and bend ($u\left(x\right)$, $w\left(x\right)$ respectively). The corresponding equilibrium equations are written as follows:

$$\frac{\mathit{dN}}{\mathit{dx}}=0,\frac{\mathit{dQ}}{\mathit{dx}}=0,\frac{\mathit{dM}}{\mathit{dx}}=Q.$$E6

The reactions $R_A$ and $R_B$ can be determined by considering force equilibrium $R_A=R_B=P/2$. The bending moments at the support points are equal to zero: $M_A=M_B=0$. The solution of the equation system (Eq. (6)) can be expressed as follows:

$$\begin{array}{l}N=0,Q\left(x\right)=\left\{\begin{array}{ll}P/2,& 0\le x\le l/2,\\ -P/2,& l/2\le x\le l,\end{array}\right.\\ M\left(x\right)=\left\{\begin{array}{ll}\mathit{Px}/2,& 0\le x\le l/2,\\ -P\left(x-l\right)/2,& l/2\le x\le l.\end{array}\right.\end{array}$$E7

Strain distribution for the beam’s thickness can be obtained from the Kirchhoff-Love kinematic hypotheses:

$$\epsilon \left(x,z\right)=e\left(x\right)+\mathit{z\kappa }\left(x\right),$$E8

$$e\left(x\right)=\frac{\mathit{du}}{\mathit{dx}},\kappa \left(x\right)=-\frac{d^{2}w}{{\mathit{dx}}^2},$$E9

where $\epsilon \left(x,z\right)$ is the strain in the beam; $e\left(x\right)$ is the median surface strain; and $\kappa \left(x\right)$ denotes changes in the median surface curvature. As mentioned earlier, the beam undergoes tension and compression strain, whose interface will be marked as $z_1$. In this case, for the section area $-h\le z\le z_1$, the strain will be negative, and for $z_1\le z\le h$ positive. At the interface of these two states, the strains $\epsilon$ vanish, so the interface itself is determined as follows:

$$z_1=-\frac{e}{\kappa },-h\le z_1\le h.$$E10

The constitutive equation can be expressed as:

$${\sigma }^{\pm }\left(x,z\right)=f_i^{\pm }\left(\epsilon \right),$$E11

where the superscript “+” refers to the areas with positive strains and “–” – to the area with negative ones; $f\left(\epsilon \right)$ denotes the approximation selected for the stress-strain curve (a linear function, a polynomial, or a combination of linear and power-law functions).

The longitudinal force $N$ and the bending moment $M$ in the beam cross section are determined by the equations:

$$\begin{array}{ll}& N=b\left(\underset{-h}{\overset{z_1}{\int }}{\sigma }^{-}\mathit{dz}+\underset{z_1}{\overset{h}{\int }}{\sigma }^{+}\mathit{dz}\right),\\ & M=b\left(\underset{-h}{\overset{z_1}{\int }}{\sigma }^{-}\mathit{zdz}+\underset{z_1}{\overset{h}{\int }}{\sigma }^{+}\mathit{zdz}\right).\end{array}$$E12

Having substituted (Eq. (12)) with the relations (Eq. (8)), (Eq. (10)), (Eq. (11)) into and integrated it over the beam thickness, one obtains a system of equations to determine $\kappa$ and $e$:

at $0\le x\le l/2$

$$\left\{\begin{array}{l}N\left(\kappa ,e,x\right)=0,\\ M\left(\kappa ,e,x\right)=\mathit{Px}/2,\end{array}\right.$$E13

at $l/2<x\le l$

$$\left\{\begin{array}{l}N\left(\kappa ,e,x\right)=0,\\ M\left(\kappa ,e,x\right)=-P\left(x-l\right)/2.\end{array}\right.$$E14

The system of equations (Eq. (13)) and (Eq. (14)) in general case is nonlinear, but in the case of piecewise linear constitutive equations which take into account different strength and stiffness behavior in tension and compression expressed as follows:

$${\sigma }^{\pm }\left(x,z\right)=E^{\pm }\epsilon ,$$E15

it can be solved analytically. In the nonlinear case, the Newton method is applied to solve the equations (Eq. (13)) and (Eq. (14)), and then, the linearized system

$$\begin{array}{l}N\left({\epsilon }_0,{\kappa }_0\right)+\frac{\partial N\left({\epsilon }_0,{\kappa }_0\right)}{\partial \epsilon }\left(\epsilon -{\epsilon }_0\right)+\frac{\partial N\left({\epsilon }_0,{\kappa }_0\right)}{\partial \kappa }\left(\kappa -{\kappa }_0\right)=0,\\ M\left({\epsilon }_0,{\kappa }_0\right)+\frac{\partial M\left({\epsilon }_0,{\kappa }_0\right)}{\partial \epsilon }\left(\epsilon -{\epsilon }_0\right)+\frac{\partial M\left({\epsilon }_0,{\kappa }_0\right)}{\partial \kappa }\left(\kappa -{\kappa }_0\right)=M\left(x\right)\end{array}$$

can be solved for unknown values

$$\kappa =F\left({\epsilon }_0,{\kappa }_0,M\left(x\right)\right),\epsilon =G\left({\epsilon }_0,{\kappa }_0\right),$$E16

where ${\epsilon }_0$ and ${\kappa }_0$ are the initial approximations, and $M\left(x\right)$ is determined from (Eq. (7)).

As the initial approximation at small values of the load $P$, the solutions obtained for the linear constitutive equations (Eq. (15)) were used. Since the computation is performed with a relatively small increment of $P$, in case of big values of $P$, one can use the computation results acquired at a previous step as the initial approximation for current step.

Having determined the change of median surface curvature from the equations (Eq. (13)) and (Eq. (14))

$$\kappa \left(x\right)=\left\{\begin{array}{ll}{\kappa }_1\left(x\right),& \mathrm{at}x\in \left[0,l/2\right),\\ {\kappa }_2\left(x\right),& \mathrm{at}x\in \left[l/2,l\right],\end{array}\right.$$

one can write down a differential equation to determine the beam bend. For that purpose, the bend function is expressed as follows:

$$w\left(x\right)=\left\{\begin{array}{ll}{w}_1\left(x\right),& \mathrm{at}x\in \left[0,l/2\right),\\ w_2\left(x\right),& \mathrm{at}x\in \left[l/2,l\right].\end{array}\right.$$

Using the equation (Eq. (9)) and the beam’s fixing conditions, a system of equations can be derived:

$$\begin{array}{l}\frac{d^2{w}_1}{{\mathit{dx}}^2}=-{\kappa }_1,\frac{d^2{w}_2}{{\mathit{dx}}^2}=-{\kappa }_2,\\ w_1\left(0\right)=w_2\left(l\right)=0,w_1\left(l/2\right)=w_2\left(l/2\right),\\ \frac{{\mathit{dw}}_1\left(l/2\right)}{\mathit{dx}}=\frac{{\mathit{dw}}_2\left(l/2\right)}{\mathit{dx}}.\end{array}$$

The solution of these equations can be obtained using the methods of solving boundary-value problems for systems of ordinary differential equations. For that purpose, the modified collocation and least-residuals method [19, 20, 21] were applied.

Numerical analysis of deformation processes in VSE-1212 polymer matrix and VKU-28 structured carbon fiber is based on approximation of stress-strain curves and the three-point bending model.

Further, three different specimens with the geometrical sizes $l\times 2h\times b$ are considered:

1. specimen 1—VSE-1212 polymer matrix, $75\times 4.78\times 10.05\mathrm{mm}$,

2. specimen 2—VKU-28 structured carbon fiber (the specimen was cut out along the reinforcement), $90\times 3.45\times 9.85\mathrm{mm}$;

3. specimen 3—VKU-28 structured carbon fiber (the specimen was cut out perpendicular to the reinforcement), $90\times 3.40\times 9.95\mathrm{mm}$.

In Figure 2, one can see the simulation results for beam three-point bending, obtained through different approaches to approximation of the constitutive equations, and their comparison with the experimental data.

Figure 2.

Experimental (solid curves) and dependencies of beam-deflection and load obtained in simulation: linear approximation (1); quadratic approximation by a polynomial of the second degree (2); cubic approximation (3); linear and power-law approximation (4); a–c are specimens 1 − 3 respectively; d—the solution to a three-point bending problem without account for the different strength and stiffness behavior in tension (curve 1) and compression (curve 2). The solid line shows the results of mechanical tests.

Applying the linear dependencies to tension and compression has not resulted in adequate approximation even for 30% of curve. Using more complex than quadratic approximation laws at first led to a significant deviation from the experimental curve and then to divergence of the Newton method iteration process. This is explained by the fact that in tension tests, due to specimens’ fragility, the strain range for the polymer matrix specimens was limited to 2%, while in the bending tests, the strains in tension zone reached 4–5%.

Thus, to solve the bending problem, the tension curve was extrapolated into the domain of high strains. The extrapolations obtained using a polynomial of the third degree, and by linear and power-law function reached the maximum too quickly and then started to decrease, which is against the physics behind the deformation process. A similar effect was observed when calculating the bending of the carbon-fiber specimens cut out along direction of reinforcement filler.

The calculations using quadratic approximation and extrapolation of tension curves and approximation of compression curves within a short (up to 6%) segment have turned out to be best for qualitative and quantitative description of the nonlinear character of VSE-1212 polymer matrix bending. In the case of the specimen cut out perpendicular to direction of its reinforcement, all the approximations have shown the results close to experiment. At the same time for the test with the maximum load, the best option has still been application of quadratic approximations.

Taking different strength and stiffness behavior in tension and compression into account has an essential effect. As it was demonstrated earlier, the tension tests of VKU-28 specimens produced nonlinear stress-strain curves, while the difference of characteristics between tension and compression reached 5–7% for the longitudinal reinforcements and 12–15%—for the transverse ones (see Tables 1 and 2).

Approximation type	Approximation coefficients			MSD
	$a_1\cdot {10}^{-9}$	$a_2\cdot {10}^{-9}$	$a_3$
Tension of VKU-28 CFRP, $\epsilon \in \left[0;0.015\right]$
A1	160.8			1.4e – 2
A2	144.9	1.44e + 3		5.7e – 4
A3	144.0	1.66e + 3	–1.14e + 13	4.3e – 4
A4	143.0	8.87e + 2	1.87	4.1e – 4
Compression of VKU-28 CFRP, $\epsilon \in \left[0;0.0018\right]$
A1	155.4			5.8e – 3
A2	160.2	–3.33e + 3		3.4e – 3
A3	155.9	4.31e + 3	−3.00e + 15	2.9e – 3
A4	157.7	−5.81e + 8	3.98	3.0e – 3

Table 1.

Approximation coefficients for stress-strain curves of VKU-28 carbon fiber specimens reinforced in longitudinal direction and mean square deviation (MSD) of $f\left(x\right)$.

Approximation type	Approximation coefficients			MSD
	$a_1\cdot {10}^{-9}$	$a_2\cdot {10}^{-9}$	$a_3$
Tension of VKU-28 CFRP, $\epsilon \in \left[0;0.0076\right]$
A1	7.37			1.1e – 2
A2	7.89	−9.22e + 1		6.4e – 4
A3	7.87	−8.97e + 1	2.54e + 11	4.3e – 4
A4	7.82	−2.23e + 2	2.20	4.8e – 4
Compression of VKU-28 CFRP, $\epsilon \in \left[0;0.0034\right]$
A1	8.90			6.7e – 3
A2	9.21	−1.20e + 2		4.1e – 3
A3	8.96	1.16e + 2	−5.09e + 13	3.7e – 3
A4	9.04	−5.16e + 6	3.95	3.7e – 3

Table 2.

Approximation coefficients for stress-strain curves of carbon fiber specimens reinforced in transverse direction and mean square deviation (MSD) of $f\left(x\right)$.

For the polymer matrix this difference exceeded 15% (see Tables 3 and 4).

Approximation type	Approximation coefficients			MSD
	$a_1\cdot {10}^{-9}$	$a_2\cdot {10}^{-9}$	$a_3$
Strain of VSE-1212 polymer matrix (Constant cross section) $\epsilon \in \left[0;0.018\right]$
A1	3.30	—	—	2.9e − 2
A2	3.90	−4.38e + 1	—	1.5e − 3
A3	3.83	−3.17e + 1	−4.94e + 11	6.7e − 4
A4	3.80	−1.05e + 2	2.25	6.6e − 4
Strain of VSE-1212 polymer matrix (Variable cross section)$\epsilon \in \left[0;0.018\right]$
A1	3.33			2.7e − 2
A2	3.89	−4.02e + 1		1.8e − 3
A3	3.80	−2.48e + 1	−6.30e + 2	4.0e − 4
A4	3.77	−1.40e + 2	2.35	2.7e − 4

Table 3.

Approximation coefficients for tension curves of VSE-1212 polymer matrix, $\epsilon \in \left[0;0.018\right]$ and mean square deviation (MSD) of $f\left(x\right)$.

Approximation type	Approximation coefficients			MSD
	$a_1\cdot {10}^{-9}$	$a_2\cdot {10}^{-9}$	$a_3$
Compression of VSE-1212 polymer matrix (Constant cross section), $\epsilon \in \left[0;0.28\right]$
A1	0.77			2.3e − 1
A2	1.60	−3.97		8.2e − 2
A3	2.36	−1.29e + 1	2.37e + 10	1.4e − 2
A4	−5.71	5.49	0.90	3.8e − 2
Compression of VSE-1212 polymer matrix (Variable cross section), $\epsilon \in \left[0;0.28\right]$
A1	0.69			3.5e − 1
A2	1.69	−5.22		1.4e − 1
A3	2.71	−1.84e + 1	3.84e + 1	3.6e − 2
A4	−2.07	1.72	0.72	4.8e − 2
Compression of VSE-1212 polymer matrix $\epsilon \in \left[0;0.06\right]$ (Variable cross section, shortened test area)
A1	2.10			7.2e − 2
A2	3.05	−2.12e + 1		4.2e − 3
A3	3.18	−2.85e + 1	9.13e + 1	1.1e − 3
A4	3.31	−1.24e + 1	1.75	1.9e − 3

Table 4.

Approximation coefficients for compression curves of VSE-1212 polymer matrix and mean square deviation (MSD) of $f\left(x\right)$.

Tables 1–4 show approximation results for the above-presented stress-strain curves by different functions at different intervals:

1. by the linear approximation $\sigma =a_1\epsilon$ (A1),

2. by the polynomial of the second degree $\sigma =a_1\epsilon +a_2{\epsilon }^2$ (A2),

3. by the polynomial of the third degree $\sigma =a_1\epsilon +a_2{\epsilon }^2+a_3{\epsilon }^3$ (A3),

4. by a combination of linear and power-law functions $\sigma =a_1\epsilon +a_2{\epsilon }^{a_3}$ (A4).

However, if bending tests have been performed to determine an elasticity module of CFRPs, different strength and stiffness behavior in tension and compression is compensated and one obtains some averaged characteristic.

It is useful to consider the effect of the way for determining and setting of the mechanical characteristics on modeling of three-point bending of the carbon-fiber beam cutout perpendicular to its reinforcements. Figure 2d shows the solutions obtained while using a linear approximation of the constitutive equations with equal elastic moduli for tension and compression: for curve 1 the modulus was obtained from tension experiments, for curve 2—from compression ones (see Table 2).

As one can see the calculated linear results without account for the different strength and stiffness behavior in tension and compression have differed from the results of mechanical tests (the solid curve) by more than 15%.

Most of the real CFRP structures under day-to-day service conditions bear complex loads that result in formation of tension, compression and bending zones as well as their combinations in the structures. Applying the traditional methods for determination of material characteristics in combination with linear deformation models (in particular those that do not account for the different strength and stiffness behavior in tension and compression) for calculation of such structures, one risks to distort the deformation and stress pattern significantly, which, in its turn, results in either underestimation or overestimation of the structure’s strength and rigidity. Keeping in mind that carbon fibers are used for manufacturing of high-duty structures, their computation demands different strength and stiffness behavior in tension and compression to be taken into account.

4. Numerical analysis and design of pressure vessels

Composite overwrapped pressure vessels (COPV) are used in the rocket and spacecraft making industry due to their high strength and lightweight. Consisting of a thin, nonstructural liner wrapped with a structural fiber composite COPV are produced to hold the inner pressure of tens and hundreds atmospheres. COPV have been one of the most actual and perspective directions of research, supported especially by NASA [22, 23].

Designing of a highly reliable and efficient COPV requires a technology for analysis of its deformation behavior and strength assessment. This technology should allow one to obtain target COPV parameters through changing vessel’s geometry, structural and mechanical material parameters while keeping its useful load.

Application of combination mathematical modeling and numerical optimization makes it possible to reduce the cost and the duration of identifying the best parameters for a COPV. However, this approach is characterized by a number of hurdles. Overcoming these hurdles determines the success of an optimum designing of such structures.

So far, there have been two main approaches in optimization of composite structures: analytical and numerical ones.

In the first approach, the problems are solved basing on their simplified statement, for example using the momentless (membrane) shell theory and the netting model of composite material (CM) [24, 25, 26, 27]. The obtained results may be far from reality; however, they are of value for testing of numerical optimization methods.

Application of the numerical approach in designing, on the other hand, produces a number of challenges that must be overcome, for example, lack of reliable methods for global optimization; nonconvexity and nonlinearity of constraint functions; ill-conditioned boundary value problems; different scaling of optimization criteria represents just some of the obstacles that prevent from reliable optimization of COPV.

Numerical analysis is usually a computation-intensive process and takes considerable time. One way to solve this problem is approximation of the objective function using different approaches, such as response surface method [28] and neural network [29]. Some kinds of numerical analyses use a small number of design variables, functions and/or corresponding set of their discrete values (analytical geometry parametrization [30], finite set of feasible winding angles [31]).

Another way is reasonable simplification of the elasticity problem statement, for example by using the membrane theory or other shell theories [30, 32, 33], that leaves the question of results validity. This is the approach we have applied in our study. For validation, we have used the Timoshenko [34] and Andreev-Nemirovskii [35] shell theories, accounting transverse shears with different degrees of accuracy.

Of course, it should be taken into account that the computed solutions are not optimum in the strict mathematical sense. However, these solutions could provide the considerable economy of the weight while keeping the required strength, and, therefore, they have high engineering value.

4.1. The problem statement and the mathematical models

Let us consider a multilayer composite pressure vessel at a state of equilibrium under equidistributed inner pressure. We need to determine the parameters of structure and CM meeting the following requirements:

$$V\ge V_0,P\ge P_0,M\le M_0,$$E17

where $V$ is the volume of the vessel, $P$ is inner pressure and $M$ is the vessel’s mass and they are constrained by some preset values $V_0,P_0,M_0$.

We define the optimization problems the following way: to find extremum of one functional from (Eq. (17)) under other constraints.

The structural optimization problem statement includes selection of objective functional, formulation of constitutive equations and constraints on performance and design variables.

The mathematical models describing the vessel’s state are based on the following assumptions:

1. the vessel is a multilayer thin-walled structure;

2. the vessel’s layers can have different mechanical characteristics;

3. the reinforced layer’s material is quasi-homogeneous;

4. the vessel’s main loading is high inner pressure.

These assumptions allow us to reduce dimension of the corresponding mathematical problem and to build the mathematical vessel’s models based on the different theories of multilayer nonisotropic shells.

Let us consider the vessel as a shell rigidly compressed on the edge. Taking into account a symmetry plane in the middle of the vessel, it is enough to calculate and design only its one half. The type of loading and boundary conditions allows considering the axisymmetric problem statement.

The shell is set by rotation of the generatrix $r=r\left(\theta \right)$ around axis $0y$ (Figure 3) where $r$ is the current point of the shell radius, $\theta$ is the angle between the normal to the shell surface and the spin axis changing within $\left[{\theta }_0;{\theta }_1\right]$.

Figure 3.

Shell of rotation geometry.

The Kirchhoff-Love shell theory [36] (KLST) and the improved Timoshenko [34] (TiST) and Andreev-Nemirovskii [35] (ANST) theories are used to solve the direct calculation problems of multilayer composite vessels, to analyze their behavior and to verify optimization problem solutions. The full systems of equations were described in the paper [17].

Relations between stresses and strains are described by the structural models [18]. The main idea of these models is that CM parameters are calculated through matrix and fibers mechanical parameters, fibers volume content and winding angles. The stress-strain state of matrix and fibers is evaluated through stresses and strains of the composite shell. A failure criterion is applied for every component of CM. Here we use the Mises criterion to determine the first stage of failure.

The objective function whose minimum is required is the minimum mass:

$$M=2\pi {\int }_{{\theta }_0}^{{\theta }_1}{\mathit{rR}}_1\mathit{hd\theta }\left[{\rho }_m\left(1-{\omega }_r\right)+{\rho }_r{\omega }_r\right]\to \text{min},$$E18

where ${\rho }_m,{\rho }_r$ are the densities of matrix and reinforcing fibers, ${\omega }_r$ is the volume content of reinforcement.

We chose the following design functions: the curvature radius $R_1\left(\theta \right)$ to define the generatrix; the thickness of the shell $h\left(\theta \right)$; the reinforcement angle $\psi \left(\theta \right)$ (Figure 3).

The solution has to satisfy the constraints on the shell’s inner volume:

$$\pi {\int }_{{\theta }_0}^{{\theta }_1}{r}^2{R}_{1}sin\mathit{\theta d\theta }=V_0$$E19

and the strength requirement:

$$\text{max}\left\{{\mathit{bs}}_r,{\mathit{bs}}_m\right\}\le 1,$$E20

where ${\mathit{bs}}_r,$ ${\mathit{bs}}_m$ are the normalized von Mises stresses in the matrix and fibers [1]. Note that the factor of safety is widely used while solving engineering problems. It can be considered by correction of the right-hand side of the inequality (Eq.(20)).

We used the following constraints on the design functions:

$$0\le \psi \le 90,h_0^{\ast }\le h\le h_1^{\ast },R_0^{\ast }\le R_1\le R_1^{\ast }.$$E21

The method of the continuous geodesic winding has been widely used in the manufacturing of composite shells of revolutions. In this case the winding angles are defined by the Clairaut’s formula:

$$rsin\psi \left(r\right)=C,$$E22

where $C$—the constant is defined, as a rule, from the condition at the equator of the shell. The thickness equation is

$$h\left(r\right)=h_R\frac{Rcos{\psi }_R}{rcos\psi \left(r\right)},$$E23

which has the singularity at the edge where the winding angle has to be equal to ${90}^{\circ }$. The formula (Eq. (23)) is applied into practice at $r\ge r_0+r_{\omega }$, where $r_{\omega }$ is equal to the width of the reinforcement tape. As a result, the equation determining the vessel’s thickness takes the form:

$$h\left(r\right)=\left\{\begin{array}{ll}{h}_R\frac{Rcos{\psi }_R}{r_{\omega }cos\psi \left(r_0+r_{\omega }\right)},& r\le r_0+r_{\omega };\\ h_R\frac{Rcos{\psi }_R}{rcos\psi \left(r\right)},& r\ge r_0+r_{\omega }.\end{array}\right.$$E24

We did not consider the problem of fibers slippage. The main goal of the study was to demonstrate the potentials of using CM.

4.2. Direct problems: analysis of the shell theories

Estimation of composite vessel stress-strain state using offered models leads to the solution of boundary value problems for rigid systems of differential equations. These problems are ill-conditioned, and their solutions have pronounced character of thin boundary layers. Numerical analysis was performed by the spline collocation and discrete orthogonalization methods, realized in the COLSYS [37] and GMDO [38] software. These computing tools have proved to be effective in numerical solving of wide range of problems of composite shell mechanics [1].

We investigated the vessel’s deformations by computing its stress-strain state based on the different shell theories. The vessel’s shape was a part of a toroid: $R_1=2.46\mathrm{m}$, ${\theta }_0={0.108}^{\circ }$, ${\theta }_1={90}^{\circ }$ (the computed half), $r\left({\theta }_0\right)=0.04\mathrm{m}$. The carbon composite parameters were: $E_m=3\cdot {10}^9\mathrm{Pa}$, ${\nu }_m=0.34$, $E_r=300\cdot {10}^9\mathrm{Pa}$, ${\nu }_r=0.3$, ${\omega }_r=0.55$, $V_0=350$ liters where $E_m,E_r$ are the Young’s modulus of the matrix and fibers, ${\nu }_m,{\nu }_r$—their Poisson’s ratio.

Figure 4 shows the stress-strain state characteristics of the vessel with the thickness $h=0.6\mathrm{cm}$, reinforced in the circumferential direction ($\psi ={90}^{\circ }$) under the load of 170 atm. On the left, the displacements of the reference surface along the generatrix $u_1\left(r\right)$ (dashed curves) and the normal displacement of these surfaces $w\left(r\right)$ (solid curves) are shown. On the right is the distribution of normalized von Mises stress (nVMS) along the thickness in the matrix ${\mathit{bs}}_m\left(r\right)$. The solid curves correspond to a slice at the shell edge, the dashed curves — to a slice at $\theta =0.1$.

Figure 4.

The stress-strain state characteristics of the composite vessel computed using different shell theories. Longitudinal displacement $u$—dashed curves; deflection $w$—solid curves. The curves without symbols correspond to KLST simulations, the curves marked with $\Delta$—to those using TiST, and $\square$—to ANST.

It is easy to see that the basic kinematic characteristics coincide both qualitatively and quantitatively. Small differences are observed only for the stresses and deformations near the compressed edge. The maximum results and qualitative difference were obtained for ANST. This is due to accounting for the transverse shears by nonlinear distribution in a thickness of a shell. Earlier it was shown [1] that ANST’s-based results were the closest to the ones of 3D elastic theory in most cases.

The winding angle’s influence on the COPV performance was investigated using parametric analysis. Dependence of the maximum nVMS in the matrix ${\mathit{bs}}_m$ (dashed curves) and the fibers ${\mathit{bs}}_r$ (dash-dotted curves), and the maximum size of the displacement vector $\Vert \overrightarrow{v}\Vert$ (solid curves) are shown in Figure 5.

Figure 5.

The winding angle’s influence on the composite vessel stress-strain state. KLST’s results are drawn without marks, TiST — with symbols $\Delta$, ANST — with $\square$.

The calculated values are very close in the area of their minima (Figure 5 left side). The graphs of kinematic function $\Vert v\Vert$ coincide qualitatively. Some noticeable quantitative differences are revealed only for KLST’s results.

The range $\psi \in \left(42;45\right)$ corresponds to the zones of minimum values (Figure 5 right side), which practically coincide (${\text{min}}_{\psi }{\mathit{bs}}_m\approx 0.65$, ${\text{min}}_{\psi }{\mathit{bs}}_r\approx 1.05$, ${\text{min}}_{\psi }\Vert v\Vert \approx 5\cdot {10}^{-3}\mathrm{m}$), as well as the angles, where these values are obtained ($\psi \approx {43.2}^{\circ }$ for ${\mathit{bs}}_m$ and ${\mathit{bs}}_r$, $\psi \approx {43.8}^{\circ }$ for $\Vert v\Vert$).

It was revealed that the winding angles of minimum stresses values were almost insensitive to the thickness variation. The change of $h$ from 0.6 to 1.6 cm corresponded to the angle’s change about ${0.2}^{\circ }$.

Additionally, we investigated stress–strain state of the vessel (the thickness $h=0.6\mathrm{cm}$, the winding angles at $\psi =\pm 43.2$), when nVMS in the matrix and the fibers were near their minimum (Figure 6). The adopted notation is the same as in Figure 4.

Figure 6.

The stress-strain state of the vessel ($\psi =\pm 43.2$), computed using the three shell theories.

Again the difference is visible only in a very small region near the edge, but now this difference is small enough to be neglected. Moreover, the displacement values of the reference surface, the efforts and the moments completely coincide for all the theories.

All the theories (KLST, TiST, ANST) provided similar estimated characteristics of stress-strain state. This vessel was characterized not only by essential decrease of the maximal nVMS in the matrix and the fibers, but also by their uniform distribution along the generatrix. At the same time, the values of bending moments significantly reduced bringing vessel’s stress-strain state close to momentless.

The performed analysis showed that the optimization problem can be solved using rather simple shell theories (KLST, TiST). These theories are characterized by lower computational complexity of corresponding boundary value problem if compared to ANST. It takes from 10 to 20 times less resources.

One can see that the winding angle as a design parameter gives an opportunity to increase the vessel’s strength significantly. The difference between the “best” and “worst” designs can reach 20–35 times comparing their nVMS in the matrix and fibers. The “worst” designs have the winding angle close to ${90}^{\circ }$. In this case are considerable transverse shears near the compressed edge, and the loading is redistributed to a rather weak matrix while the fibers remain unloaded.

4.3. Inverse problems: optimization of the vessel

Inverse problems involve not only numerical methods for fast and reliable solving of direct boundary value problems, but also require numerical optimization methods for finding design parameters.

Here we considered conditional optimization problem, including direct constraints on design functions and trajectory constraints on the solution imposed at the end of the interval. The sequential unconstrained optimization is one of the most widespread approaches to solution of such problems. The main idea of the method is terminal functional convolution and multiple solutions of one-criterion problem using different optimization methods [39]. In our study, the modified Lagrange function was used for the convolution.

Hence we sought for solution of a nonconvex problem of finite-dimensional optimization [40] by discretization of design functions. The methods realized in the OPTCON-A software [41] were used to get the corresponding solution.

The considered design with the continuous geodesic winding has been one of COPV widely used in practice [42, 43].

Important additional design characteristic is its “adaptability in manufacturing.” For example, the 5–10 times difference of thickness along the meridian would become a serious obstacle for vessels manufacturing. Thus, designs of nearly minimum mass possessing good properties and satisfying to the given technological constraints could be of great value than optimum without them.

According to Amelina et al. [44], the design with the geodesic continuous winding has the thickness ratio about 10 and large gradient near the edge.

We verified the solutions of optimization problem by substituting the obtained design parameters into the direct problem. In [44] shown that all three theories yielded close results (Figure 7). The difference is noticeable only for ANST in narrow zones (less than $1\%$ of all area of calculation) at the edges, where non-linear accounting for transverse shear gives difference of about $5\%$. At the same time, the estimated efforts and bending moments are very close for all the theories, and the bending moments are very small.

Figure 7.

The stress-strain state characteristics of the vessel with the optimized design functions based on the three shell theories. Longitudinal force $T_{11}$—dashed curves; bending moment $M_{11}$—solid curves.

Thus, it is possible to use the simplest shell theory to solve such optimization problem and the estimation of stress-strain state will be close to those obtained using more complex theories.

5. Conclusions

• Mathematical models for nonlinear flexural deformation of CFRPs and polymer matrices with account for their different strength and stiffness behavior in tension and compression have been built. A satisfactory match with the results of mechanical tests has been obtained. The study has proved that the nonlinear properties of polymer matrices and carbon fibers should be taken into account when calculating and designing real structures.

• The technology of optimization of COPV has been developed. It makes possible to obtain high pressure vessel designs that not only meet such requirements as minimum mass, preset volume and strength, but also possess a number of additional valuable engineering characteristics including stress-strain state close to momentless and almost equally stressed fibers.

• Nonconstant design parameters, such as thickness, winding angles and curvature radius of composite shell give the possibility for additional reduction of COPV mass while keeping its strength. The solutions of the optimization problem have been verified by solving the direct problems with obtained design parameters using the classical and improved shell theories.

• The study has demonstrated acceptability and convenience of using simple mathematical models based on Kirchhoff—Love and Timoshenko shell theories for numerical solving optimization problems.

Acknowledgments

The study was supported by project 18-13-00392 of Russian Science Foundation.