To simplify the process of modeling the 3D path of a wound yarn, a reasonable procedure is to first calculate the trajectory of the yarn centerline and then assign a yarn width and thickness to the centerline. The process of calculating the yarn centerline involves calculating the geodesic and semi-geodesic on the mandrel surface, as well as the winding pattern generation. As shown in Figure 2A, the yarn centerline consists of three kinds of trajectories, which are circular arc, helix, and semi-geodesic. When the cylindrical surface is expanded into a rectangle (Figure 2B), both the arcs and the helixes are shown as straight lines, for they are both geodesics (the geodesic is a straight line after the developable surface is spread). In contrast, the semi-geodesics are shown as curves, which means that its winding angle, i.e., the angle from the meridional direction to the yarn direction (Koussios, 2004), varies according to the z-coordinate.

A parametric representation of the helix and circle can be made by the following: { x = R c o s ( ϕ ) , y = R s i n ( ϕ ) , z = C ϕ           h e l i x x = R c o s ( ϕ ) , y = R s i n ( ϕ ) , z = C             c i r c l e (1) where R is the cylinder radius, φ is the rotation angle from the x-axis, and C is a constant. The winding angle, α, can be given by the following: { α = a r c t a n ( R / C )               h e l i x α = π / 2                                             c i r c l e (2) The semi-geodesic has a variable winding angle and is used to reverse yarn motion direction at both ends of the mandrel. However, compared with the geodesic, the semi-geodesic trajectory can lead to yarn slippage on the mandrel. In order to maintain the mechanical stability of the yarn, it is essential to ensure that the ratio of the frictional force and the normal pressure acting on the yarn is less than the friction coefficient μ between the yarn and the mandrel. This equilibrium condition could be expressed by differential geometric quantities, geodesic curvature, and the normal curvature, which can be written as follows (De Carvalho et al., 1995; Koussios, 2004; Koussios et al., 2004; Kim et al., 2005; Vargas Rojas et al., 2014): d α d z = η sin 2 α R cos α (3) where η is slip coefficient, it must satisfy 0 ⩽ η ⩽ μ .

Through integral ∫ z 0 z d t = ∫ α 0 α R cos α d α η sin 2 α , the analytical solution of Eq. 3 can be written as follows: z − z 0 = R η ( 1 sin α 0 − 1 sin α ) (4) where α 0 and z 0 are the winding angle and axial position of the start point, respectively.

In order to maintain a constant length increment ds along the semi-geodesic, the step size of z is calculated in the following way: Δ z = z i + 1 − z i = d s ⋅ c o s ( α ) (5) The relationship between rotation angle ϕ i + 1 and ϕ i can be expressed as follows: ϕ i + 1 = ϕ i + ∫ z i z i + 1 tan [ α ( z ) ] R d z (6) To have the yarn cover the mandrel uniformly, the total rotation angle experienced during a cycle must be designed according to the particular winding pattern. The winding pattern is formed by the yarn going back and forth along the mandrel and crossing and overlapping itself at certain points. Once the helix and semi-geodesic trajectories are determined, the yarn usually dwells in both ends of the mandrel when α = π/2 for a certain angle of Δ ϕ d w e l l , forming a circular arc. Δ ϕ d w e l l can be expressed as follows: { Δ ϕ d w e l l = π − 1 2 m o d ( Δ ϕ n − 2 π d N ,   2 π ) N = c e i l ( 2 π R b cos ( α ) ) (7) where Δ ϕ n is the rotation angle experienced by helixes and semi-geodesics in a cycle; N is the number of cycles required to completely cover the mandrel; b is yarn width; d governs the pattern formed during the winding (Johansen et al., 1998), it must satisfy 1 ⩽ d < N , and the ratio d / N is non-reducible. After N cycles, the yarn will eventually cover the mandrel (Figure 3).

Many different yarn cross-sections have been presented over the years to indicate the actual condition of the yarn before and after tension is applied. These studies have mainly been carried out for representative volume elements (RVE) or advanced fiber placement (AFP) techniques (Potluri and Manan, 2007; Ansar et al., 2011; Kim et al., 2015b; Vernet and Trochu, 2016). Commonly used yarn cross-section shapes include rectangular, ellipsoidal, lenticular, and so on. However, in the filament winding process, the yarn is expanded into a flat shape with a rectangular cross-section. Hence, our work considers the yarn cross-section as a rectangular cross-section, and the shape remains consistent throughout the winding process.

As shown in Figure 4A, based on the centerline of the yarn path T ( 0 ) , multiple parallel trajectories are obtained by geodesic offsets on the mandrel surface to describe the yarn width (Fu et al., 2016). The multiple parallel trajectories are represented by a set of B = { T ( i ) ∣ i = 0,1 , ... , 2 m } , where m is the number of equal parts of the half yarn width. Each trajectory T ( i ) = { P ( i , j ) ∣ j = 0,1 , ... , n } is composed of a sequence of discrete points of P ( i , j ) , where n is the number of all discrete points along the centerline T ( 0 ) . All parallel trajectories are then raised by one yarn thickness along the outer normal of the mandrel surface to describe the yarn thickness (Figure 4B), so the number of trajectories in set B is doubled to 4m + 2.

On the cylindrical surface, the geodesic offsets are easy to calculate. T ( i ) = { i = 1 , ... , 4 m + 1 } can be obtained from T ( 0 ) according to the winding angle, the distance from the mandrel axis, the width, and thickness of the yarn. As shown in Figure 4C, assuming the cylindrical coordinate of a point P ( 0 , j ) on T ( 0 ) is ( ρ ( 0 , j ) ,   φ ( 0 , j ) ,   z ( 0 , j ) ) , the variation of the rotation angle φ and the axial coordinate z corresponding to the half yarn width can be written as Δ φ = w ⋅ cos α 2 ρ and Δ z = w ⋅ sin α 2 , respectively, where w is the yarn width. The points P ( 2 i − 1 , j ) ( i = 1 , ... , m − 1 ) are uniformly distributed over the width between P ( 2 m − 1 , j ) and P ( 0 , j ) , its cylindrical coordinates can be expressed as follows: P ( 2 i − 1 , j ) = ( ρ ( 2i−1 , j ) ,   φ ( 0 , j ) − Δ ϕ i m ,   z ( 0 , j ) − Δ z i m ) (8) Points P ( 2 i + 2 m , j ) ( i = 1 , ... , m ) have the same ϕ and z with P ( 2 i − 1 , j ) ( i = 1 , ... , m ) . Compared with point P ( 2 i − 1 , j ) , point P ( 2 i + 2 m , j ) is t further from the mandrel surface, where t is the yarn thickness. The cylindrical coordinates of P ( 2 i + 2 m , j ) can be expressed as follows: P ( 2 i + 2 m , j ) = ( ρ ( 0 , j ) + t ,   φ ( 0 , j ) − Δ ϕ i m ,   z ( 0 , j ) − Δ z i m ) (9) Similarly, for points P ( 2 i , j ) and P ( 2 i + 2 m + 1 , j ) , we can get the following: { P ( 2 i , j ) = ( ρ ( 0 , j ) ,   φ ( 0 , j ) + Δ ϕ i m ,   z ( 0 , j ) + Δ z i m ) P ( 2 i + 2 m + 1 , j ) = ( ρ ( 0 , j ) + t ,   φ ( 0 , j ) + Δ ϕ i m ,   z ( 0 , j ) + Δ z i m ) (10)

As soon as the yarn crosses another part of itself in the winding process, overlap will occur. Assuming that the yarn is always in contact with the mandrel surface or a previous part of itself, a point on the current yarn path obtains a new radial distance ρ n e w = ρ + t in the cylindrical coordinate whenever the current part of the yarn encounters a previous part of itself. Then, if the current part of the yarn crosses with a new part of the yarn later, the radial distance in cylindrical coordinate continues to increase, i.e., ρ n e w = ρ + 2 t .

Figure 5 shows a schematic representation of the overlap processing of two yarns (yarn A and yarn B). Yarn A can be regarded as the part of the yarn wound first and yarn B as the part of the yarn wound later. If there is no overlap processing, yarn B interferes with yarn A at the intersection (Figure 5A), the paths of yarn A and B are on the original mandrel surface. It fails to represent the yarn accumulation on the mandrel surface. After eliminating interference, the 3D path of yarn B is shown in Figure 5B. Intuitively, it is clear that the behavior of yarn B does not conform to physical reality; its distance from the mandrel surface changes abruptly before and after crossing yarn A and therefore, its path appears discontinuous. The underlying reason for the discontinuity is that the radial distance of points on yarn B in the cylindrical coordinate is updated by integer multiples of the yarn thickness during the overlap process.

During the winding process, the yarn is subjected to winding tension. As a result, the tensioned yarn B temporarily leaves the mandrel surface before and after crossing the underlying yarn A. In a stable winding process, the mechanical behavior of the yarn should conform to the principle of minimum potential energy. It is important to note that the yarn sliding in the tangential direction receives a tension constraint and can be considered as having no tendency to slip (the friction force in this direction is zero), so the tension of the yarn can be considered as equal everywhere. Then the elongation of the yarn, ε , can be defined as follows: ε = F E A (11) where F, E, and A denote the tension, the modulus in the length direction of the yarn, and the cross-section area of the yarn, respectively.

Ignoring first the cross-section of the yarn and denoting P ( i , j ) ( i = 0,1 , ... , n ) by P j , the deformation energy, U , of yarn B, can be expressed as follows: U = 1 2 E A ε ∑ j = 0 n − 1 | P j P j + 1 | (12) For the displacement constraint of P j ( j = 1 , ... , n − 1 ) , its coordinates φ and z are obtained based on the previous calculation of the centerline trajectory and cannot be changed. Otherwise, it deviates from the intended winding path. And the coordinate ρ is obtained based on a simplified overlap processing. In fact, the current value of ρ is the minimum value that guarantees the yarn B will not interfere geometrically with the mandrel and yarn A. It does not ensure that the 3D yarn path is physically real, and the actual ρ can be larger than that value. In a word, yarn B acts as a discrete connected system, P j can only move to the outer normal direction of the mandrel, and its final radial coordinate ρ ˜ j must satisfy ρ ˜ j ⩾ ρ j . In the above displacement adjustment process, the support reaction force does not do any work. Hence, the total potential energy of the whole system is equal to the total elastic energy and it can be given as follows: Π = 1 2 E A ε ∑ j = 0 n − 1 | P j P j + 1 | (13) where Π denotes the total potential energy of yarn B.

According to the principle of minimum potential energy (Chandrupatla et al., 2002), under the displacement constraint, the potential energy can be used as the objective function, and the discrete connected system of yarn B will remain in mechanical equilibrium only when the potential energy of the whole system is minimized. Considering that the coefficient, E, A and ε are constants, so the objective function can be simplified as follows: Q ( ρ ˜ j ) = ∑ j = 0 n − 1 | P j P j + 1 | (14) The physical meaning of Eq. 14 is that the equilibrium state of the yarn B, determined according to the principle of minimum potential energy, is in fact, also intuitively the shape of the shortest path of the yarn B.

It can be proven that the shortest total path is equivalent to the smallest sum of squared path lengths under the current constraint (the proof process is omitted here). Therefore, the objective function can be transformed into the following: Q ˜ ( ρ ˜ j ) = ∑ j = 0 n − 1 | P j P j + 1 | 2 (15) The new objective function, Q ˜ , can be understood as the sum of the elastic energy of a zero-length spring with a stiffness of 2 between P j and P j + 1 . It is more suitable for optimal solving than Q because the expression does not contain a root sign, which is beneficial to improve the computational efficiency.

Yarn B has a rectangular cross-section in practical optimization, and its 3D path consists of many trajectories T ( i ) = { i = 0 , ... , 4 m + 1 } . Take, for example, one of these trajectories, T ( i ) ( 0 ⩽ i ⩽ 2 m ) , which is on the lower surface of yarn B and composed of ordered points P ( i , j ) ( i = 0,1 , ... , n ) with cylindrical coordinates of ( ρ ( i , j ) ,   φ ( i , j ) ,   z ( i , j ) ) . It is assumed that the points P ( i , j ) ( j = 1,2 , … , 2 m − 1 ) will be raised to position P ˜ ( i , j ) , and the cylindrical coordinates of P ˜ ( i , j ) can be written as ( ρ ( i , j + 1 ) + δ ( i , j + 1 ) 2 ,   φ ( i , j + 1 ) ,   z ( i , j + 1 ) ) , which naturally satisfy the displacement constraint.

The square of the length of the line segment P ˜ ( i , j ) P ˜ ( i , j + 1 ) can be written as follows: | P ˜ ( i , j ) P ˜ ( i , j + 1 ) | 2 = ( ( ρ ( i , j ) + δ ( i , j ) 2 ) ⋅ cos ( φ ( i , j ) ) − ( ρ ( i , j + 1 ) + δ ( i , j + 1 ) 2 ) ⋅ cos ( φ ( i , j + 1 ) ) ) 2 + ( ( ρ ( i , j ) + δ ( i , j ) 2 ) ⋅ sin ( φ ( i , j ) ) − ( ρ ( i , j + 1 ) + δ ( i , j + 1 ) 2 ) ⋅ sin ( φ ( i , j + 1 ) ) ) 2 + ( z ( i , j ) − z ( i , j + 1 ) ) 2 (16) The objective function Q ˜ ( ρ ˜ j ) , which becomes a quadratic polynomial in the variable δ ( i , j ) ( j = 1,2 , … , 2 m − 1 ) , which can be written as follows: Q ˜ ( δ ( i , j ) ) = ∑ j = 0 n − 1 | P ˜ ( i , j ) P ˜ ( i , j + 1 ) | 2 (17) The value of δ ( i , j ) corresponding to the minimum value of this objective function is the solution to this optimization problem. We use the L-BFGS (Nocedal, 1980) method to solve this optimization problem, which is the most common method for solving unconstrained nonlinear programming problems and has the advantages of fast convergence and low memory overhead often found in various machine learning algorithms. In order to apply the L-BGFS algorithm, the gradient function of the objective function Q ˜ ( δ ( i , j ) ) with respect to the variable δ ( i , j ) must be provided. ∂ Q ˜ ∂ δ ( i , j ) = 4 [ 2 ( ρ ( i , j ) + δ ( i , j ) 2 ) − ( ρ ( i , j − 1 ) + δ ( i , j − 1 ) 2 ) cos ( φ ( i , j − 1 ) − φ ( i , j ) ) − ( ρ ( i , j + 1 ) + δ ( i , j + 1 ) 2 ) cos ( φ ( i , j + 1 ) − φ ( i , j ) ) ] δ ( i , j ) (18)

The optimized 3D path of the yarn is shown in Figure 5C. It can be seen that the shape of the yarn B becomes smooth and leaves the surface of the mandrel before and after crossing yarn A, conforming to physical reality.

As shown in Figure 6, the 3D path modeling strategy for wound yarn was applied to the fuel cladding tube. The main geometrical and processing parameters are listed below:

• The width and thickness of the SiC_f yarn are 1.1 and 0.1 mm, respectively;

The simulated 3D path of yarn shows the “bamboo knot” structure of the wound part (Figure 6B). The winding path in the turnaround zone is semi-geodesic and the winding angle changes from the 45°to 90°, where yarn enrichment occurs and the thickness is, therefore, thicker than that in the other areas, which is also consistent with reality.

As seen in Figure 7A, the tube has a length of 100 mm, an inner diameter of 20 mm, and a wall thickness of 0.5 mm. The wound yarns have 153,504 C3D8R elements, and the matrix has 77,952 elements (Figure 7B), while the laminated tube has 103,936 elements (Figure 7C). The global fiber volume fraction is 28% for both wound and laminated tubes. The input material properties are listed in Table 1, calculated from the micro-UD model (Feng et al., 2020). The SiC matrix has Young’s modulus of 300 GPa and a Poisson ratio of 0.14.

The embedded element technique is used to specify that an element or group of elements is embedded in “host” elements. As shown in Figure 8A, the matrix mesh is the host mesh, and the yarn mesh is the embedded mesh. Abaqus searches for the geometric relationships between nodes of the embedded elements and the host elements. If a node of an embedded element lies within a host element, the translational degrees of freedom at the node are eliminated and the node becomes an “embedded node.” The translational degrees of freedom of the embedded node are constrained to the interpolated values of the corresponding degrees of freedom of the host element (Tabatabaei et al., 2014). Figure 8B shows the displacement boundary conditions at both ends of the tube. The load is applied directly to the master node at the end of the tube. Then the master node transfers the load to the slave nodes at both ends of the tube through continuum distributing constraints by weight. Take the axial tension condition as an example, the three translational degrees of freedom of the left master node are set to zero, and the tension load is applied through the axial displacement of the right master node.

With the obtained mesh (Figure 7B) and material properties (Table 1), the meso-scale FE model of the cladding tube can be applied to structural optimization. Generally, the model of the laminated tube with identical fiber orientations are often used to predict the major mechanical properties, e.g., apparent modulus and strength, as a simplified substitute. However, the weakening effect of the yarn waviness and interlocks caused by the overlap effect has not yet been specifically investigated in the previous work of cladding tube design. The mechanical stresses on the cladding tubes are small, and the main risk is the mechanical failure during a sudden axial movement (Sauder, 2014), which will result in axial tensile or compressive stresses. In addition, the weakness of the SiC_f/SiC composite is the low shear strength, which is reflected in the axial tensile condition where the fiber direction is exactly the maximum shear stress. Therefore, in the following discussion, an equivalent tensile loading will be applied to the tube; the apparent moduli of the wound tube and laminated tube will be calculated and compared.

Figure 9A shows that the global axial tensile moduli of the wound and laminated tubes are quasi equivalent. The global tensile modulus of the laminated tube is 322.0 GPa, and the global tensile modulus of the wound tube is 322.2 GPa. The deviation is negligible. In such a circumstance, the meso-scale FE model can be used to predict the apparent moduli of the cladding tube, especially under complex service loadings.

Although the calculated moduli of the wound tube and the laminated tube present good agreement, the local stress states are totally different since the waviness and interlock of the yarn cause stress concentration, which will introduce an earlier damage initiation and final failure of the cladding tube. Shear stress is the typical stress state of the tube under axial tension and the shear stress safety factor is demonstrated in Figure 9B. The safety factor is defined as the ratio of the shear stress and shear strength (the shear strengths of SiC_f/SiC are 60.1 and 78.8 MPa when the fiber volume fraction is 70 and 28% (Feng et al., 2020)). For the wound tube, the shear stress is extracted from the nodes along the yarn on which the maximum stress is located, and for the laminated tube, the shear stress is extracted from the nodes along the axial direction. According to our assumption, the safety factors in the wound yarns are almost 36.7% lower than that of the laminated tube. The lowest safety factor occurs at the “bamboo knot” in the yarn, which means that the damage initiation will occur here due to stress concentration caused by the overlap effect, which cannot be reflected in the laminated tube model. Figure 10A shows the global shear stress just before the damage of the yarn, and it can be seen that the maximum shear stress in the wound tube yarn is 20% greater than the stress in the laminated tube. At the same time, a higher fiber volume fraction will lead to a reduction in shear strength, so the difference in safety factor between the two models is greater than the difference in shear stress. In such a context, with the same fiber volume fraction, the wound tube and the laminated one have similar values of moduli but quite different apparent strengths. The stress concentration and higher local fiber volume fraction result in a lower strength for wound tubes. This result confirms that a precise wound structure on meso-scale is indispensable for the property prediction of the cladding tube.