Sage Journals: Discover world-class research

Abstract

Fiber path design plays a pivotal role in determining the forming quality of composite pressure vessels. Current methodologies for asymmetric dome structures predominantly rely on constant slippage coefficient non-geodesic trajectories, which significantly constrain design flexibility. To overcome this limitation, this study proposes a fiber path design method for asymmetric dome composite pressure vessels based on a variable slippage coefficient non-geodesic trajectories. The research framework comprises three principal components: First, development of a novel fiber trajectory model through variable slippage coefficient non-geodesic formulation. Second, systematic establishment of five distinct fiber path design patterns with corresponding implementation schemes for cylindrical sections. Additionally, simulation analyses are conducted on the winding angle, slippage coefficient, and the five fiber path winding types using MATLAB. Finally, under internal pressure strength requirements, Barlow’s formula and netting theory are applied to design the liner and fiber layers, respectively. ABAQUS is employed for modeling and mechanical performance analysis. The results demonstrate the feasibility of the proposed variable slippage coefficient non-geodesic filament winding method for asymmetric dome pressure vessels. Finite element analysis confirms the structural reliability of the composite pressure vessel under internal pressure. This study provides valuable reference for fiber path design in asymmetric dome composite pressure vessels.

Keywords

Fiber path composite pres sure vessels variable slippage coefficient non-geodesic winding types mechanical performance

Introduction

Filament winding is a technique used to fabricate composite structures by winding reinforcing fibers around the surface of a mandrel.^1,2 Due to the excellent strength and stiffness of the products it produces, it is widely used in the manufacturing of high-pressure vessels such as solid rocket motors, liquefied natural gas (LNG) tanks, and oxygen storage equipment.

It is well known that fiber path design is the foundation of composite pressure vessel design and plays a crucial role in ensuring process stability and the high load-carrying capacity of the dome. Currently, research on pressure vessels mainly focuses on: concave composites,^3,4 toroidal pressure vessels,^5,6 bends,⁷ symmetric domes,^8,9 and asymmetric domes.¹⁰ However, symmetric and asymmetric domes are relatively more popular among researchers. In terms of fiber path design for symmetric domes, Taner et al.¹¹ designed a new composite hydraulic cylinder by constructing a constant curvature meridian profile using geodesics. Numerical analysis was conducted using finite element software, followed by response surface methodology to optimally determine the design parameters. The results showed a 53.78% weight reduction compared to conventional steel hydraulic cylinders. Li et al.¹² proposed the concept of constant angle winding on curved surfaces to overcome the limitations of geodesic and semi-geodesic lines in traditional winding methods. By deriving the anti-slip algorithm for constant angle winding curves on curved surfaces, they determined the constant angle winding area and optimized the winding mode of the pressure vessel by combining it with semi-geodesic lines. This approach allowed for constant angle winding in the transition area between the cylindrical and dome sections, ensuring that the fiber completely covers the core axis. Zu et al.¹³ introduced an optimal winding parameter design method for composite pressure vessels based on non-geodesic trajectories. The non-geodesic equations of the pressure vessel were derived using differential theory and winding principles. After determining the polar opening radius, the optimal winding parameters were obtained using an interior-point penalty function to ensure that the fiber met the stability and uniform coverage requirements on the mandrel. Burst pressure was also predicted using finite element analysis and compared with experimental results. Zu et al.¹⁴ further investigated the stress/strain distribution and burst pressure analysis of composite spherical pressure vessels based on non-geodesic winding paths. Zhou et al.¹⁵ proposed a dome shape optimization design method that considers both geodesic and non-geodesic winding methods. Using the shape factor as the objective and the aspect ratio and power index in the hyperelliptic function as design variables, they optimized the relevant parameters using a particle swarm algorithm. The aforementioned scholars have conducted extensive research on fiber path design for symmetric domes, head shape optimization, and finite element modeling. However, there remains a lack of a comprehensive path design and simulation system for asymmetric domes.¹⁶

According to the available literature, although only a few scholars have studied the asymmetric dome pressure vessel, there are usually too many different unequal polar holes in solid rocket motors in practical applications. Asymmetric domes are classified into three cases: unequal polar hole radii with equal ellipsoidal rates on both sides, unequal ellipsoidal rates with equal polar hole radii on both sides, and unequal ellipsoidal rates and polar hole radii on both sides. In terms of the fiber path design for asymmetric domes, Zu et al.¹⁶ have investigated a design method to achieve a continuous transition of the winding angle of the dome with unequal pole holes in the cylindrical section using non-geometric. The calculation of linear parameters was performed by applying the principle of continuous fraction, and the feasibility of this method was verified by simulation and experimentation. Zu et al.¹⁷ provided a design method for determining an isotonic structure with unequal polar openings and evaluated the effect of the non-geodesic on the geometry and performance of an isotonic dome. Guo et al.^18,19 proposed a set of constant slippage coefficient non-geodesic (the ratio of geodesic curvature to normal curvature) design methods for asymmetric fiber-reinforced polymer domes. The method considers five design cases and defines the corresponding design regions. It then establishes model relationships for the line shape, the number of tangent points, and fiber width, which were subsequently verified by a simulation software system developed using MATLAB (MathWorks, USA). However, the research conducted by the aforementioned scholars on the design of dome fiber paths, using constant slippage coefficient non-geodesic lines, has not been extended to the variable slippage coefficient fiber paths used in wrapping asymmetric composite pressure vessels.

Although the above discussion addresses fiber path design methods for composite pressure vessels, the winding and forming process is the key to the practical implementation of these designs. Depending on the properties of the fiber layers, the winding process can be categorized into thermosetting filament winding and thermoplastic filament winding.

In terms of thermosetting filament winding, existing studies have mainly focused on forming characteristics and process optimization. Madhavi et al.²⁰ conducted experimental investigations on the forming behavior of non-geodesic wound composite pressure vessels, analyzing the available coefficient of friction between the fibers and mandrel surface and its variation over time, and established a correlation model between the coefficient of friction and material adhesion. İpekçi et al.²¹ addressed the limitations of conventional filament winding methods, which rely on mold-based forming and make it difficult to achieve orthogonal ply placement, thereby restricting the design freedom of the structure. To overcome these issues, they proposed combining photopolymerizable resin with robotic additive manufacturing techniques to develop a mold-free ultraviolet curing winding process. Overall, thermosetting filament winding is a mature technology that enables strong interlaminar bonding and excellent mechanical performance, with good resin flowability, suitable for complex geometries and thick-walled components, and relatively low material cost. However, its curing is irreversible, thick components require long curing cycles, temperature and pressure must be strictly controlled, and prepregs require refrigerated storage, which limits production efficiency and process flexibility.

Regarding thermoplastic filament winding, this process achieves substantial bonding between the fiber tape and substrate through localized heating, combined with applied pressure or roller compaction, forming a near-infinite friction/no-slippage condition, which allows for complex non-geodesic fiber path winding. Smith et al.²² experimentally evaluated the polymerization behavior of low-exotherm liquid thermoplastic resins in thick-section laminate manufacturing to address critical knowledge gaps in thick composite production. Xu et al.²³ investigated the fiber layer relaxation issues that may arise during high-tension winding of thermoplastic composites. Based on thick-walled cylinder theory and small deformation assumptions, they established a mathematical model for fiber stress distribution and mandrel surface radial compressive stress considering the coupling effects of tension and temperature fields, and proposed a method to balance thermal stresses by applying controlled tension. The results showed that thermal stress can effectively reduce inner layer fiber relaxation, and finite element simulations validated the reliability of the model. Thermoplastic filament winding offers fast forming speeds and good toughness, making it suitable for automated in-situ forming and resulting in fewer voids in thick-section structures. However, thermoplastic resins exhibit high viscosity, interlaminar bonding depends on applied pressure or rolling, thick components require precise heat management, material costs are relatively high, and thick-walled parts may develop internal stresses or warpage. In comparison, thermosetting and thermoplastic filament winding processes each have their advantages and limitations; notably, thermoplastic filament winding can effectively eliminate interlaminar slippage, enabling the successful implementation of non-geodesic fiber paths with high slippage coefficients.

In summary, previous studies have focused on geodesic or constant slippage coefficient non-geodesic path designs for both symmetric and asymmetric domes. However, no research has integrated variable slippage coefficient non-geodesic trajectories for fiber path design in asymmetric dome composite pressure vessels. To address this gap, this study proposes a variable slippage coefficient fiber path design method for asymmetric dome composite pressure vessels, building on the latest advancements in variable slippage coefficient non-geodesic trajectories from our research group.²⁴ The details are as follows: firstly, the dome variable slippage coefficient non-geodesic control equation is given. Secondly, design equations for the liner thickness and fiber layer thickness under internal pressure strength are established. Subsequently, the existing fiber path types at the cylindrical section of the asymmetric composite pressure vessel are discussed, and detailed implementation schemes are provided. Finally, MATLAB-based simulations are conducted to analyze the five fiber path schemes, and ABAQUS (Dassault Systèmes, France) is employed for mechanical performance analysis.

Methods

To facilitate understanding of the subsequent formula derivations, the principal symbols employed in this paper are defined in Table 1.

Table 1.

Principal symbols.

Symbol	Meaning	Unit
$R$	Sphere radius	$m m$
$r_{c}$	Polar hole radius	$m m$
$a$	Ellipsoid major semi-axis	$m m$
$b$	Ellipsoid minor semi-axis	$m m$
$α_{s}$	Initial winding angle of the non-geodesic path on the cylindrical section	°
$α_{e}$	Terminal winding angle of the non-geodesic path on the cylindrical section	°
$θ$	Central rotation angle	°
$L$	Total length of the cylindrical section	$m m$
$α_{1}$	Initial winding angle of the left dome	°
$α_{2}$	Geodesic winding angle of the cylindrical section	°
$α_{3}$	Initial winding angle of the right dome	°
$H_{1}$	Length of the left non-geodesic region of the cylindrical section	$m m$
$H_{2}$	Length of the right non-geodesic region of the cylindrical section	$m m$
$L_{c}$	Length of the geodesic region of the cylindrical section	$m m$

Variable slippage coefficient non-geodesic

By the latest research results from the studio, Wang et al.²⁴ proposed a variable slippage coefficient non-geodesic generated through meridian rotation mapping. The meridian rotation mapping mechanism works as follows: First, the sphere is used as the matrix to calculate the sphere’s meridian contour; second, the meridian is rotated to become tangent to its edges at the polar hole using the rotation matrix; finally, it is converted into a fiber path on the ellipsoid according to the mapping relationship between the parametric equations of the sphere and the ellipsoid. Figure 1 shows the meridian rotation mapping model, with the outermost transparent hemisphere and the inner half-ellipsoid. According to the rotation mapping mechanism, it can be seen that curve 2 is obtained from the meridian 1 through the X-axis rotation matrix transformation.

Figure 1.

Meridian rotation mapping model²⁰.

Curve 2 is generated by applying the X-axis rotation matrix to meridian 1, making it tangent to the polar hole of the sphere. According to the mapping relationship between the sphere and the ellipsoid, curve 3 is then obtained on the ellipsoid. The X-axis rotation matrix can be expressed as²⁵:

R_{x}^{'} = [\begin{array}{l} 1 & 0 & 0 \\ 0 & \cos α^{'} & \sin α^{'} \\ 0 & - \sin α^{'} & \cos α^{'} \end{array}]

(1)

α^{'} = \frac{π}{2} - α

(2)

According to the trigonometric relationship $α$ can be obtained as:

α = a r c \cos (\frac{r_{c}}{R})

(3)where

r_{c}

is the radius of the polar aperture, and the final rotation matrix can be obtained by bringing equation (3) into (1) as:

R_{x} = [\begin{array}{l} 1 & 0 & 0 \\ 0 & \sin α & \cos α \\ 0 & - \cos α & \sin α \end{array}]

(4)

It is well known that the parametric equations for the sphere and ellipsoid, respectively, are expressed as follows:

\vec{G} = [\begin{array}{l} R \cos v_{2} \cos v_{1} & R \cos v_{2} \sin v_{1} & R \sin v_{2} \end{array}]

(5)

\vec{P} = [\begin{array}{l} a \cos u_{2} \cos u_{1} & a \cos u_{2} \sin u_{1} & b \sin u_{2} \end{array}]

(6)where

R

is the radius of the sphere,

v_{1}

is the latitudinal direction angle coordinate, and

v_{2}

is the longitudinal direction angle coordinate.

a

is the long axis of the ellipsoid,

b

is the short axis of the ellipsoid,

u_{1}

denotes the latitudinal direction angle coordinate, and

u_{2}

denotes the longitudinal direction angle coordinate. Since

v_{1}

v_{2}

and

u_{1}

u_{2}

the same range of values, let

v_{1} = u_{1}, v_{2} = u_{2}

. According to equations (5) and (6), the mapping matrix of the two can be known as:

A = [\begin{array}{l} \frac{a}{R} \\ \frac{a}{R} \\ \frac{b}{R} \end{array}]

(7)

According to the meridian rotation mapping mechanism, the final ellipsoidal winding trajectory equation is obtained as:

P^{G} = G • R_{x} • A

(8)

It should be noted that: (1) the ellipsoid is made of a rotatable ellipse; (2) the sphere and the ellipsoid polar aperture radius size is set to the same.

Fiber path equations at the cylinder

The structure of the cylinder section is simple and the fiber path design can be based on both geodesic and non-geodesic methods. Geodesic exhibit a constant winding angle on the cylinder, so the equations for the fiber paths and center angle equations can be obtained according to reference,²⁴ As shown in Figure 2, $L$ and $D$ represent the length and diameter of the cylindrical section, respectively.

\begin{array}{l} w = \frac{π D}{\tan α_{s}} \\ θ = \frac{L}{w} \times 360 ° \end{array}

(9)

Figure 2.

Schematic diagram of fiber trajectories on the cylindrical cross-section.

The non-geodesic needs to be calculated by the Runge-Kutta method. According to Ref. 19:

L = L_{0} + \frac{R}{λ} (\frac{1}{\sin α_{s}} - \frac{1}{\sin α_{e}})

(10)

\frac{d θ}{d z} = \frac{\tan α}{R}

(11)where

α_{s}

and

α_{e}

are the initial and termination winding angles of the cylinder non-geodesic, respectively, and

θ

is the mandrel center angle.

Barlow’s equation

Barlow’s equation is employed to determine the liner thickness of pressure vessels, and the detailed formulation is presented as follows²⁶:

t_{θ 1} = ξ_{1} \frac{PR}{σ_{1}}

(12)where

ξ_{1}

is the proportional coefficient, set to 0.2,

σ_{1}

is the liner’s ultimate strength, and

P

is the burst pressure.

Netting theory

In academic and engineering applications, netting theory is widely used as a simple and suitable method for the design of filament-wound pressure vessels.²⁶ Assuming that the helical and hoop fiber layers fail simultaneously under burst pressure, the thickness formulas for the helical and hoop fiber layers in the cylindrical section under burst pressure are as follows:

t_{a} = K \frac{R P}{2 δ σ \cos^{2} α}

(13)

t_{θ} = \frac{R P}{2 σ} (2 - \tan^{2} α)

(14)where

t_{a}

and

t_{θ}

represent the total thickness of the helical and hoop layers, respectively.

K

indicates the strength reinforcement coefficient, and the value range is 1.05–1.4.

R

denotes the outside radius of the cylinder section,

δ

is generally set to 0.65–0.85,

P

refers to the burst pressure,

σ

is fiber strength, and

α

is the winding angle of the cylinder section.

Dome thickness prediction

According to the existing thickness prediction algorithm, the cubic spline function method is considered to be the most accurate.^27,28 It is divided into two parts, as follows:

(1) In the two bandwidth ranges, fiber thickness is calculated as follows:

t (r) = m_{1} + m_{2} \times r + m_{3} \times r^{2} + m_{4} \times r^{3}, (r_{0} \leq r \leq r_{2 b})

(15)

(2) Beyond the two bandwidths, the fiber thickness is calculated as follows:

t (r) = \frac{m_{R} n_{R}}{π} [\arccos (\frac{r_{0}}{r}) - \arccos (\frac{r_{b}}{r})] t_{p}, (r_{2 b} \leq r \leq R)

(16)where the coefficients m₁, m₂, m₃, and m₄ should be obtained by four boundary conditions,

t (r)

indicates the fiber thickness at any point on the dome,

r

denotes the radius of any parallel circle on the dome,

r_{0}

refers to the polar hole radius,

b

represents the width of the fiber band,

t_{p}

is the thickness of a single-fiber band,

r_{b}

indicates the radius of a parallel circle at one bandwidth of the dome,

r_{2 b}

is the radius of a parallel circle at two bandwidths of the dome,

α_{0}

denotes the angle of winding at the dome’s equator,

m_{R}

represents the number of fiber bands of the cylindrical section,

n_{R}

indicates the number of helically wound single layers of the cylindrical section,

m_{0}

is the number of fiber bands of the polar hole,

n_{0}

refers to the number of winding single layers at the polar hole, and

t_{R}

is the fiber thickness of the helical layer in the cylindrical section.

Linear design of asymmetric dome composite material pressure vessels

Type of line on the cylinder

When designing fiber paths for pressure vessels, the winding angle of the cylinder section is typically set to a constant value. However, this is not always feasible when the two sides of the dome have unequal ellipsoidal rates or polar hole radii. Therefore, it becomes necessary to design the cylinder section winding angle as variable to bridge the initial winding angle of the asymmetric domes on both sides. In this paper, the dome fiber path is determined to follow a variable slippage coefficient non-geodesic design mode, so only the cylinder section needs to be designed with a linear winding scheme. Figure 3 illustrates the design of the linear winding scheme for the asymmetric dome composite pressure vessel. According to the literature,^18,19 there are five types of linear designs for the cylinder, which are classified as follows: (1) full geodesic; (2) left geodesic, right non-geodesic; (3) left non-geodesic, right geodesic; (4) left and right non-geodesic, and geodesic in the middle; and (5) full non-geodesic. The corresponding design scheme is determined according to the actual requirements and linear parameters.

Figure 3.

Linear design of asymmetric dome composite pressure vessel.

The following conditions must be met during the filament winding process: (1) the endpoint of the fiber after a complete cycle should be offset from the starting point by one fiber bandwidth. (2) the slippage coefficient value corresponding to each point on the fiber path should not exceed the specified maximum friction coefficient value. (3) the winding angle should be continuous at the connection between the head and the cylinder, without abrupt changes.

The important parameters at the cylinder are $α_{1}, α_{2}, α_{3}, H_{1}, H_{2}, L$ . $α_{1}, α_{2}, α_{3}$ are the initial winding angle of the left dome, the winding angle of the middle geodesic of the cylinder body section, and the initial winding angle of the right dome, respectively. $H_{1}, H_{2}, L$ are the length of the left non-geodesic region of the cylinder section, the length of the right non-geodesic region of the cylinder section, and the total length of the cylinder section, respectively. In conjunction with Figure 4, these five linear design options and the related parameter solution process can be described.

(1) When $α_{1} = α_{2} = α_{3}$ , the cylinder section becomes a full geodesic section, i.e., the length of the non-geodesic region on the left and right sides is zero $(H_{1} = H_{2} = 0)$ . At this point, it is necessary to solve the geometric parameters (ellipsoidal rate and polar hole radius) of the domes on both sides, ensuring that the domes on both sides are asymmetric. According to Figure 3, the bandwidth is used as a constraint to solve the polar hole radius of the left dome, and then the ellipsoidal rate and the radius of the polar hole of the right head are calculated with the constraint that the initial winding angles of the two domes are equal and ensure that the two domes are asymmetric. Finally, when the values of fiber path slippage coefficients of the left and right domes are less than the maximum friction coefficient, the geometrical parameters of the left and right domes and the number of tangent points are output.

(2) When $α_{1} = α_{2} \neq α_{3}$ , at this time, the left side of the cylinder section is a geodesic section, and the right side is a non-geodesic section, i.e., $H_{1} = 0, H_{2} \neq 0$ . Preset the geometric parameters of the left dome (ellipsoidal rate and polar hole radius) and the $u_{r}$ , and the length of the right side of the cylinder section, non-geodesic length $H_{2}$ can be calculated according to equation (17), then the length of the left side of the cylinder section, geodesic length, can be obtained by the calculation of equation (18). Taking the bandwidth and stable winding as constraints, the polar hole radius of the right dome and the number of tangent points are output.

H_{2} = | \frac{R}{u_{r}} (\frac{1}{\sin α_{2}} - \frac{1}{\sin α_{3}}) |

(17)

L_{c} = L - H_{2}

(18)where

u_{r}

is the set value of the right non-geodesic slippage coefficient.

(3) When $α_{1} \neq α_{2} = α_{3}$ , the left side of the cylinder section is a non-geodesic section, and the right side is a geodesic section, i.e., $H_{1} \neq 0, H_{2} = 0$ . Preset the left side dome geometrical parameters and $u_{l}$ , according to equation (19) can calculate the non-geodesic length $H_{1}$ of the left side of the cylinder section, then the geodesic length of the right side of the cylinder section can be calculated by equation (20). The bandwidth and stable winding are used as constraints to output the polar hole radius of the right side dome and the number of tangent points.

H_{1} = | \frac{R}{u_{l}} (\frac{1}{\sin α_{1}} - \frac{1}{\sin α_{2}}) |

(19)

L_{c} = L - H_{1}

(20)where

u_{l}

is the value of the left non-geodesic slippage coefficient.

(4) When $α_{1} \neq α_{2} \neq α_{3}$ , at this time, the left and right sides of the cylinder section are non-geodesic sections, and the middle section is geodesic section. The geometrical parameters of the left dome and the values of the constant slippage coefficients ( $u_{l}$ and $u_{r}$ ) on the left and right sides of the cylinder section are preset, and the non-geodesic lengths $H_{1}$ and $H_{2}$ on the left and right sides of the cylinder section can be calculated according to equations (19) and (17), and the geodesic lengths of the middle section of the cylinder section are calculated by equation (21). With the bandwidth and stable winding as constraints, the polar hole radius of the right dome and the number of tangent points are determined. Then the length of the middle geodesic section $L_{c}$ is:

L_{c} = L - H_{1} - H_{2}

(21)

(5) When $α_{2}$ does not exist, there is no geodesic section in the cylinder section at this time, and all are non-geodesic sections. At this time, in order to avoid the complexity of the calculation, the value of the fiber path slippage coefficient of the cylinder section is directly set to a constant value u. According to the preset geometrical parameters of the left dome in the bandwidth and the stable winding as the constraints, the polar hole radius of the right dome and the number of the tangent points are output.

Figure 4.

Flowchart of winding parameter solution.

Constraint relationship

In this paper, the left and right dome fiber paths are designed with variable slippage coefficient non-geodesic and when the dome geometries are given, the filament winding pattern is already determined. At this point, it is necessary to determine whether the design width and the actual width in the winding pattern are equal, where the design width is based on the literature,²⁴ and the expression can be obtained as:

B = (4 γ + 2 θ - (\frac{1}{n} + N) \times 2 π) n R \cos α_{s}

(22)where

γ

is the central rotation angle of the dome fiber path,

θ

is the rotation angle of the cylinder fiber path,

n

is the number of tangent points, and

N

is a non-negative integer. When evenly distributed, the number of filament winding back and forth is:

{\begin{cases} M = c e i l (\frac{2 π R \cos α_{s} + B}{n B}) \\ S = n \cdot M \end{cases}

(23)where

M

is the number of complete cycles when the fibers are uniformly covered, and

S

is the total number of cycles when the fibers are uniformly covered.

Simulation and analysis

Slippage coefficient analysis

A further analysis of the slippage coefficient distribution of the shell under different geometric configurations in Section 2.1 helps determine whether the slippage coefficient corresponding to the optimal winding parameters generated in Figure 4 exceeds the allowable range. According to the slippage coefficient calculation formula in the literature,²⁴ the slippage coefficient distribution as shown in Figure 4 is obtained. Let $R = a = 1$ , $m = 0.85$ , the slippage coefficient distribution of Figure 5(a) is obtained for different polar hole radii under the parameter radian $u_{1}$ . From the figure, it can be seen that the slippage coefficients on the fiber path are continuously changing and centrosymmetric on the left and right sides, and the slippage coefficients are 0 when a is equal to 0 or $π$ . With the increase of polar hole radius, the tip arcs on both sides are decreasing, and the intermediate curves show an increasing trend. Figure 5(b) shows the maximum slippage coefficient values under different polar hole radii, where the maximum value is 1.2138 and the minimum value is 0.0922. As the polar hole radius increases, the maximum slippage coefficient value decreases nonlinearly. Therefore, the optimum polar hole radius can be determined based on the maximum slippage coefficient value. Let $R = a = 1$ , $r_{c} = 0.4$ , the distribution of slippage coefficients for different ellipsoidal rates under parameter radian $u_{1}$ in Figure 5(c) is obtained. From the figure, it can be seen that the slippage coefficients on the fiber paths are continuously changing, and the same slippage coefficient is 0 when $u_{1}$ is equal to 0 or $π$ . When $m$ is in the range of 0.7∼1, the variable slippage coefficient curves show a decreasing trend as a whole. m shows an increasing trend as a whole when $m$ is in the range of 1∼1.3. Figure 5(d) shows the maximum slippage coefficient values under different ellipsoidal rates, and the overall trend of the variable slippage coefficient curve in Figure 5(c) is consistent with the first decline and then increase. It should be noted that the slippage coefficient of the fiber path on the dome is constant 0 at $m = 1$ , and the fiber path is constructed by the geodesic. As shown in Figure 5(b) and (d), once the dome geometric parameters are determined, the continuous variation curve of the slip coefficient along the fiber path is also fixed.

Figure 5.

Slippage coefficient distribution of dome.

Winding angle analysis

It is well known that the distribution of the winding angle in the dome portion of a pressure vessel is very important for its axial and radial load carrying capacity. The size of the winding angle indicates the degree of deviation of the fiber path from the meridian. The larger the winding angle, the greater the deviation, and consequently, the greater the radial load capacity of the pressure vessel. Similarly, the smaller the winding angle and the smaller the degree of deviation, the stronger the axial load capacity. Let $R = a = 1$ , $m = 0.85$ and $R = a = 1$ , $r_{c} = 0.4$ . Figure 6(a) and (b) show the distribution of the dome winding angle under different radii and ellipsoidal rates. When the polar hole radius or ellipsoidal rate is fixed, the winding angle increases monotonically with the increase of the dome height, and the increase of the winding angle near the polar hole slows down. When the value of dome height is fixed, the winding angle shows a decreasing (or increasing) trend with the increase of the ellipsoidal rate (or polar hole radius). From the figures, it can be observed that the effect of the polar hole radius on the initial winding angle is larger. It can also be concluded that when the geometrical parameters of the dome are determined, the distribution curve of the winding angle of the fiber path on the dome is also fixed.

Figure 6.

Dome winding angle for variable geometrical parameters.

Simulation of winding forming process

In order to visualize the winding process, the first design path is chosen as a case study, and the specific winding parameters are shown in Table 2. Figure 7 shows the effect of pressure vessel filament winding under different numbers of complete cycles. When M = 1 is a complete cycle, the number of tangent points is 5, which are tangent to the polar hole radius and are in equal position. As shown in Table 2, the bandwidth error of the beginning and end fiber spacing after one cycle is 0.3%, which is in agreement with the actual design requirements. With the increase of the number of complete cycles, the fiber is gradually and uniformly wrapped all over the surface of the mandrel. When M = 15, the fiber has been completely and uniformly covered the surface of the mandrel. From the figure, it can be seen that there is a serious accumulation of fiber near the dome polar hole region, which is in line with the actual objective facts.^29,30 This winding effect proves the correctness of the procedure.

Table 2.

Main winding parameters of five fiber paths.

	1	2	3	4	5
Left dome
$m_{1}$	1.1	0.8	0.85	0.9	0.9
$r_{1}$	0.3877 R	0.65 R	0.7 R	0.7 R	0.6 R
$α_{1}$	20.92°	46.91°	49.07°	47.44°	39.71°
$λ_{1 \max}$	0.3530	0.1646	0.1151	0.0809	0.1071
Cylinder section
$u_{l}$	0	0	0.1570	0.26	0.3906
$α_{2}$	20.92°	46.91°	40.36°	45°	/
$u_{r}$	0	0.28	0	0.2	0.3906
Right dome
$m_{2}$	0.9314	0.8	0.85	0.9	0.9
$r_{2}$	0.3355 R	0.5187 R	0.5856 R	0.5960 R	0.4148 R
$α_{3}$	20.92°	37.18°	40.36	39.51°	26.86°
$λ_{2 \max}$	0.2368	0.2416	0.1570	0.1084	0.2153
$n$	5	2	2	2	3
Error	0.3%	0.11%	0.026%	0.098%	0.3%

Figure 7.

Filament winding effect of pressure vessel under different complete cycles.

Five winding pattern simulation

Firstly, the basic parameter $R = 60 m m, L = 100 m m, b = 5 m m, u_{s} = 0.4$ is preset. The winding parameters shown in Table 2 can be obtained by the winding parameter calculation process, where the errors in the table refer to the differences between the theoretical and actual bandwidths. The axial distribution of the winding angle under the five fiber paths shown in Figure 8 is drawn by MATLAB. When $α_{2}$ is the first design case, the cylindrical section represents a geodesic segment. Thus, by setting the ellipsoidal rate of the left dome to $m_{1} = 1.1$ and the polar opening radius to $r_{1} = 0.5 R$ , the pressure vessel can be regarded as a symmetric shell at this stage. Based on this, the fminsearch function is used to solve for the left polar opening radius $r_{1} = 0.3877 R$ so that the winding principle is satisfied (i.e., the difference between the theoretical and actual bandwidths is minimized). Then, using the determined geometric parameters of the left dome, the initial winding angle is obtained. Taking this initial winding angle as the target, the fminsearch function is further employed to determine the ellipsoidal rate and polar opening radius $m_{2} = 0.9314, r_{2} = 0.3355 R$ of the right dome, minimizing the difference from the initial winding angle. As shown in the figure, the winding angle on the left dome gradually decreases along the axial direction, remains constant along the cylindrical section, and gradually increases along the axial direction on the right dome. The winding parameters in Table 2 clearly show that in the first fiber path design $m_{1} \neq m_{2}, r_{1} \neq r_{2}$ , the left and right sides of the dome are asymmetric. When $α_{2}$ is the second design case, the left side of the cylindrical section represents a geodesic segment, while the right side represents a non-geodesic segment. The ellipsoidal rate and polar hole radius of the left dome are set as constant values ( $m_{1} = 0.8, r_{1} = 0.65 R$ ), the non-geodesic slippage coefficient on the right side of the cylindrical section is $u_{r} = 0.28$ , the ellipsoidal rate of the right dome is $m_{2} = 0.8$ , and the polar hole radius of the right dome is preset as $r_{2} = 0.4 R$ . Based on the geometric parameters of the left and right domes, the initial winding angles and central rotation angles for both sides can be calculated. The geodesic winding angle on the left side of the cylindrical section is consistent with the initial winding angle of the left dome. Using equations (17) and (18), the lengths of the left geodesic segment and the right non-geodesic segment of the cylindrical section can be determined. Equations (9) and (11) are then used to calculate the overall central rotation angle of the cylindrical section. These calculation results serve as internal parameter inputs for the bandwidth objective function. Therefore, by setting the bandwidth as the target and optimizing under stability constraints, the polar hole radius of the right dome is obtained as $r_{2} = 0.5187 R$ . Combined with Table 2, it can be concluded that the left and right sides of the dome ellipsoidal rate is equal, and the polar hole radius is not the same size. At this time, the left side of the cylinder section is a geodesic section, and the right side is a non-geodesic section. The winding angle of the cylinder section is first kept constant and then gradually reduced. When $α_{2}$ is the third design case, the left side of the cylinder section is the non-geodesic section, and the right side is the geodesic section, the optimization details of the winding parameters are similar to those in the second design case. At this time, the winding angle of the cylinder section gradually decreases and then remains unchanged. When $α_{2}$ is the fourth design case, the left side of the cylinder section is the non-geodesic section, the middle is the geodesic section, and the right side is the non-geodesic section. At this time, the winding angle of the barrel section first decreases, then remains unchanged and finally gradually decreases. The left dome’s ellipsoidal rate and polar hole radius are set as constant values ( $m_{1} = 0.9, r_{1} = 0.7 R$ ), while the cylindrical section has a geodesic winding angle $α_{2} = 45^{\circ}$ . The non-geodesic slippage coefficients on the left and right sides of the cylindrical section are $u_{l} = 0.26$ and $u_{r} = 0.2$ , respectively. The right dome shares the same ellipsoidal rate as the left dome, and its polar hole radius $r_{2} = 0.5 R$ is preset. Based on the geometric parameters of both domes, the initial winding angles and central rotation angles on each side can be determined. Combined with the geodesic winding angle of the cylindrical section and differential geometry theory, the non-geodesic lengths on both sides of the cylindrical section and the overall central rotation angle can be further obtained. By taking the bandwidth as the objective function and applying stability constraints, the optimized right dome polar hole radius $r_{2} = 0.5960 R$ is determined. As shown in the figure, the winding angle of the cylindrical section first decreases, then remains constant, and finally decreases again. When $α_{2}$ does not exist, the cylinder section is a non-geodesic section. By setting the ellipsoidal rate and pole hole radius of the left dome as constant values ( $m_{1} = 0.9, r_{1} = 0.6 R$ ), and assuming that the ellipsoidal rate of the right dome is identical to that of the left dome, with a preset pole hole radius of the right dome $r_{2} = 0.4 R$ , the initial winding angles on both sides can be determined based on the geometric parameters of the two domes. The slippage coefficient of the non-geodesic section in the cylindrical part can then be calculated using equation (10). Following the same winding parameter optimization procedure, the pole hole radius of the right dome $r_{2} = 0.4148 R$ can be further optimized. The winding angle on the left dome and the cylinder section gradually decreases, and the winding angle on the right dome gradually increases to 90°.

Figure 8.

Axial distribution of winding angle under five fiber paths.

From the above five path design cases, it can be concluded that: First, both the left and right side domes are asymmetric. Second, the slippage coefficient on the left and right sides of the dome varies continuously along the fiber path, that is, variable slippage coefficient non-geodesic. Third, the slippage coefficient of the composite pressure vessel at any point along the fiber path is always less than 0.4, meaning that slipping will not occur during actual production. Of course, determining which specific linear design is more suitable for engineering applications requires considering process parameters, dome segment variable thickness, curing processes, and many other important factors, which will be studied in the next step.

Finite element analysis

Since the liner thickness is calculated using Barlow’s formula, while the fiber layer thickness of the cylindrical section and the dome is determined via netting theory and cubic spline functions, respectively, the design theory can be applied to all composite pressure vessels with linear configurations shown in Figure 8. This study selects the case in Figure 8(a) for internal pressure strength design and finite element analysis. The simulated working internal pressure is 30 MPa, with T700 carbon fiber as the fiber layer material and aluminum alloy 6061 as the liner material, as shown in Tables 3 and 4.

Table 3.

Mechanical parameters of aluminum alloy 6061³¹.

Material	Elastic modulus (GPa)	Poisson’s ratio	Yield strength (MPa)	Tensile strength (MPa)
6061	69	0.324	298	330

Table 4.

Material properties of T700/epoxy composite²⁷.

Property	T700/epoxy composite
Elastic modulus $E_{X}$	134000 MPa
Elastic modulus $E_{Y}$	7420 MPa
Elastic modulus $E_{Z}$	7420 MPa
Shear modulus $G_{12}$	3710 MPa
Shear modulus $G_{23}$	3710 MPa
Shear modulus $G_{13}$	4790 MPa
Poisson’s ratio $u_{12}$	0.28
Poisson’s ratio $u_{23}$	0.28
Poisson’s ratio $u_{13}$	0.3
Longitudinal tensile strength $X_{t}$	2300 MPa
Longitudinal compressive strength $X_{c}$	1250 MPa
Transverse tensile strength $Y_{t}$	74 MPa
Transverse compressive strength $Y_{c}$	180 MPa
Shear strength $S$	50 MPa
Fiber strength $σ_{f}$	2800 MPa
Fiber band thickness	0.227 mm

According to DOT-CFFC requirements, the burst pressure is set to 2.8 times the working pressure (84 MPa). For the cylindrical section in Figure 8(a), with a winding angle of 20.92°, equation (12) yields a liner thickness of approximately 3 mm. Using equations (13) and (14), the thicknesses of the hoop and helical layers are calculated as 1.67 mm and 1.38 mm, respectively. Given a single fiber layer thickness of 0.227 mm, the required numbers of hoop and helical winding layers are 7.356 and 6.079, respectively. Considering that filament winding typically employs cross-ply double layers, both the hoop and helical layers are rounded up to 8 layers. The thickness of the left and right dome sections is calculated using cubic spline functions, as shown in Figure 9.

Figure 9.

Thickness distribution of dome fiber layers.

Due to the structural complexity and material anisotropy of composite pressure vessels, mechanical performance analysis is challenging. This study employs the finite element software ABAQUS to evaluate the mechanical behavior under both working and burst pressures. A full-scale finite element model of the composite pressure vessel was established using three-dimensional modeling in ABAQUS. The liner was modeled as a 3D solid structure, meshed with hexahedral elements, totaling 10,900 elements of type C3D8R. The fiber layers were modeled as 3D shell structures, meshed with quadrilateral elements, totaling 7254 elements of type S4R. The filament winding angles calculated from Figure 8(a) were applied to the corresponding axial positions in the fiber layer model. The boundary conditions of the finite element model were determined according to the actual working conditions: one end of the vessel was fully fixed, while the other end was axially constrained. A tie constraint was applied between the liner and the fiber layers. The material properties of the liner and fiber layers were assigned according to Tables 3 and 4, respectively. Figure 10 illustrates the stress distribution of the composite pressure vessel under variable internal pressures. At the working pressure of 30 MPa, the maximum von Mises stress in the liner is 263 MPa, and the maximum principal stress in the fiber layer is 359.6 MPa. Under the burst pressure of 84 MPa, the maximum von Mises stress in the liner increases to 318 MPa, while the maximum principal stress in the fiber layer rises to 1421 MPa. As shown in the figure, at the working pressure, the maximum von Mises stress in the liner does not reach its yield strength, whereas at the burst pressure it exceeds the yield strength, satisfying the relevant DOT-CFFC requirements. It can also be observed that as the internal pressure increases, the maximum principal stress in the fiber layer at the burst pressure becomes 2.95 times higher than that at the working pressure. Therefore, once the liner’s von Mises stress reaches a certain level, the load is gradually transferred from the liner to the fiber layer.

Figure 10.

Stress distribution in composite pressure vessel under variable internal pressure.

During the bursting process, the fracture strain of the hoop winding layer is 85% of that of the pure fiber, and that of the helical-wound layer is 75% of the pure fiber fracture strain. Given that the fracture strain of the pure fiber is 0.018, the fracture strain of the hoop winding layer and the helical-wound layer are 0.0153 and 0.0135, respectively. Figure 11 illustrates the strain distribution of the hoop layer under variable internal pressure. As shown in the figure, the maximum strain of the fiber layer under the burst pressure of 84 MPa is 0.01026. Therefore, according to the maximum strain criterion, the composite pressure vessel does not experience burst failure under the burst pressure. This validates the effectiveness and reliability of the design for the liner thickness and fiber layer thickness under internal pressure strength requirements.

Figure 11.

Strain distribution in the hoop layer under variable internal pressure.

Conclusion

In this paper, a new fiber path design method for asymmetric dome composite pressure vessels is proposed, which uses variable slippage coefficient non-geodesics, constant slippage coefficient non-geodesics, and geodesics.

First, the left and right dome fiber paths are determined as variable slippage coefficient non-geodesics, and the corresponding control equations are presented. Secondly, five fiber path design schemes for the cylindrical section are proposed, along with a calculation workflow for winding parameters under stability and bandwidth constraints. Additionally, MATLAB-based simulations are conducted to analyze the winding angle and slippage coefficient distribution of the dome under variable geometric parameters. Simulations of axial winding angles and winding patterns are also performed for cases derived from the five fiber path design modes. Finally, using the optimized winding parameters and the liner/fiber layer thicknesses designed for internal pressure strength, a 3D model of the composite pressure vessel is constructed in ABAQUS. Mechanical performance analysis and strength validation are carried out for the liner and fiber layer under variable internal pressures. The results indicate that the fiber path simulation process exhibits no anomalies such as sudden changes in winding angles, irregular fiber distribution, or localized severe overlaps, thereby validating the rationality of the fiber path design. According to the maximum strain criterion, the composite pressure vessel does not experience burst failure under the burst pressure, confirming the reliability of its internal pressure strength design. This study provides valuable insights for the fiber path design and internal pressure strength optimization of asymmetric dome composite pressure vessels in practical engineering applications.

Footnotes

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This research is supported by the Fundamental Research Funds for the Central Universities (Grant No. 2232024G-14), Opening Foundation of Shanghai Collaborative Innovation Center for High Performance Fiber Composites and Natural Science Foundation of Shanghai (Grant No. 22ZR1401200, Grant No. 23ZR1400500).

ORCID iD

Jun Hu

References

Zhang

, et al. Design and production of filament-wound composite square tubes. Compos Struct 2018; 191: 202–208.

Peters

. Composite filament winding. ASM International, 2011.

Qiaoguo

Lei

Qian

, et al. Filament winding analysis and experimental verification of a combined revolution body with a concave surface. Compos Struct 2022; 280: 114951.

Min

Geng

. Stability of filament-wound hyperbolic flexible pipes under internal pressure based on non-geodesic winding. Sci Rep 2021; 11(1): 6043.

. Design and analysis of winding parameters of composite annular high pressure vessel. J Astro 2006; 06: 1350–1355.

. Equicohesive winding linear design and stability analysis on toroidal pressure vessel. J Solid Rocket Technol 2006; 02: 146–149.

Zhang

, et al. Design of filament-wound composite elbows based on non-geodesic trajectories. Compos Struct 2018; 189: 635–640.

Rong

Bin

. Fiber trajectories design of ellipsoid component based on topological mapping methodology. Appl Compos Mater 2017; 24(3): 661–674.

Wang

Song

, et al. Determining the optimal dome shape using a novel variable slippage coefficient non‐geodesic. Polym Compos 2024; 45(13): 11730–11742.

10.

Wang

. Stress analysis and optimization for filament wound ellipsoidal pressure vessels with unequal polar openings. Aeros Mat Tech 2014; 44(04): 25–30.

11.

Taner

Sinan

. Design of the composite hydraulic cylinder with geodesic dome trajectory: a numerical study. Polym Compos 2022; 43(9): 5894–5907.

12.

Haisheng

Mingkun

. Constant winding angle curve on revolution surface and its application. Comput Aided Des 2022; 144: 103160.

13.

Wang

, et al. Design and analysis of filament-wound composite pressure vessels based on non-geodesic winding. Compos Struct 2019; 207: 41–52.

14.

Zhang

, et al. Design of filament-wound spherical pressure vessels based on non-geodesic trajectories. Compos Struct 2019; 218: 71–78.

15.

Zhou

Chen

Zheng

, et al. Dome shape optimization of filament-wound composite pressure vessels based on hyperelliptic functions considering both geodesic and non-geodesic winding patterns. J Compos Mater 2017; 51(14): 1961–1969.

16.

Xiao

Zhang

, et al. Research on filament winding line path and technology of unequal pole hole. Comp Sci Eng. 2021; 07: 48–54.

17.

Koussios

Beukers

. Design of filament-wound isotensoid pressure vessels with unequal polar openings. Compos Struct 2009; 92(9): 2307–2313.

18.

Guo

Wang

Wen

, et al. Winding pattern design of fiber reinforced resin polymer composites winding vessels with unequal polar openings based on non-geodesics. Acta Mater Compos Sin. 2019; 36(05): 1189–1199.

19.

Guo

Wen

Xiao

, et al. Design of winding pattern of filament-wound composite pressure vessel with unequal openings based on non-geodesics. J Eng Fib Fabr 2020; 15: 1558925020933976.

20.

Madhavi

Rama Swamy

Rao

, et al. Experimental evaluation of available surface friction factor on a hoop substrate for a filament-wound composite pressure vessel. J Thermoplast Compos Mater 2011; 24(2): 279–297.

21.

Ahmet

Bülent

. An innovative composite elbow manufacturing method with 6-axis robotic additive manufacturing for fabrication of complex composite structures. J Thermoplast Compos Mater 2024; 37(1): 402–425.

22.

Smith

ILM

Lykkeberg

Garden

, et al. Experimental analysis of alternative thermoplastic resins for manufacturing thick-section composite laminates in tidal turbine blades. J Thermoplast Compos Mater. 2025; 1–19.

23.

Wang

Song

, et al. Analysis of fiber stresses during high-tension winding molding of thermoplastic composites. J Thermoplast Compos Mater 2025; 38(10): 3802–3824.

24.

Wang

Song

, et al. Design method of variable slippage coefficient non-geodesic filament-wound composite pressure vessels. Polym Compos 2024; 45(05): 3950–3964.

25.

Smith

. Engineering computation with MATLAB. Addison Wesley, 2007.

26.

Sharma

Bera

Semwal

, et al. Theoretical analysis of design of filament wound type 3 composite cylinder for the storage of compressed hydrogen gas. Int J Hydrogen Energy 2020; 45(46): 25386–25397.

27.

Wang

Song

, et al. Layup design and mechanical properties of an optimal dome profile. Polym Compos 2024; 45(16): 14993–15004.

28.

Jin

, et al. Failure behavior and influencing factors of volume characteristics of filament wound pressure vessel based on refined model. Comp Sci Eng 2021; 8: 40–47.

29.

Wang

Weicheng

Wenbo

, et al. A new method for predicting dome thickness of composite pressure vessels. J Reinforc Plast Compos 2010; 29(22): 3345–3352.

30.

Zhang

Jia

, et al. Design of a 70 MPa type IV hydrogen storage vessel using accurate modeling techniques for dome thickness prediction. Compos Struct 2020; 236: 111915.

31.

Zhang

Jiang

. Study on structural design and mechanical properties of 70 MPa carbon fiber wound hydrogen storage cylinder. Press Vessel Technol 2023; 40(03): 28–37.

Fiber path design and mechanical performance analysis of asymmetric dome composite pressure vessels

Abstract

Keywords

Introduction

Methods

Variable slippage coefficient non-geodesic

Fiber path equations at the cylinder

Barlow’s equation

Netting theory

Dome thickness prediction

Linear design of asymmetric dome composite material pressure vessels

Type of line on the cylinder

Constraint relationship

Simulation and analysis

Slippage coefficient analysis

Winding angle analysis

Simulation of winding forming process

Five winding pattern simulation

Finite element analysis

Conclusion

Footnotes

Declaration of conflicting interests

Funding

ORCID iD

References