Design of winding pattern of filament-wound composite pressure vessel with unequal openings based on non-geodesics

In this article, we proposed a new approach to design the winding patterns of filament-wound composite pressure vessel with unequal polar openings with non-geodesics. To ensure the continuity of winding angles between trajectories along the cylinder and the dome, the non-geodesics for cylindrical part were used. The developed winding patterns of the vessels were simulated using the MATLAB software to verify the feasibility of the acquired trajectories. To demonstrate the performance in designing the winding path for big polar ratios, we analyzed vessels with polar pole ratios of 1:2–1:4, respectively. The developed winding patterns have successfully achieved uniform fiber distributions along the mandrel without severe overlap, except for the polar pole regions. To avoid the severe overlap between filament bands, we further studied the relationship among the winding pattern, bandwidth, and the number of tangent points, and derived a suitable bandwidth based on the winding pattern. These simulated results proved the effectiveness of the developed method in design of winding pattern with unequal polar openings.


Introduction
Filament winding is a technology for producing composite parts by winding reinforcing filaments on the mandrel surface of a desired shape. [1][2][3] Products manufactured by filament winding technology have an excellent strength-to-weight and stiffness-to-weight ratio and can be used as high-pressure vessels in solid rocket motors, fuel tanks of fighters, liquefied natural gas tanks, torpedo launchers, and so on. [4][5][6][7][8] In fabrication, the continuous fibers can be oriented to match the direction of major stress in composite vessels to improve the carrying capacity. However, these applications in the aerospace field are very sensitive to weight. To lower the launch costs, it is urgent to reduce the weight of each subsystem in the rocket. To produce the lightweight solid rocket, digital design of the winding pattern is the first step, and plays a key role in weight reduction. In a filament winding system, fiber tows under tension can be wound on the mandrel surface in designed patterns to form a uniform composite structure. 9 As a result, the winding pattern determines winding angle distribution, which has a considerable effect on the mechanical properties of composite vessels. [10][11][12] Therefore, an optimized winding pattern can improve the carrying capacity, and simultaneously reduce the weight. 11,[13][14][15] The method of winding pattern design can be split into two ways: by using an analytical or a discrete method. For the analytical method, the winding patterns are calculated by differential equations. 16,17 In the early stage, due to the conciseness, geodesics were widely used in the composite industry. 18 However, the geodesic winding parameters are unique when the mandrel geometry is determined, which restricts the design scope of winding patterns. To broaden the designability, the non-geodesics winding is more and more used in the filament winding on complex parts 19,20,21 By employing the friction between the filament bundles and its supporting surface, the non-geodesics significantly enlarge the available winding path. Koussios 22,23 derived the basic equation for non-geodesic filament winding on generic shells of revolution, which laid the foundation for winding patterns design. For filament winding of complex shape mandrels, a unified approach 24 was developed to simulate the filament winding process. However, this method is more suitable for convex and concave geometries. The spline method 25 for winding pattern design was proposed for non-axisymmetric winding in the composite industry. Zu et al. 26,27 proposed the winding pattern design of pressure vessels, and further improved the dome section performance. The discrete method subdivided the whole revolved surface into a sequence of sections to calculate the fiber trajectories. An algorithm was proposed for nonaxisymmetric structures comprising rectangular facets that compute the filament paths. 28 However, this method is limited to a polygonal cross-section mandrel. To overcome the defect, Fu et al. 29 used a triangular mesh termed as the stereolithography (STL) model to calculate winding paths for both axisymmetric and non-axisymmetric mandrels. Furthermore, Fu et al. 9,30 developed some different methods for generation of winding path for abnormal mandrel based on STL model, including the method of winding path generation based on the inverse process of stability analysis and that based on the principal stress. In the principal stress field method, variable winding directions as close as possible to the major principal stress directions of the product are calculated to generate filament-winding paths that can endure the maximum load and satisfy the requirements of the winding manufacture. In the method of winding path generation based on the inverse process of stability analysis, according to the inverse process of stability analysis, the next path point satisfied the stability conditions. This method can be applied to both axisymmetric and non-axisymmetric mandrel models represented by triangular patches. However, these methods require a lot of calculation effort and are not suitable for fast winding pattern designs of filament-wound composite pressure vessel with two unequal polar poles in the early stages. 25 Overall, most of these methods focused on the subpart sections or the full part of the mandrel for winding pattern design and overlooked the transition zones of fiber trajectories between the dome and the cylinder. In essence, the transition zones played an important role in designing the winding pattern for a vessel, especially in the solid rocket motor case field with a big polar pole ratio. For example, the ideal winding angle of a helicoidally geodesic path over the cylindrical body is 55 o . However, based on the requirement of continuity in the filament winding process, 55 o could not be obtained. In reality, the solid rocket motor case usually has too much different unequal openings due to the requirement of the general structure design, so the non-geodesic transition between trajectories along the cylinder and the dome must be considered. 24 Therefore, the more different the two polar poles, the more difficult the winding pattern design is. 21 The purpose of the present article is to develop a digital method to design the winding patterns of a filament-wound composite pressure vessel under the limitation of unequal polar openings and the requirement of continuity in winding. Furthermore, to obtain a suitable bandwidth in winding, we proposed a model to find the relationship among the winding patterns, bandwidth, and the number of tangent points. This article is arranged as follows: in section "Differential equations for winding angle and rotating angle of mandrel," we reviewed differential equations of non-geodesics, which are the fundamentals of the winding pattern design. In section "Winding pattern design for composite pressure vessel," we demonstrated the detailed method to design the winding pattern of a vessel with unequal polar poles. In section "Computer implementation of the designed winding pattern," the detailed procedures of computer implementation of the designed winding pattern were given to verify the feasibility of the designed winding pattern. In sections "Results and discussion" and "Conclusion," some results and discussions were conducted.

Winding angle equation
where r is the radial coordinate, z the axial coordinate, α the winding angle, and λ the slippage coefficient, respectively, as shown in Figure 1. When λ = 0, equation (1) becomes the well-known Clairaut equation. In the other case (λ ≠ 0), the fiber path will deviate from the geodesic, which is termed as the non-geodesic winding. One should note that there is no analytical solution of equation (1) for λ = 0. Therefore, the winding angle along the axial direction could be computed by the Runge-Kutta method based on some constraint conditions.
The rotation angle of mandrel Figure 1 shows the geometrical relationship of filament trajectory on a revolved surface. A-P-B is the fiber path, R cylindrical radius, and θ angular coordinate in the circumferential direction. ϕ is the angle between radial direction and the curvature radius direction at point P. The relationship between ds and dϕ is where r 1 is the curvature of radius at point P. According to the geometrical relationship of dr and dϕ depicted in Figure 1(b), we have By the geometrical relationship of d s , ds and rdθ depicted in Figure 1 Using the trigonometric function relationship, the cosine value of angle ϕ can be derived as  where r′ is the first derivative of r with respect to z. By substituting equation (5) into equation (4), the first-order differential equation of angular coordinate θ is obtained as d dr Simultaneous solution of the system of differential equations (1) and (6) will finally yield the winding angle and angular coordinate distribution along the axial direction, so the fiber trajectory on the revolution surface will be determined using the angular coordinate.
The filament winding process requires perfect tangent placement of the roving when passing the polar pole of the dome to make filament continuous to a next wound circuit. Therefore, the winding angle at the polar opening should equal 90°. It should be emphasized that a slightly reduced value for α at tangent point of polar pole is desirable (herein α π = 0 4999 . ) to avoid singularity in the numerical solution procedure of equation (1). The slippage coefficient belongs to a feasible interval − { } µ µ , where µ is the friction coefficient between the filament bundle and its supporting surface.

Winding pattern design for composite pressure vessel
To design the winding pattern for a vessel, it is necessary to study the pathways for both the dome and cylindrical sections. This section will give the distribution of the fiber paths for different λ on its supporting surface such as dome surface and cylinder surface, then draw forth how to design a winding pattern for a whole vessel. 3. This range of lambda is selected to demonstrate the effect of slippage coefficient on the distribution of fiber trajectories and winding angle curves of the dome section. Furthermore, filament trajectories could be determined using the development method according to the given λ. Thus, it is easy to change λ to design a winding pattern while the mandrel profile is given. While the axial coordinate z is fixed, the winding angle decreases with decrease in λ. When the λ is negative, the winding angle decreases monotonously with the increase in axial coordinate z. By contrast, while the λ becomes positive, with increase of z, the development of the winding angle along the axial direction can be split into three sections: quick section, slow down section, and slow up section. There is a remarkable difference between the maximum and the minimum winding angles in the domecylinder conjunction corresponding to λ = 0 3

Fiber pathways on dome section
. and λ = −0 3 . , and this provides the scope for winding pattern of a vessel. Figure 2(b) plots the fiber trajectories distribution with different dome heights and slippage coefficients on an ellipsoidal dome. To compare the distribution of these pathways, the initial points that are located at the domecylinder conjunction of all these trajectories are set to the same point. For λ = 0, the paths are geodesics and for λ ≠ 0 the paths stand for non-geodesics. Fiber trajectories corresponding to different slippage coefficient signs are distributed on both sides of geodesics. The deviation of the trajectory from the geodesics increases with the absolute value of the slippage coefficient of the non-geodesics.

Fiber pathways for cylinder section
For the cylinder section, considering the meridian equation is r R = , equations (1) and (6) can be rewritten as where α 0 and α are the initial winding angle and termination winding angle for cylindrical non-geodesics, and L 0 and L are the initial axial coordinate and final axial coordinate. By solving equation (7), α can be obtained, and θ can be determined by equation (8). Here, we acquire the pathway coordinate data. Figure 3 shows the fiber trajectories on the cylindrical surface with various λ . This cylinder radius and height are 100 and 300 mm, respectively. Figure 3(a) and (b) shows positive λ paths and negative λ paths. For comparative analysis of these trajectories, the start points of these paths are fixed at one point and their start angles are fixed as 60 o . In Figure 3, positive λ means the increase in winding angle, and negative λ means a decrease. Different λ at the cylindrical section determines its winding pattern and its winding angle change, while the initial point is selected. This law provides the possibility of designing various winding patterns for a vessel. Figure 4(a) illustrates the winding pattern scheme for a vessel with unequal openings. The vessel is divided into three parts: left dome, cylinder, and right dome. Each part has two possibilities: geodesic and non-geodesic. Owing to this division, there are eight winding schemes for the vessel. Previously, the winding pattern design for a vessel usually adopted a constant winding angle for the cylindrical section. Nevertheless, this is not always feasible, especially when the radius of the polar pole is much bigger than another one. Consequently, it is necessary to employ variational winding angle for the cylinder section to design the winding pattern for a vessel with two unequal openings, and this scenario is more complicated.

Winding pattern design for a vessel with unequal polar openings
There are five patterns for the cylindrical section, as shown in Figure 4(b). The five scenarios are (1) left nongeodesic, middle geodesic, and right non-geodesic; (2) left non-geodesic and middle geodesic; (3) left geodesic and right non-geodesic; (4) geodesic for entire cylinder section; and (5) non-geodesic for whole cylinder section.
To ensure fiber stability and full coverage on the mandrel, four conditions must be satisfied, as follows: (1) all tangent points of a circuit evenly divide the cross-section; (2) the circumferential distance between adjacent rovings at a proper situation is the bandwidth; (3) the winding angle is contiguous; and (4) the slippage coefficient λ is less than the static coefficient of friction µ.
The dome end angle interval is defined as the winding angle interval at the dome-cylinder conjunction when the λ ranges from -µ to µ. In terms of the intersection between the left dome end angle interval and the right one, there are two cases: (1) intersection for the two intervals and (2) no intersection for the two intervals. Figure 5 illustrates the scheme of winding pattern design of both the dome end angle intervals having no intersection. x 1 represents the x-axis of the left dome and is the axial distance between a point of the left dome and the left polar pole. x 2 represents the x-axis of the right dome and is the axial distance between a point of the right dome and the right polar pole. Five sections are employed to plan the winding patterns. Five parameters α α α  of the cylinder which is designated by the user. H 1 and H 2 denote the non-geodesic region lengths for both sides of the cylinder section. α 2 represents the winding angles at the left dome-cylinder conjunction. α 3 represents winding angles at the right dome-cylinder conjunction.
For convenience, let α LeftMax and α LeftMin stand for the maximum and minimum winding angles for the left sides of cylinder section, respectively. Let α RightMax and α RightMin be the maximum and minimum winding angles for the right sides of cylinder section, respectively. In sum, the  procedure of the winding pattern design for both sides of the dome end angle intervals having no intersection is given as follows: = . The computational method of H 2 and α 3 is identical to the one when α 1 is in the first design section.

α 1 is in the third design section
In this case, α 1 is restrained to an interval α α = . The computational method of H 1 and α 2 is identical to that when α 1 is in the third design section. ) for this case and both sides of the cylinder section have a non-geodesic area. A similar process related to α 1 in the third design section is carried out to calculate H 1 , α 2 and H 2 , α 3 . In order for better understanding, Table 1 gives the condition summary when both the dome end angle intervals have no intersection.
For the case when both sides of the dome end angle intervals have no intersection and α RightMin is less than zero, all the four design sections are considered separately to calculate the winding patterns. The computing method of parameters α 1

< < H H mid
When both sides of dome end angle intervals have an intersection as shown in Figure 6, a similar process is carried out to design the winding pattern. Under this scenario, the third design section has a constant winding angle in its cylindrical section. But in any other design sections, it is unable to use the geodesic pattern in the cylindrical section to wind. To satisfy the non-slipping condition and winding angle continuum condition, it is necessary to adopt the proposed method of winding pattern design.

Computer implementation of the designed winding pattern
Before filament winding, simulations prove a tool to verify whether the predesigned path is appropriate at a low cost. Here, we show the simulated winding process using the computer. The proposed algorithm for the filament winding pattern simulation was implemented using the MATLAB software platform in a Windows 7 environment. To avoid producing a gap that leaves the mandrel exposed, the method of full coverage of filament band on mandrel needs to be developed.

Searching slippage coefficients
A key factor in this design is the determination of a pair of slippage coefficients for both the dome parts, to provide non-geodesic winding with unequal polar openings. The filament winding process requires perfect tangent placement of the rovings when passing the polar areas of the dome structure, to make it continuous to the next wound circuit. As a result, by application of the non-geodesic roving trajectories, the requirement must be met. The winding angle at the polar area should reach 90 o . Moreover, to ensure the continuous variation of winding angle when passing the dome-cylinder conjunction, the winding angle of dome and cylinder section at the dome-cylinder conjunction must have the same value. Thus the following constraints at the location of the dome-cylinder conjunction can be given as To prevent slipping of the rovings on the supporting surface, the slippage coefficient λ should meet the stability condition as follows where µ is the friction coefficient between the fiber bundle and the supporting surface. Here the slippage coefficient is given by a piecewise function

Left dome
Right dome (18) The target is to find a pair of slippage coefficients ( , ) λ λ

The process of computer implementation
During the filament winding process, the fibers are wound on the mandrel surface in a continuous manner and are tangent to the opening at its polar radius. Therefore, the mandrel angular θ monotonically increases from zero to a value when the winding is over. The winding coordinate system is defined as shown in Figure 8. The fiber trajectories are created using its Cartesian coordinates. The Cartesian coordinates of the fiber paths can be given as where r z ( ) is the radius of the mandrel at the location on z-axis, and α is the winding angle.
To visualize the designed winding patterns, the process of fiber paths implementation for a wound cycle is separated into six portions: cylindrical winding section, right dome up section, right dome down section, cylindrical back winding section, left dome up section, and left dome down section, as shown in Figure 8. After a cycle winding, the mandrel angular θ is updated by the end value of the previous cycle winding. Repeat the cycle winding until a certain number, then the mandrel could be in full coverage. The detailed description of the method of computer implementation for the designed winding pattern is depicted in Figure 9.

Full coverage
To ensure the mandrel is uniformly covered by fibers, the number of wound circuits (a wound circuit is defined as the winding process in which continuous fibers start at a certain point on the mandrel surface and wind several cycles before returning to the starting point.) must be acquired. In a filament winding system, fiber tows under tension are wound onto the mandrel surface side by side in designed patterns to form a uniform composite shell that is precisely manipulated by computer numerical control. Consequently, the circumferential distance between each circuit path at the proper situation should be a bandwidth to avoid producing a gap that leaves the mandrel exposed. The distance is defined by a pitch angle δ as where b is the width of the filament band, D is the diameter of the mandrel at the location on the z-axis, and α ( ) z is the winding angle. While getting the minimum pitch angle for the whole mandrel, the number of the wound circuit and cycle of the mandrel is uniformly covered with fibers as Ci r z b n Cy n ceil Ci where n is the number of tangents; and Ci and Cy are the number of circuits and cycles when the mandrel is uniformly covered by filaments. Ci is a real number but Cy is an integer number. The circumferential width of the individual filament band placed on the mandrel is not necessarily an integer multiple of the cross-sectional circumference of the mandrel. As a result, the overlap of filament bands may exist in the last wound circuit, as shown in Figure 10. It demonstrates the process of a mandrel uniformly covered by fibers, which includes the first circuit, the sixth circuit, and the final circuit winding, and the results show that the mandrel can be coated uniformly with filament bands, which indicates the effectiveness of the proposed method.
In order to quantify the cover factor and overlap of filament band, C f and Lap are defined as cover factor and overlap factor of filament band. C f represents the cover degree of filament band after a layer of finished winding. Combining equation (21), cover factor C f is expressed as Combining equation (22), the overlap factor of filament band is obtained as The relationship among the number of tangent points, roving bandwidth, and winding pattern Once the mandrel geometry and winding parameters of composite vessel are determined, the winding pattern is simultaneously given and the critical roving width is a desirable parameter for manufacturing the vessel. When the real bandwidth is less than the critical roving width, a gap will appear between two adjacent rovings in some area of the vessel. By contrast, while the real bandwidth is greater than the critical roving width, the overlap of rovings exists between adjacent yarns in all areas of the wound product. Here, let θ n denote the practical angular of the mandrel for a cycle of pay-out eye. The parameter θ n ′ is defined the angle of a mandrel for a pay-out eye cycle without considering the bandwidth. The parameter ∆θ represents the tiny angle of a mandrel for a pay-out eye cycle caused by a bandwidth. The relationship for θ n , θ n ′ , and ∆θ is given as For a mandrel, the angle of the mandrel for a pay-out eye cycle is a constant while the winding pattern is determined. The critical bandwidth is related to the angle θ n ′ , tiny angle ∆θ , and the number of tangent points n. The number of tangent point n is a key parameter. Once the number of tangent points is determined, the angle θ n ′ and ∆θ then are two constants. Hence, there are two uncorrelated parameters in the angle θ n ′ , the angle ∆θ , and n. The constraint condition that angle ∆θ is limited to a predetermined value is adopted to calculate θ n ′ , ∆θ, and n. Figure  11(a) depicts a cross-section for a wound circuit for θ n = 453 74 . o and ∆θ < 15 o . Small circles in Figure 11(a) and (b) stand for a roving and the number represents the order of roving wound on the mandrel. Figure 11(b) shows the cross-section of a wound circuit for θ n = 453 74 . o and ∆θ < 10 o . It is observed from Figure 11 that the angle between the start point and the nth point equals that between the first point and the (n + 1)th point. In Figure  11(a), the number of tangent points is 4 for ∆θ = 14 9 . o . In Figure 11(b), the number of tangent points is 23 for ∆θ = 3 89 . o . The smaller the angle ∆θ , the greater the number of tangent points is. The detailed flowchart of calculating bandwidth and the number of tangents is given in Figure 12.

Results and discussion
The simulation system has the ability of designing the winding pattern for a vessel with unequal polar openings.
The simulation system has the ability of rotation, panning, zooming, and view switching, which enable the designer to observe the winding process by dragging the mouse. If there are any abnormal phenomena during the filament winding process, the engineer can stop and adjust the winding pattern until the filament bands cover the mandrel uniformly. Therefore, the simulation is convenient and effective in producing a composite vessel.  To validate the ability of generating winding patterns of the developed simulation system for vessels with different ratio of polar poles, Figure 13 displays non-geodesic winding patterns for revolution shells with different ratios of polar openings for 1 2 1 4 :~: . As followed from this picture, for each ratio of polar pole, the proposed method always has the ability to cover the mandrel uniformly with fibers even if the ratio of polar poles is greater at 1:4. Figure 13 also indicates that as the number of the tangent point increases, the number of the wound circuit decreases after the mandrel is coated with an individual layer. When the tangent number of points is 3, 5, and 6, respectively, the wound circuits are 37, 21, and 19, respectively. There may have been overlap of filament bands in the last wound circuit, as observed in Figure 13. This is caused when the circumferential width of the individual filament band placed on the mandrel is not an integer multiple of the cross-sectional circumference of the mandrel, and for this scenario the overlap is unavoidable.
For the sake of displaying winding patterns for different design sections for a certain vessel, a revolution shell is selected and the geometrical parameters are depicted in Figure 14. The semi-minor axis and the semi-major axis are 330 and 385 mm, respectively, fiber bandwidth b = 20mm,  thickness of a filament band t = 0 2 . mm, a helical geodesic winding, and the maximum friction coefficient µ is set as 0.3. The calculation result indicates that for this shell there are four design sections of winding pattern. Figure 15 gives the winding angle curves along the axial direction for four design sections of winding pattern of the vessel. When α 1 is in the first design section, as the axial coordinate z increases the winding angle of the left dome decreases and then increases. In its cylindrical section, both the sides have a non-geodesic region. The winding angle increases to 60 o , then keeps this value and finally decreases. As for the right dome, the winding angle decreases, then increases to 90 o , which satisfies the requirement of being tangent at its polar radius. When α 1 is in the second design section, the variations of winding angle in the left and right domes are similar to those when α 1 is in the first design section. In the cylindrical section for this case, non-geodesic exists only on the right side of the cylinder section, and for any other position the winding angle is a constant 40 o . While α 1 is in the third design section, for the left dome section, the winding angle decreases monotonously. No non-geodesic region exists in its cylindrical section, in other words, it is geodesic winding for the cylindrical section. For the right dome, the winding angle decreases, and then increases to 90 o . In the last case while α 1 is in the fourth design section, the winding angle decreases in the left dome. The non-geodesic winding area lies in the left part of the cylindrical section. As can be seen from the figure, the winding angle is continuous along the axial direction and the slippage coefficient at any position of the shell is less than the predetermined value 0.3, which accounts for the feasibility of the proposed method of winding pattern design.  Table 2 shows that numerous tangent points and bandwidths for each design section of winding pattern could be selected. For the first and second design sections, as the tiny angle ∆θ is a certain value, so the number of tangents and critical bandwidth are invariable for different ∆θ max . By contrast, for the third and fourth design sections, ∆θ decreases as ∆θ max decreases, thus the number of tangents and critical bandwidth are variable. As the tiny angle ∆θ decreases, the number of tangent point increases and the critical bandwidth decreases.
In other words, the more the number of tangent points, the more the number of wound circuits when the mandrel is covered uniformly with a single layer. In general, it was found that composite vessel strength was significantly affected by the laminate stacking sequence, winding angle (associated with winding pattern), winding tension, winding-tension gradient, and winding time. 31,32 Thus, it is tough work to find the optimal winding pattern under the given vessel geometry and this will be the future work.  n is the number of tangent point, b is the critical bandwidth, ∆θ represents the tiny angle of a mandrel caused by a bandwidth, ∆θ max is the maximum limitation of ∆θ .

Conclusion
In this article, the method of winding pattern design of filament wound composite pressure vessel with unequal openings was proposed using geodesics and non-geodesics, which also satisfy the winding principles. A model was established to determine the relationship among the number of tangent points, bandwidth, and winding pattern, which aim to ensure the fiber stability and full coverage on the mandrel to satisfy the winding process requirements. The winding patterns for composite pressure vessel were simulated on the MATLAB software platform to verify the feasibility of the trajectory deign. The simulation results indicate that the mandrel was coated with fibers uniformly and no severe overlap regions except the polar pole regions exist on the mandrel surface, which proved the effectiveness of the method. The method introduced in this article can be used in the early stages of designing a composite pressure vessel, especially in a solid fuel rocket case.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.