Method for calculating element stiffness matrix of winding layer seal head section of composite material pressure vessel
By calculating the azimuth angle and winding angle in the head section of the winding layer of the composite pressure vessel, combined with coordinate transformation, the problem of large calculation amount and insufficient accuracy in the traditional method is solved, and efficient and accurate unit stiffness matrix calculation is achieved, supporting design and analysis.
Patent Information
- Application Number
- CN202510834002.9
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2025-06-20
- Publication Date
- 2025-07-22
- Estimated Expiration
- Not applicable · inactive patent
AI Technical Summary
Traditional methods use large amounts of calculation and insufficient accuracy when calculating the head section stiffness matrix of the composite pressure vessel winding layer, which affects design and safety analysis.
Gmsh software is used to extract node coordinates, establish a global coordinate system, and accurately calculate the unit stiffness matrix by calculating azimuth and winding angles, combining coordinate transformation and finite element methods.
Improves computational accuracy and efficiency, ensures correct expression of material properties, and supports structural design and performance analysis.
Smart Images
Figure CN120354557A_ABST
Abstract
Description
Technical Field
[0001] The present invention belongs to the technical field of composite material winding optimization, and specifically discloses a method for calculating the element stiffness matrix of the head section of a composite pressure vessel winding layer. Background Art
[0002] Composite pressure vessels are widely used in the fields of aerospace, petrochemical, deep-sea exploration, etc. due to their excellent mechanical properties, light weight, and corrosion resistance. As a key structural part of composite pressure vessels, the stiffness characteristics of the head section of the winding layer directly affect the load-bearing capacity and safety performance of the entire pressure vessel. The accurate solution of the stiffness matrix of the head section of the winding layer is crucial for the design and safety analysis of pressure vessels. However, due to the complex geometric shape of the head section of the winding layer and the mechanical properties of multi-layer composite materials, traditional methods have problems of large computational amount and insufficient accuracy when calculating its element stiffness matrix. Summary of the Invention
[0003] The present invention provides a method for calculating the element stiffness matrix of the head section of a composite pressure vessel winding layer, which improves the problems of large computational amount and insufficient accuracy existing in traditional methods when calculating the element stiffness matrix of the head section.
[0004] The above method for calculating the element stiffness matrix of the head section of a composite pressure vessel winding layer includes the following steps: S1, Element information acquisition and preprocessing: t1, Extract the coordinate information of each node in the three-dimensional space in the head section from the design model or simulation model, number each node, determine the geometric shape of each element in the head section, the numbers of the connected nodes in each element, and the winding layer-related parameters corresponding to each element; t2, Taking the geometric center point of the pressure vessel as the central coordinate according to the coordinate information of the nodes o , taking the central axis of the pressure vessel as y axis, taking the initially selected generatrix of the pressure vessel as z axis, taking the straight line perpendicular to y axis and z axis and passing through the central coordinate o as x axis, establish a global coordinate system, and calculate the centroid coordinates of all elements in the head section according to the coordinate information of the nodes; S2, Azimuth angle and winding angle calculation: Calculate the distance y in the radial direction between the parallel circle where the centroid of each element is located and the r axis, and use r and the reference value of the winding angle of the cylinder section to calculate the winding angle α' of each element in the head section; Through the relative positions of the body center coordinates and the center coordinates of each unit o , calculate the rotation angle of each unit around the y axis φ' , and determine the azimuth angle of each unit in the head section in the global coordinate system according to the quadrant of φ' in the global coordinate system β' ; S3, Coordinate transformation and elastic matrix conversion: According to the azimuth angle β' and the winding angle α' calculated in step S2, construct the rotation matrices of the unit around the z axis and the y axis. The rotation matrices include the forward coordinate transformation matrix and the reverse coordinate transformation matrix. Using the coordinate transformation matrix, convert the elastic matrix of the winding layer material in the local coordinate system to the global coordinate system to obtain the global elastic matrix; S4, Element stiffness matrix calculation: According to the global elastic matrix and the geometric shape parameters of the unit, calculate the element stiffness matrix using the standard formula of the finite element method.
[0005] Step S1 is carried out using Gmsh software.
[0006] In step S1, the geometric shape of the unit is a tetrahedron.
[0007] In step S1, the parameters related to the winding layer include the type, number of layers and thickness of the fiber.
[0008] In step S2, take the center coordinate o as ( x, y, z ), where x, y, z are all taken as 0. The body center coordinates of all units in the head section are ( x 1, y 1, z 1), then ; The winding angle α' = sin -1 of each unit in the head section is r e sin α e / r ); r e sin α e is the reference value of the winding angle of the cylinder section, r e is the radius of the cylinder section, α e is the winding angle of the cylinder section.
[0009] In step S2, each unit of the head section rotates around the y axis by an angle of tan φ' = ( x 1 - x ) / ( z 1 - z ). The position of each unit of the head section in the global coordinate system is expressed as ( y 1, φ' ); If z 1 > 0, then the azimuth angle β' = φ' of each unit of the head section; If z 1 < 0, x 1 > 0, then the azimuth angle β' = 180° + φ' ; If z 1 < 0, x 1 < 0, then the azimuth angle β' = - 180° + φ' .
[0010] In step S3, first rotate around the z axis to obtain the forward coordinate transformation matrix , then rotate around the y axis to obtain the forward coordinate transformation matrix . The elastic matrix of the winding layer material in the global coordinate system after forward coordinate transformation is: ; where D is the elastic matrix of the winding layer material in the local coordinate system, and are respectively the transposes of and ; Calculate the reverse coordinate transformation matrices and . The elastic matrix of the winding layer material in the global coordinate system after reverse coordinate transformation is: ; where D is the elastic matrix of the winding layer material in the local coordinate system, and are respectively the transposes of and ; The global elastic matrix is: .
[0011] In step S3, ; ; Among them sin α' = S α , cos α' = C α ; sin β' = S β , cos β' = C β ; ; ; Among them sin ( -α' ) = S -α , cos ( -α' ) = C -α ; sin ( -β' ) = S -β , cos ( -β' ) = C -β .
[0012] In step S4, for each tetrahedral element, the geometric shape parameters of the element include the volume and the gradient matrix of the shape function B ; According to the node numbers and node coordinate information of the four vertices in each element obtained in step S1, the volume of each element is calculated as: ; Among them, r 12 is the vector from vertex 1 to vertex 2, r 13 is the vector from vertex 1 to vertex 3, r 14 is the vector from vertex 1 to vertex 4; The element stiffness matrix is: ; Among them, B is the gradient matrix of the shape function, V is the volume of the element.
[0013] Compared with the prior art, the present invention has the following beneficial effects.
[0014] High calculation accuracy: By accurately calculating the fiber orientation (i.e., azimuth angle and winding angle) of the winding layer and combining coordinate transformation technology, the elastic matrix in the local coordinate system is accurately transformed into the global coordinate system in the present invention, and then the stiffness matrix of each element is accurately calculated, ensuring the correct expression of material properties under different loading conditions, thereby improving the calculation accuracy.
[0015] Efficient method: Compared with traditional methods, the present invention provides a systematic solution process, including steps such as node coordinate extraction, element information determination, azimuth angle and winding angle calculation, coordinate transformation and elastic matrix transformation, and element stiffness matrix calculation. These steps are interconnected to form an efficient solution system, reducing redundancy and repetition in the calculation process, improving the solution efficiency, and facilitating subsequent finite element calculations.
[0016] Strong applicability: The present invention is not only applicable to the winding layer head section of composite pressure vessels, but can also be extended to other fields with similar structures or materials. Its method has generality and portability, providing strong support for research and applications in related fields.
[0017] Support for design and analysis: The solution method provided by the present invention can accurately calculate the element stiffness matrix of composite pressure vessels, which is crucial for structural design and performance analysis. Designers can perform more accurate structural optimization and performance prediction based on the results of this method, improving the overall performance of the product. Description of the drawings
[0018] In order to more clearly illustrate the specific embodiments of the present invention or the technical solutions in the prior art, the following will briefly introduce the drawings required for use in the description of the specific embodiments or the prior art. Obviously, the following drawings are some embodiments of the present invention. For those of ordinary skill in the art, without creative efforts, other drawings can be obtained based on these drawings.
[0019] Figure 1 It is a flowchart of the method for calculating the element stiffness matrix of the winding layer head section of a composite pressure vessel; Figure 2 It is a definition diagram of fiber orientation; Figure 3 It is a position diagram of the head section element in the global coordinate system; Figure 4 For Figure 3 Top view; Figure 5 For the rotation angle φ' Schematic diagrams when located in different quadrants; Figure 6 For the rotation angle φ' And the azimuth angle β'Schematic diagrams in different quadrants. Detailed implementation manners
[0020] The technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are some, but not all, of the embodiments of the present invention. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.
[0021] This embodiment provides a method for calculating the stiffness matrix of the winding layer head segment unit of a composite pressure vessel. This method accurately calculates the fiber orientations (i.e., azimuth angle and winding angle) of the winding layer, and combines coordinate transformation to convert the elastic matrix in the local coordinate system to the global coordinate system, thereby accurately calculating the stiffness matrix of each unit. The specific steps are as follows.
[0022] S1, Obtaining and preprocessing unit information.
[0023] t1, Use Gmsh software to extract the coordinate information of each node in the three-dimensional space in the head segment from the design model or simulation model, number each node, determine the geometric shape of each unit in the head segment, the numbers of the connected nodes in each unit, and the relevant parameters of the winding layer corresponding to each unit.
[0024] t2, Taking the geometric center point of the pressure vessel as the central coordinate o and taking the central axis of the pressure vessel as the y axis, taking the initially selected bus bar of the pressure vessel as the z axis, taking the straight line perpendicular to the y axis and the z axis and passing through the central coordinate o as the x axis, establish a global coordinate system, and calculate the centroid coordinates of all units in the head segment according to the coordinate information of the nodes.
[0025] In step S1, the geometric shape of the unit is a tetrahedron.
[0026] In step S1, the relevant parameters of the winding layer include the type, number of layers, and thickness of the fiber.
[0027] S2, Calculating the azimuth angle and winding angle.
[0028] As Figure 2 shown, in the three-dimensional space, use a method based on the azimuth angle β and the winding angle (also called elevation angle) αThe coordinate system defines the fiber orientation; the winding angle defines the directionality of the fiber on the unit surface and directly affects the stiffness characteristics of the unit; the azimuth angle is used to describe the direction of the unit in the winding layer and is usually expressed in radians.
[0029] Calculate the distance in the radial direction between the parallel circle passing through the center of the body of each unit (the parallel circle obtained by rotating the center of the body of the unit around the y axis) and the y axis r , and use r and the reference value of the winding angle of the cylindrical section to calculate the winding angle α' of each unit in the head section.
[0030] As Figure 3 and Figure 4 shown, specifically, take the central coordinate o as ( x, y, z ), where x, y, z are all taken as 0, and the central coordinates of all units in the head section are ( x 1, y 1, z 1), and the distance between the central coordinate of each unit in the head section and the z axis is ( x 2 - x ), then ; The winding angle α' = sin -1 of each unit in the head section is ( r e sin α e / r ); r e sin α e is the reference value of the winding angle of the cylindrical section, r e is the radius of the cylindrical section, α e is the winding angle of the cylindrical section.
[0031] By the relative position of the central coordinate of each unit and the central coordinate o , calculate the angle y rotated by each unit around the φ' axis (i.e., the rotation angle), and determine the azimuth angle φ' of each unit in the head section in the global coordinate system according to the quadrant of β' in the global coordinate system.
[0032] Specifically, the angle rotated by each unit in the head section around the y axis is tan φ' = ( x1 - x ) / ( z 1 - z ), the position of each unit in the head section is represented as ( y 1, φ' ); As Figure 5 and Figure 6 shown, φ' if it is in the fourth quadrant, the value is positive, β' = φ' , in the figure, φ' 4 and β' 4 represent φ' and β' the relationship when in the fourth quadrant; φ' if it is in the first quadrant, the value is negative, β' = 180° + φ' , in the figure, φ' 1 and β' 1 represent φ' and β' the relationship when in the first quadrant; φ' if it is in the second quadrant, the value is positive, β' = -180° + φ' , in the figure, φ' 2 and β' 2 represent φ' and β' the relationship when in the second quadrant; φ' if it is in the third quadrant, the value is negative, β' = φ' , in the figure, φ' 3 and β' 3 represent φ' and β' the relationship when in the third quadrant; Combining the third and fourth quadrants, and keeping the first and second quadrants unchanged, we get: If z 1 > 0, then the azimuth angle β' = φ' ; If z 1 < 0, x 1 > 0, then the azimuth angle β' = 180° + φ' ; If z 1 < 0, x 1 < 0, then the azimuth angle β' = - 180° + φ' .
[0033] S3, Coordinate transformation and elastic matrix transformation.
[0034] The azimuth angle calculated according to step S2 β' and the winding angle α' , a rotation matrix of the building unit about the z axis and the y axis is constructed. The rotation matrix includes a forward coordinate transformation matrix and a reverse coordinate transformation matrix. Using the coordinate transformation matrix, the elastic matrix of the winding layer material in the local coordinate system is transformed into the global coordinate system to obtain the global elastic matrix.
[0035] Specifically, first rotate about the z axis to obtain the forward coordinate transformation matrix , ; Then rotate about the y axis to obtain the forward coordinate transformation matrix , ; where sin α' = S α , cos α' = C α ; sin β' = S β , cos β' = C β ; The elastic matrix of the winding layer material in the global coordinate system after forward coordinate transformation is: ; where D is the elastic matrix of the winding layer material in the local coordinate system, and are respectively the transposes of and ; Since one - turn spiral winding includes two times of forward and reverse winding, calculate the reverse coordinate transformation matrices and , ; ; where sin ( -α' )= S -α , cos ( -α' )= C -α ; sin ( -β' )= S-β , cos ( -β' )= C -β ; The elastic matrix of the winding layer material in the global coordinate system after reverse coordinate transformation is: ; where D is the elastic matrix of the winding layer material in the local coordinate system, and are respectively and transpose; Since one - turn helical winding contains two times of positive and negative winding angles, after the three - dimensional stiffness matrices of positive and negative angles are transformed respectively and then averaged, the global elastic matrix is: .
[0036] Through the positive and negative coordinate transformation matrices, the local elastic matrix of the winding layer head section is transformed into the global elastic matrix. This step ensures the accuracy of the stiffness calculation of each element in the global coordinate system.
[0037] S4, Element stiffness matrix calculation.
[0038] According to the global elastic matrix and the geometric shape parameters of the element, the element stiffness matrix is calculated using the standard formula of the finite - element method. The stiffness matrix describes the deformation behavior of the element under different loading conditions.
[0039] Specifically, for each tetrahedral element, the geometric shape parameters of the element include the volume and the gradient matrix of the shape function B , which involves the coordinates of the element vertices and vector operations.
[0040] According to the node numbers and node coordinate information of the four vertices in each element obtained in step S1, the volume of each element is obtained as: ; where, r 12 is the vector from vertex 1 to vertex 2, r 13 is the vector from vertex 1 to vertex 3, r 14 is the vector from vertex 1 to vertex 4; For each element, the element stiffness matrix K can be calculated by the standard formula of the finite - element method, which involves the global elastic matrix D e and the geometric shape of the element, and the element stiffness matrix is: ; Among them, B is the gradient matrix of the shape function, V is the volume of the element.
[0041] Finally, it should be noted that: The above embodiments are only used to illustrate the technical solutions of the present invention, rather than limiting them; Although the present invention has been described in detail with reference to the foregoing embodiments, those of ordinary skill in the art should understand that: They can still modify the technical solutions recorded in the foregoing embodiments, or perform equivalent replacements on some or all of the technical features; And these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the scope of the technical solutions of the embodiments of the present invention.
Claims
1. A calculation method for the stiffness matrix of the winding layer head section unit of a composite material pressure vessel, characterized in that It includes the following steps: S1, Unit information acquisition and preprocessing: t1, Extract the coordinate information of each node in the three-dimensional space in the head section from the design model or simulation model, number each node, determine the geometric shape of each unit in the head section, the numbers of the connected nodes in each unit, and the winding layer related parameters corresponding to each unit; t2, with the geometric center point of the pressure vessel as the central coordinate according to the coordinate information of the nodes o , with the central axis of the pressure vessel as the y axis, with the generatrix of the initially selected pressure vessel as the z axis, with the line perpendicular to the y axis and the z axis and passing through the central coordinate o as the x axis, establish a global coordinate system, and calculate the centroid coordinates of all elements in the head section according to the coordinate information of the nodes; S2, Azimuth angle and winding angle calculation: Calculate the distance in the radial direction between the parallel circle where the body center of each unit is located and the y axis r , and use the r and the reference value of the winding angle of the barrel section to calculate the winding angle of each unit in the head section α' ; Based on the relative positions of the body-centered coordinates and the central coordinates of each cell o calculate the angle y by which each cell rotates about the φ' axis, and determine the azimuth angle φ' in the global coordinate system of each cell in the head section according to the quadrant β' in the global coordinate system; S3, Coordinate transformation and elastic matrix conversion: The azimuth angle calculated according to step S2 β' and the winding angle α' , construct the rotation matrix of the unit around the z axis and the y axis. The rotation matrix includes a forward coordinate transformation matrix and a reverse coordinate transformation matrix. Using the coordinate transformation matrix, the elastic matrix of the winding layer material in the local coordinate system is transformed to the global coordinate system to obtain the global elastic matrix; S4, Element stiffness matrix calculation: According to the global elastic matrix and the geometric shape parameters of the element, use the standard formula of the finite element method to calculate the element stiffness matrix.
2. The calculation method of the stiffness matrix of the winding layer head section unit of the composite material pressure vessel according to claim 1, wherein, Step S1 is carried out using Gmsh software.
3. The calculation method of the stiffness matrix of the winding layer head section unit of the composite material pressure vessel according to claim 2, characterized in that, In step S1, the geometric shape of the unit is a tetrahedron.
4. The calculation method of the stiffness matrix of the winding layer head section unit of the composite material pressure vessel according to claim 3, wherein In step S1, the winding layer related parameters include the type, number of layers and thickness of the fiber.
5. The calculation method of the unit stiffness matrix of the winding layer head section of the composite material pressure vessel according to claim 4, characterized in that, In step S2, take the central coordinates o as ( x, y, z ), where x, y, z are all taken as 0, and the body center coordinates of all units in the head segment are ( x 1, y 1, z 1), then ; Winding angle of each unit in the head section α' =sin -1 ( r e sinα e / r ) r e sinα e is the reference value of the winding angle of the cylinder section, r e is the radius of the cylinder section, α e is the winding angle of the cylinder section.
6. The calculation method of the stiffness matrix of the winding layer head section unit of the composite material pressure vessel according to claim 5, characterized in that In step S2, the angle by which each unit in the head section rotates around y the axis is tanφ' = ( x 1 - x ) / ( z 1 - z ). The position of each unit in the head section is represented in the global coordinate system as ( y 1, φ' ); If z 1 > 0, the azimuth angle of each unit in the head section β'=φ' ; If z 1 < 0, x 1 > 0, then the azimuth angle of each unit in the head section β' = 180° + φ' ; If z 1 < 0, x 1 < 0, then the azimuth angle of each unit in the head section β' = - 180° + φ' .
7. The calculation method of the stiffness matrix of the winding layer head section unit of the composite material pressure vessel according to claim 6, characterized in that, In step S3, first rotate around z axis to obtain the forward coordinate transformation matrix , then rotate around y axis to obtain the forward coordinate transformation matrix . The elastic matrix of the winding layer material in the global coordinate system after forward coordinate transformation is: ； Among them D is the elastic matrix of the winding layer material in the local coordinate system, and are respectively and the transposes of; Calculate the inverse coordinate transformation matrix and , the elastic matrix of the winding layer material in the global coordinate system after inverse coordinate transformation is: ； wherein D is the elastic matrix of the winding layer material in the local coordinate system, and are respectively and the transposes of; The global elastic matrix is: 。 8. The calculation method of the stiffness matrix of the winding layer head section unit of the composite material pressure vessel according to claim 7, characterized in that, In step S3, ； ； Among them sinα' = S α , cosα' = C α ; sinβ' = S β , cosβ' = C β ; ； ； Among them sin ( -α' ) = S -α , cos ( -α' ) = C -α ; sin ( -β' ) = S -β , cos ( -β' ) = C -β .
9. The calculation method of the stiffness matrix of the winding layer head section unit of the composite material pressure vessel according to claim 8, characterized in that, In step S4, for each tetrahedral element, the geometric shape parameters of the element include the volume and the gradient matrix of the shape function B ; According to the node numbers and node coordinate information of the four vertices in each unit obtained in step S1, calculate the volume of each unit as: ； Among them, r 12 is the vector from vertex 1 to vertex 2, r 13 is the vector from vertex 1 to vertex 3, r 14 is the vector from vertex 1 to vertex 4; The element stiffness matrix is: ； Among them, B is the gradient matrix of the shape function, V is the volume of the element.