Composite pressure vessel filament winding angle sensitivity analysis and optimization method

Abstract

The application provides a composite pressure vessel fiber winding angle sensitivity analysis and optimization method, and belongs to the technical field of composite pressure vessel design and manufacturing. The method comprises the following steps: S1, model partitioning; S2, manufacturing constraint; S3, establishing a constitutive relation; S4, design variable parameterization; S5, mathematical expression of the optimization problem; S6, sensitivity analysis solution; S7, design variable updating; S8, iteration cycle and convergence determination. The method converts the design variables of the optimization problem from the unmanufacturable "local layer angle" to the actual "winding path parameters" that determine the fiber space angle, thereby strictly embedding the manufacturing constraint in the optimization process. Through the construction of a high-fidelity finite element model and the use of efficient analytical sensitivity analysis, efficient exploration of the design space is realized, and finally a winding angle distribution scheme that seamlessly connects with the numerical control winding equipment and can be directly used for multi-beam non-crossing production is output, thereby realizing the technical closed loop from optimization design to manufacturing planning.

Description

Technical Field

[0001] This invention belongs to the technical field of composite material pressure vessel design and manufacturing, and specifically discloses a method for analyzing and optimizing the sensitivity of fiber winding angle in composite material pressure vessels. Background Technology

[0002] Composite material pressure vessels, as high-performance storage and transportation structures, have been widely used in aerospace, rail transportation, clean energy, chemical storage and transportation, fire protection, and medical fields. For example, in the aerospace field, their weight directly affects carrying capacity and mission costs; in hydrogen fuel cell vehicles and compressed natural gas vehicles, their lightweight and high safety requirements directly affect driving range and operating economy; in the chemical and energy fields, they are used to store and transport high-pressure, flammable, or corrosive media, requiring excellent pressure resistance and corrosion resistance.

[0003] Fiber winding is a mainstream process for manufacturing composite pressure vessels. Multi-bundle winding, as an advanced process, significantly improves winding efficiency and material property uniformity by simultaneously controlling the paths and tensions of multiple fiber bundles. In particular, the non-crossing winding mode within the same layer ensures that multiple fiber bundles are arranged in parallel and do not cross each other through precise path planning. This effectively avoids local thickness abrupt changes and stress concentrations caused by fiber crossings, thus potentially further improving the overall structural performance and fatigue life. The winding angle, as a core parameter of the process, directly determines the fiber distribution relative to the principal stress direction, thereby affecting structural stiffness, strength, stress distribution pattern, and fatigue life.

[0004] However, the design of the winding angle is subject to strict process constraints. In multi-bundle winding, to achieve stable, non-crossing layup within the same layer, the spatial trajectory of the fibers on the three-dimensional curved surface must satisfy more complex geometric coordination conditions. Once the core geometry and winding pattern (such as the number of tangent points and winding links) are determined, the spatial angle of each fiber bundle at any point on the shell is locked to a fixed value, rather than the locally adjustable layup angle assumed in conventional optimization. Many existing winding angle optimization design methods are based on the assumption of a single bundle or equivalent single layer, and do not rigorously consider the strong constraint of fixed spatial angles for each fiber bundle, especially in multi-bundle, non-crossing winding within the same layer. They often simplify the winding angle into an independently designable variable. This neglect of actual process constraints results in optimization results that, although numerically optimal, cannot be implemented on multi-bundle non-crossing winding equipment, lacking manufacturability.

[0005] Therefore, a new winding angle optimization method is urgently needed for the advanced manufacturing mode of multi-bundle winding and non-intersecting winding within the same layer. This method needs to be able to strictly embed the core manufacturing constraint of "fixed azimuth angles of multiple fiber bundles" within a high-precision finite element analysis framework, treating complex winding path parameters, rather than fictitious free angles, as design variables. Simultaneously, it must efficiently handle the large number of variables introduced by multi-bundle winding through analytical sensitivity calculations, thereby significantly improving computational efficiency and design accuracy while ensuring that the optimization results can be directly used for multi-bundle non-intersecting winding path planning. This will not only reduce the need for extensive physical prototyping relying on trial and error, but also provide key technical support for the rapid and accurate design of high-performance composite material pressure vessels, shortening the R&D cycle, reducing costs, and overcoming existing performance bottlenecks. Summary of the Invention

[0006] The purpose of this invention is to overcome the limitations of existing optimization methods in handling multi-bundle non-crossing winding process constraints, and to propose a sensitivity analysis and optimization method for fiber winding angle in composite pressure vessels that considers fixed spatial angle constraints. The core concept of this method is to transform the design variables of the optimization problem from unmanufacturable "local layup angles" into actual "winding path parameters" that determine the fiber spatial angles. This allows for the strict embedding of manufacturing constraints during the optimization process. By constructing a high-fidelity finite element model and employing efficient analytical sensitivity analysis, the design space is explored efficiently. The final output is a winding angle distribution scheme that seamlessly integrates with CNC winding equipment and can be directly used for multi-bundle non-crossing production, achieving a technical closed loop from optimization design to manufacturing planning.

[0007] The above-mentioned method for analyzing and optimizing the sensitivity of fiber winding angle in composite pressure vessels includes the following steps:

[0008] S1, Model Partitioning:

[0009] Composite material pressure vessel domain Oh Divided into mutually exclusive subfields:

[0010] (1);

[0011] in, Oh L This indicates the lining area and the end cap area. Oh H Indicates the region of the composite material winding layer;

[0012] Using the finite element method, the composite material winding layer region was analyzed. Oh H Discrete N e A system composed of individual units, targeting Oh H Each unit within e ,e =1,2, ..., N e It imparts homogenized orthotropic anisotropic material properties to arbitrary elements. e The fiber orientation is determined by a pair of angles ( β , i (Together described) β It is the azimuth angle. i For the wrapping angle;

[0013] S2, Manufacturing Constraints:

[0014] Constraints are created by embedding a fixed azimuth angle principle into the finite element model. Under the constraint of the fixed azimuth angle principle, the azimuth angle of each element is... β The winding angle is considered a fixed parameter predetermined by the geometry of the composite pressure vessel, rather than a design variable. i As design variables, the optimization process only occurs when the variables are fixed. β Design variables under distribution i Make adjustments;

[0015] S3, Establish constitutive relations:

[0016] The liner and end cap are made of isotropic materials, and the stress vector is... With strain vector The constitutive equation between them is expressed as:

[0017] (3);

[0018] (4);

[0019] in, D L Here is the stiffness matrix. E For Young's modulus, Poisson's ratio;

[0020] At the microscale, the composite winding layer is equivalent to a macroscopic orthogonal anisotropic material, and a local coordinate system is established with three mutually orthogonal principal directions of the material as coordinate axes. In the local coordinate system, the fiber axis is represented as direction 1, and the fiber radial direction is represented as directions 2 and 3. Assuming that the fiber axis is consistent with the global coordinate system, the stress vector... With strain vector The constitutive relation between them is expressed as:

[0021] (5);

[0022] (6);

[0023] in, It is the elastic matrix, serving as the original constitutive matrix and the basis for subsequent coordinate transformations to any direction;

[0024] ; ; ;

[0025] ; ;

[0026] ;

[0027] ; ; ;

[0028] The Young's modulus in the 1st direction. The Young's modulus in two directions. It represents the Young's modulus in three directions.

[0029] For the Poisson's ratio in the 1-2 plane, For the Poisson's ratio in the 1-3 plane, For the Poisson's ratio in the 2-3 plane,

[0030] The shear modulus of the 1-2 plane. The shear modulus of the 1-3 plane. For the shear modulus of the 2-3 plane, ;

[0031] Constrained by constitutive compatibility relations, according to ( i , j = 1,2,3; i ≠ j ),available , , ;

[0032] S4, Parameterization of Design Variables:

[0033] For the original constitutive matrix Applying a rotational transformation yields a result containing fiber orientation ( β , i elasticity matrix D e :

[0034] (7);

[0035] in, T ( i ') andT ( β ')for i and β The rotation tensor;

[0036] S5, Mathematical formulation of the optimization problem:

[0037] The total strain energy of a composite pressure vessel under a given load. Represented as:

[0038] (13);

[0039] in, J L This represents the strain energy of the lining and end caps, independent of fiber orientation. J H ( i ) indicates the angle of wrapping. i The strain energies of the relevant composite winding layers are as follows:

[0040] (14);

[0041] in, This indicates the strain energy density of the lining and end caps. This represents the strain energy density of the composite winding layer;

[0042] The fiber orientation optimization problem can be expressed in the following mathematical form:

[0043] (15);

[0044] in, Describe the objective function. K ( i , β This indicates dependence on the fiber orientation field. i , β The structural stiffness matrix, u Represents the displacement vector. F Represents the load vector;

[0045] The integral form of the optimized objective function:

[0046] (19);

[0047] in, = Bu For the true strain vector, the direction component is obtained through the rotation tensor. T ( i ') and T ( β ') is directly related to the fiber orientation angle;

[0048] S6, Sensitivity Analysis Solution:

[0049] Due to the azimuth angle of each element after geometric discretization β This has been determined, therefore the sensitivity analysis only considers the winding angle. i The sensitivity was derived using the adjoint variable method, and a virtual strain tensor was introduced under the assumption of material isotropy. Alternative For the winding angle i Perform sensitivity analysis:

[0050] (20);

[0051] The sensitivity is made dimensionless by using maximum value normalization. For each unit, take... i Maximum directional sensitivity And define the gradient direction normalized sensitivity as :

[0052] (twenty three);

[0053] in, J = The normalized sensitivity is limited to the interval [-1, 1].

[0054] S7, Design Variable Update:

[0055] Wrapping angle i and azimuth β Expressed in spherical coordinates:

[0056] (twenty four);

[0057] in, and They represent i and β Gradient direction normalized sensitivity, or Indicates the step size coefficient. x Indicates the basic angular update magnitude;

[0058] The gradient descent-based iterative process updates directly within a spherical coordinate system. The current iteration ( t The angle vector value is based on the previous iteration ( t The gradient information of -1) is calculated;

[0059] S8, Iteration Loop and Convergence Determination:

[0060] The fiber orientation field updated in step S7 ( i ,β The fiber orientation field is smoothed by performing a smoothing process. i , β Re-enter the finite element analysis module and repeat steps S4–S7 until the convergence condition is met. Finally, output the optimized fiber orientation distribution to achieve the design of maximizing structural stiffness.

[0061] In step S2, for the composite material winding layer including the cylinder section and the end cap section, the winding angle of the cylinder section is... i a constant, i a = i The end cap section is wrapped at the corner i e Following the axisymmetric variable angle distribution, and based on the process geometry constraints of the winding process, the winding angle of the cylinder section is established. i a Wrapping angle with end cap section i e Mapping relationship:

[0062] (2);

[0063] in, r Indicates the reference radius of the cylinder section. r e This indicates the radial distance between the head section unit and the axis of the composite pressure vessel.

[0064] In step S4,

[0065] (8);

[0066] (9);

[0067] in, sinth = S θ ,cosθ = C θ , sinβ = S β ,cosβ = C β .

[0068] In step S4, T ( i ')and T ( β ') By ​​the formula of the rotation axis T The formula for the axis of rotation is derived. T Combining the loop z Axis rotation Rz ( i ) and around y Axis rotation Ry ( β Transformation of )

[0069] (10);

[0070] (11);

[0071] (12);

[0072] Among them, coefficient a ij ( i = 1, 2, 3; j = 1, 2, 3) correspond to rotation matrices respectively. Rz ( i ) and Ry ( β Substituting the elements of (11) into (10) yields (8), and substituting (12) into (10) yields (9).

[0073] In step S5, to achieve the numerical solution of equation (15), the finite element method is used to analyze the composite material winding layer region. Oh H Discretize the liner region and the head region independently. Oh L Then the structural stiffness matrix K ( i , β ) is represented as:

[0074] (16);

[0075] Among them, the stiffness matrix of the lining and head unit and the stiffness matrix of composite winding layer unit They are defined as follows:

[0076] (17);

[0077] B Represents the element strain-displacement matrix. D e For the rotated tensor T ( i ') and T ( β The transformed elasticity matrix;

[0078] Within the finite element framework, the strain energy density is expressed as:

[0079] (18);

[0080] in, = Bu For the true strain vector, the direction component is obtained through the rotation tensor. T ( i ') and T ( β ') is directly related to the fiber orientation angle;

[0081] Substituting equation (18) into equation (15), we obtain equation (19).

[0082] In steps S1 and S5, the composite material winding layer region Oh H Discretization is performed using four-node tetrahedral elements.

[0083] In step S6, in order to derive the transformed elasticity matrix... D e For the winding angle i The partial derivatives of are used to differentiate equation (7):

[0084] (twenty one);

[0085] in,

[0086] (twenty two);

[0087] in, sinth = S θ ,cosθ = C θ , sinβ = S β ,cosβ = C β .

[0088] In step S7, the step size coefficient or A dynamic adjustment mechanism is adopted: in the initial stage, i.e. t When the value is ≤50, a larger value is chosen to accelerate gradient exploration, and then gradually decreased to stabilize the convergence process.

[0089] Compared with the prior art, the present invention has the following beneficial effects.

[0090] 1. High manufacturing feasibility: Fixed azimuth constraints are introduced in the finite element modeling stage, and geometric and winding process constraints are directly embedded into the calculation process, ensuring that the optimization results can be directly used for CNC winding machine processing without secondary adjustments.

[0091] 2. High computational efficiency: The analytical sensitivity of the winding angle is derived by direct differentiation, avoiding dependence on finite difference step size and a large number of repeated calculations, and maintaining high efficiency even when there are many design variables.

[0092] 3. Good optimization stability: The sensitivity results are normalized to eliminate the influence of differences in the magnitude of performance indicators in different regions on the optimization, thereby improving the stability and convergence speed of the optimization algorithm.

[0093] 4. Integrated Design: Geometric modeling, stiffness matrix calculation, finite element analysis, sensitivity calculation and optimization solution are integrated into a unified framework, achieving seamless connection from design to manufacturing. Attached Figure Description

[0094] To more clearly illustrate the specific embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the specific embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained from these drawings without creative effort.

[0095] Figure 1 A flowchart illustrating the sensitivity analysis and optimization method for fiber winding angle in composite pressure vessels;

[0096] Figure 2 This is a two-dimensional 1 / 4 liner model;

[0097] Figure 3 The 3D model after meshing;

[0098] Figure 4 The fiber orientation of the winding layer before optimization;

[0099] Figure 5 Optimized fiber orientation for the cylinder section;

[0100] Figure 6 The optimized fiber orientation for the end cap section. Detailed Implementation

[0101] The technical solution of the present invention will now be clearly and completely described with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0102] Example 1

[0103] like Figure 1As shown in the figure, this embodiment provides a method for analyzing and optimizing the fiber winding angle sensitivity of composite pressure vessels, including the following steps.

[0104] S1, Model Partitioning:

[0105] Composite material pressure vessel domain Oh Divided into mutually exclusive subfields:

[0106] (1);

[0107] in, Oh L This indicates the lining area and the end cap area. Oh H Indicates the region of the composite material winding layer;

[0108] Using the finite element method, the composite material winding layer region was analyzed. Oh H Discrete N e A system composed of 100 elements (using four-node tetrahedral elements), targeting Oh H Each unit within e , e =1,2, ..., N e It imparts homogenized orthotropic anisotropic material properties to arbitrary elements. e The fiber orientation is determined by a pair of angles ( β , i (Together described) β It is the azimuth angle. i The wrapping angle.

[0109] Composite material pressure vessels also include liners and heads. Oh L However, the fibers and their orientation, as design variables, only exist in the composite winding layer. Oh H middle.

[0110] The spatial variation of the fiber orientation angle defines the azimuth angle. β , winding angle i A two-parameter design scheme. Among them, β Indicates the fiber in the global coordinate system z - x Projection on a plane and z The angle between the axes, i Indicates the fiber axis relative to z - x The angle of inclination of a plane.

[0111] S2, Manufacturing Constraints:

[0112] Constraints are created by embedding a fixed azimuth angle principle into the finite element model. Under the constraint of the fixed azimuth angle principle, the azimuth angle of each element is... β The winding angle is considered a fixed parameter predetermined by the geometry of the composite pressure vessel, rather than a design variable. i As design variables, the optimization process only occurs when the variables are fixed. β Design variables under distribution i Adjustments will be made.

[0113] In step S2, for the composite material winding layer including the cylinder section and the end cap section, the winding angle of the cylinder section is... i a constant, i a = i The end cap section is wrapped at the corner i e Following the axisymmetric variable angle distribution, and based on the process geometry constraints of the winding process, the winding angle of the cylinder section is established. i a Wrapping angle with end cap section i e Mapping relationship:

[0114] (2);

[0115] in, r Indicates the reference radius of the cylinder section. r e This indicates the radial distance between the head section unit and the axis of the composite pressure vessel.

[0116] This mapping relationship originates from the geodesic principle of fiber paths, ensuring no fiber slippage during winding. Through this functional expression, design variables... i The winding angle is transferred from the cylinder section to each unit of the head section, achieving unified control of the entire layer winding angle while meeting manufacturing process constraints.

[0117] S3, Establish constitutive relations:

[0118] The liner and end cap are made of isotropic materials, and the stress vector is... With strain vector The constitutive equation between them is expressed as:

[0119] (3);

[0120] (4);

[0121] in, D L Here is the stiffness matrix. E For Young's modulus, Poisson's ratio;

[0122] To accurately characterize the anisotropic behavior of composite materials, the composite winding layer is equivalent to a macroscopically orthogonal anisotropic material at the microscale, and its mechanical behavior is completely determined by three mutually orthogonal principal material directions. A local coordinate system is established using the three mutually orthogonal principal material directions as coordinate axes. In the local coordinate system, the fiber axis is represented as direction 1, and the fiber radial direction is represented as directions 2 and 3. Assuming that the fiber axis is consistent with the global coordinate system, the stress vector... With strain vector The constitutive relation between them is expressed as:

[0123] (5);

[0124] (6);

[0125] in, It is the elastic matrix, serving as the original constitutive matrix and the basis for subsequent coordinate transformations to any direction;

[0126] ; ; ;

[0127] ; ;

[0128] ;

[0129] ; ; ;

[0130] The Young's modulus in the 1st direction. The Young's modulus in two directions. It represents the Young's modulus in three directions.

[0131] For the Poisson's ratio in the 1-2 plane, For the Poisson's ratio in the 1-3 plane, For the Poisson's ratio in the 2-3 plane,

[0132] The shear modulus of the 1-2 plane. The shear modulus of the 1-3 plane. For the shear modulus of the 2-3 plane, ;

[0133] Constrained by constitutive compatibility relations, according to ( i , j= 1,2,3; i ≠ j ),available , , .

[0134] S4, Parameterization of Design Variables:

[0135] For the original constitutive matrix Applying a rotational transformation yields a result containing fiber orientation ( β , i elasticity matrix D e :

[0136] (7);

[0137] in, T ( i ') and T ( β ')for i and β The rotation tensor;

[0138] (8);

[0139] (9);

[0140] in, sinth = S θ ,cosθ = C θ , sinβ = S β ,cosβ = C β ;

[0141] T ( i ') and T ( β ') By ​​the formula of the rotation axis T The formula for the axis of rotation is derived. T Combining the loop z Axis rotation Rz ( i ) and around y Axis rotation Ry ( β Transformation of )

[0142] (10);

[0143] (11);

[0144] (12);

[0145] Among them, coefficient a ij ( i = 1, 2, 3; j = 1, 2, 3) correspond to rotation matrices respectively. Rz ( i ) and Ry ( β Substituting the elements of (11) into (10) yields (8), and substituting (12) into (10) yields (9).

[0146] S5, Mathematical formulation of the optimization problem:

[0147] fiber orientation field ( i , β By incorporating an orientation-dependent material stiffness matrix into the finite element model, multi-scale coupling between the microscopic orientation of the material and the macroscopic structural response is achieved, thereby optimizing the winding angle of composite materials with the goal of maximizing structural stiffness.

[0148] The total strain energy of a composite pressure vessel under a given load. Represented as:

[0149] (13);

[0150] in, J L This represents the strain energy of the lining and end caps, independent of fiber orientation. J H ( i ) indicates the angle of wrapping. i The strain energies of the relevant composite winding layers are as follows:

[0151] (14);

[0152] in, This indicates the strain energy density of the lining and end caps. This represents the strain energy density of the composite winding layer;

[0153] The fiber orientation optimization problem can be expressed in the following mathematical form:

[0154] (15);

[0155] in, Describe the objective function. K ( i , β This indicates dependence on the fiber orientation field. i , β The structural stiffness matrix, u Represents the displacement vector. F Represents the load vector;

[0156] To achieve a numerical solution for equation (15), the finite element method is used to analyze the composite material winding layer region. Oh H Discretize the material (using four-node tetrahedral elements), and simultaneously discretize the inner liner region and the head region independently. Oh L Then the structural stiffness matrix K ( i , β ) is represented as:

[0157] (16);

[0158] Among them, the stiffness matrix of the lining and head unit and the stiffness matrix of composite winding layer unit They are defined as follows:

[0159] (17);

[0160] B Represents the element strain-displacement matrix. D e For the rotated tensor T ( i ') and T ( β The transformed elastic matrix allows for the modulation of the material's anisotropic stiffness by fiber orientation.

[0161] Within the finite element framework, the strain energy density is expressed as:

[0162] (18);

[0163] in, = Bu For the true strain vector, the direction component is obtained through the rotation tensor. T ( i ') and T ( β ') is directly related to the fiber orientation angle;

[0164] Substituting equation (18) into equation (15), the integral form of the optimized objective function is:

[0165] (19);

[0166] Equation (19) reveals the fiber orientation field of the composite winding layer. i , β The direct coupling relationship between the structure and its stiffness can be adjusted by... i The distribution of this distribution allows for quantitative optimization of structural stiffness while maintaining manufacturing constraints and geometric consistency.

[0167] S6, Sensitivity Analysis Solution:

[0168] Sensitivity analysis fundamentally establishes the objective function. With fiber orientation field ( i , β The analytical gradient relationship between them, due to the azimuth angle of each element after geometric discretization. β This has been determined, therefore the sensitivity analysis only considers the winding angle. i conduct.

[0169] The sensitivity derivation is performed using the adjoint variable method, which introduces a virtual strain tensor under the assumption of material isotropy. Alternative For the winding angle i Sensitivity analysis is performed to avoid the tedious chain-reaction process:

[0170] (20);

[0171] To derive the transformed elasticity matrix D e For the winding angle i The partial derivatives of are used to differentiate equation (7):

[0172] (twenty one);

[0173] in,

[0174] (twenty two);

[0175] in, sinth = S θ ,cosθ = C θ , sinβ = S β ,cosβ = C β ;

[0176] In structural optimization, magnitude differences in sensitivity values ​​can lead to difficulties in iterative convergence or ill-conditioned optimization. Therefore, maximum value normalization is used to make the sensitivity dimensionless. For each element, a value is taken as... i Maximum directional sensitivity And define the gradient direction normalized sensitivity as :

[0177] (twenty three);

[0178] in, J = The normalized sensitivity is limited to the interval [-1, 1].

[0179] S7, Design Variable Update:

[0180] Wrapping angle i and azimuth β Expressed in spherical coordinates:

[0181] (twenty four);

[0182] in, and They represent i and β Gradient direction normalized sensitivity, or Indicates the step size coefficient. x This indicates the basic angular update magnitude.

[0183] Step size coefficient or A dynamic adjustment mechanism is adopted: in the initial stage, i.e. t When the value is ≤50, a larger value is chosen to accelerate gradient exploration, and then gradually decreased to stabilize the convergence process.

[0184] Basic angular update amplitude x This ensured a physically reasonable change in orientation.

[0185] The gradient descent-based iterative process updates directly within a spherical coordinate system. The current iteration ( t The angle vector value is based on the previous iteration ( t The gradient information of -1) is calculated, and this method preserves the geometric interpretability of angle updates.

[0186] S8, Iteration Loop and Convergence Determination:

[0187] The fiber orientation field updated in step S7 ( i , β The fiber orientation field is smoothed by performing a smoothing process. i , βRe-enter the finite element analysis module and repeat steps S4–S7 until the convergence condition is met. Finally, output the optimized fiber orientation distribution to achieve the design of maximizing structural stiffness.

[0188] Example 2

[0189] This embodiment uses a three-dimensional example of a single-layer composite material pressure vessel for detailed explanation.

[0190] A composite pressure vessel with an inner diameter of 132.6 mm (h), a semi-length of 182.5 mm (j), and a hemispherical head is used as the research object. The vessel consists of an aluminum alloy liner and a carbon fiber / epoxy resin winding layer. The liner thickness is 6.1 mm (i), and the winding layer thickness is designed to be 2 mm. A two-dimensional 1 / 4 model of the liner is shown below. Figure 2 As shown. The 3D model after meshing is as follows. Figure 3 As shown. The orientation of the composite winding layer before optimization is as follows. Figure 4 As shown in the figure, different colors represent transitional colors. Fibers biased towards the x-axis are displayed in red, fibers biased towards the y-axis are displayed in green, and fibers biased towards the z-axis are displayed in blue. The initial fiber direction (whether it is the cylinder section or the end cap section) is parallel to the x-axis, so it is displayed in red.

[0191] The aluminum alloy liner uses an isotropic material model with an elastic modulus of 74.12 GPa and a Poisson's ratio of 0.28. The carbon fiber composite material uses an orthotropic material model with a fiber axial Young's modulus of... The radial elastic modulus is 141 GPa. = The shear modulus is 9.14 GPa. = The Pa is 3.57 GPa. The Pa is 4.79 GPa, and the Poisson's ratio is... =0.28, =0.3, =0.3.

[0192] In the process of calculating the element stiffness matrix, the azimuth angle is... β and wrapping angle i The material stiffness matrix is ​​input separately and then oriented through rotation transformation to obtain the element stiffness matrix containing the influence of geometric constraints and design variables, which is then assembled into the structural stiffness matrix.

[0193] After applying internal pressure loads and setting constraints, the displacement and stress fields of the vessel are solved using the finite element method, and performance indicators such as maximum equivalent stress, maximum radial displacement, and buckling load factor are calculated. These indicators are used as optimization objectives and constraints, and their analytical sensitivity formulas to the winding angle are derived using the direct differential method. The sensitivity value of each element is calculated and normalized.

[0194] During the optimization phase, the wrapping angle will be... i Design variables are constrained within the allowable range of the process, with the goal of minimizing strain energy, and an adaptive gradient descent algorithm is used for iterative optimization. In each iteration, the winding angle distribution is updated based on sensitivity information, and finite element calculations and sensitivity analyses are performed again until the objective function converges to a set threshold or the maximum number of iterations is reached.

[0195] After optimization, the obtained winding angle distribution is mapped to the actual fiber path on the container surface, and the feasibility of the winding trajectory is verified using the geodesic principle, ensuring that the optimization results can be directly imported into the control program of the CNC winding machine for processing. The optimized fiber direction of the cylinder section is as follows. Figure 5 As shown, the optimized fiber orientation of the end cap section is as follows: Figure 6 As shown, with Figure 4 Similarly, the different colors in the graph represent transitional colors: fibers biased towards the x-axis are displayed as red, fibers biased towards the y-axis as green, and fibers biased towards the z-axis as blue. Figure 5 and Figure 6 It can produce a non-crossing winding angle distribution, which can be used to guide actual winding.

[0196] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for analyzing and optimizing the sensitivity of fiber winding angle in composite pressure vessels, characterized in that, Includes the following steps: S1, Model Partitioning: Composite material pressure vessel domain Ω Divided into mutually exclusive subfields: （1）； in, Ω L This indicates the lining area and the end cap area. Ω H Indicates the region of the composite material winding layer; Using the finite element method, the composite material winding layer region was analyzed. Ω H Discrete N e A system composed of individual units, targeting Ω H Each unit within e , e =1,2, ..., N e It imparts homogenized orthotropic anisotropic material properties to arbitrary elements. e The fiber orientation is determined by a pair of angles ( β , θ (Together described) β It is the azimuth angle. θ For the wrapping angle; S2, Manufacturing Constraints: Constraints are created by embedding a fixed azimuth angle principle into the finite element model. Under the constraint of the fixed azimuth angle principle, the azimuth angle of each element is... β The winding angle is considered a fixed parameter predetermined by the geometry of the composite pressure vessel, rather than a design variable. θ As design variables, the optimization process only occurs when the variables are fixed. β Design variables under distribution θ Make adjustments; S3, Establish constitutive relations: The liner and end cap are made of isotropic materials, and the stress vector is... With strain vector The constitutive equation between them is expressed as: （3）； （4）； in, D L Here is the stiffness matrix. E For Young's modulus, Poisson's ratio; At the microscale, the composite winding layer is equivalent to a macroscopic orthogonal anisotropic material, and a local coordinate system is established with three mutually orthogonal principal directions of the material as coordinate axes. In the local coordinate system, the fiber axis is represented as direction 1, and the fiber radial direction is represented as directions 2 and 3. Assuming that the fiber axis is consistent with the global coordinate system, the stress vector... With strain vector The constitutive relation between them is expressed as: （5）； （6）； in, It is the elastic matrix, serving as the original constitutive matrix and the basis for subsequent coordinate transformations to any direction; ； ； ； ； ； ； ； ； ； The Young's modulus in the 1st direction. The Young's modulus in two directions. It represents the Young's modulus in three directions. For the Poisson's ratio in the 1-2 plane, For the Poisson's ratio in the 1-3 plane, For the Poisson's ratio in the 2-3 plane, The shear modulus of the 1-2 plane. The shear modulus of the 1-3 plane. For the shear modulus of the 2-3 plane, ; Constrained by constitutive compatibility relations, according to ( i , j = 1,2,3; i ≠ j ),available , , ; S4, Parameterization of Design Variables: For the original constitutive matrix Applying a rotational transformation yields a result containing fiber orientation ( β , θ elasticity matrix D e : （7）； in, T ( θ ') and T ( β ')for θ and β The rotation tensor; S5, Mathematical formulation of the optimization problem: The total strain energy of a composite pressure vessel under a given load. Represented as: （13）； in, J L This represents the strain energy of the lining and end caps, independent of fiber orientation. J H ( θ ) indicates the angle of wrapping. θ The strain energies of the relevant composite winding layers are as follows: （14）； in, This indicates the strain energy density of the lining and end caps. This represents the strain energy density of the composite winding layer; The fiber orientation optimization problem can be expressed in the following mathematical form: （15）； in, Describe the objective function. K ( θ , β This indicates dependence on the fiber orientation field. θ , β The structural stiffness matrix, u Represents the displacement vector. F Represents the load vector; The integral form of the optimized objective function: （19）； in, = Bu For the true strain vector, the direction component is obtained through the rotation tensor. T ( θ ')and T ( β ') is directly related to the fiber orientation angle; S6, Sensitivity Analysis Solution: Due to the azimuth angle of each element after geometric discretization β This has been determined, therefore the sensitivity analysis only considers the winding angle. θ The sensitivity was derived using the adjoint variable method, and a virtual strain tensor was introduced under the assumption of material isotropy. Alternative For the winding angle θ Perform sensitivity analysis: （20）； The sensitivity is made dimensionless by using maximum value normalization. For each unit, take... θ Maximum directional sensitivity And define the gradient direction normalized sensitivity as : （23）； in, J = The normalized sensitivity is limited to the interval [-1, 1]. S7, Design Variable Update: Wrapping angle θ and azimuth β Expressed in spherical coordinates: （24）； in, and They represent θ and β Gradient direction normalized sensitivity, η Indicates the step size coefficient. ξ Indicates the basic angular update magnitude; The gradient descent-based iterative process updates directly within a spherical coordinate system. The current iteration ( t The angle vector value is based on the previous iteration ( t The gradient information of -1) is calculated; S8, Iteration Loop and Convergence Determination: The fiber orientation field updated in step S7 ( θ , β The fiber orientation field is smoothed by performing a smoothing process. θ , β Re-enter the finite element analysis module and repeat steps S4–S7 until the convergence condition is met. Finally, output the optimized fiber orientation distribution to achieve the design of maximizing structural stiffness.

2. The method for analyzing and optimizing the sensitivity of fiber winding angle in composite pressure vessels according to claim 1, characterized in that, In step S2, for the composite material winding layer including the cylinder section and the end cap section, the winding angle of the cylinder section is... θ a constant, θ a = θ The end cap section is wrapped at the corner θ e Following the axisymmetric variable angle distribution, and based on the process geometry constraints of the winding process, the winding angle of the cylinder section is established. θ a Wrapping angle with end cap section θ e Mapping relationship: （2）； in, r Indicates the reference radius of the cylinder section. r e This indicates the radial distance between the head section unit and the axis of the composite pressure vessel.

3. The method for analyzing and optimizing the sensitivity of fiber winding angle in composite pressure vessels according to claim 1, characterized in that, In step S4, （8）； （9）； in, sinθ = S θ cosθ = C θ , sinβ = S β cosβ = C β .

4. The method for analyzing and optimizing the sensitivity of fiber winding angle in composite pressure vessels according to claim 3, characterized in that, In step S4, T ( θ ') and T ( β ') By ​​the formula of the rotation axis T The formula for the axis of rotation is derived. T Combining the loop z Axis rotation Rz ( θ ) and around y Axis rotation Ry ( β Transformation of ) （10）； （11）； （12）； Among them, coefficient a ij ( i = 1, 2, 3; j = 1, 2, 3) correspond to rotation matrices respectively. Rz ( θ ) and Ry ( β Substituting the elements of (11) into (10) yields (8), and substituting (12) into (10) yields (9).

5. The method for analyzing and optimizing the sensitivity of fiber winding angle in composite pressure vessels according to claim 1, characterized in that, In step S5, to achieve the numerical solution of equation (15), the finite element method is used to analyze the composite material winding layer region. Ω H Discretize the liner region and the head region independently. Ω L Then the structural stiffness matrix K ( θ , β ) is represented as: （16）； Among them, the stiffness matrix of the lining and head unit and the stiffness matrix of composite winding layer unit They are defined as follows: （17）； B Represents the element strain-displacement matrix. D e For the rotated tensor T ( θ ') and T ( β The transformed elastic matrix allows for the modulation of the material's anisotropic stiffness by fiber orientation. Within the finite element framework, the strain energy density is expressed as: （18）； in, = Bu For the true strain vector, the direction component is obtained through the rotation tensor. T ( θ ') and T ( β ') is directly related to the fiber orientation angle; Substituting equation (18) into equation (15), we obtain equation (19).

6. The method for analyzing and optimizing the sensitivity of fiber winding angle in composite pressure vessels according to claim 5, characterized in that, In steps S1 and S5, the composite material winding layer region Ω H Discretization is performed using four-node tetrahedral elements.

7. The method for sensitivity analysis and optimization of fiber winding angle in composite pressure vessels according to claim 1, characterized in that, In step S6, in order to derive the transformed elasticity matrix... D e For the winding angle θ The partial derivatives of are used to differentiate equation (7): （21）； in, （22）； in, sinθ = S θ cosθ = C θ , sinβ = S β cosβ = C β .

8. The method for analyzing and optimizing the sensitivity of fiber winding angle in composite pressure vessels according to claim 1, characterized in that, In step S7, the step size coefficient η A dynamic adjustment mechanism is adopted: in the initial stage, i.e. t When the value is ≤50, a larger value is chosen to accelerate gradient exploration, and then gradually decreased to stabilize the convergence process.

Images