Finite element-based modeling method for variable stiffness composite laminated plate

Abstract

The invention belongs to the field of composite material structure optimization design, and particularly relates to a modeling method of a variable stiffness composite material laminated plate based on a finite element. The method comprises the following steps: defining parameters of a composite laminated plate model; compiling a script file by utilizing a programming language, and dividing a quadrilateral grid; calculating geometric coordinates of the unit center points; calculating the fiber direction angle of each unit under the central point coordinate layer by layer and fiber by fiber, taking the fiber direction angle as the fiber direction angle of the corresponding unit in the current paving layer, and updating the fiber paving layer number of the unit; defining material attributes; defining boundary conditions and loads; generating a calculation input file from the finite element model, and submitting the calculation input file to corresponding finite element software for analysis so as to realize finite element automatic modeling of the variable-stiffness composite laminated plate; according to the method, the problems that variable-angle laying layers need to be input one by one, and efficiency is low are solved, and finite element automatic modeling of the porous structure and the bent plate can be achieved.

Description

Technical Field

[0001] The invention belongs to the field of composite material structure optimization design, and in particular relates to a modeling method of a variable stiffness composite material laminate based on finite elements. Background Art

[0002] Composite materials have become an important part of lightweight design of aviation structures due to their low density, excellent mechanical properties and good designability. Traditional composite laminates adopt the method of laying fibers along a straight line. On the one hand, the stiffness is evenly distributed, resulting in low structural load-bearing efficiency, and on the other hand, it limits the flexible design capability of composite materials. The development of automatic fiber placement technology allows the fiber direction angle to change continuously within the layer, thereby changing the stiffness distribution without changing the thickness of the laminate, realizing the design process of adjusting the fiber path according to the actual load-bearing requirements, greatly improving the material utilization rate and the load-bearing capacity of the structure. Therefore, the research on the mechanical properties of variable stiffness laminates has received more extensive attention, but the design idea of ​​curved fibers increases the design variables, which brings certain difficulties to the modeling and simulation analysis process of variable stiffness laminates. Especially in the finite element analysis commonly used in engineering calculations, the existing simulation software is difficult to realize the automatic input of the continuously changing curved fiber direction angle, resulting in cumbersome operation and low efficiency of the finite element modeling process. In addition, the overlap or gap phenomenon between adjacent fibers that may occur during the curved fiber laying process further increases the difficulty of modeling. Furthermore, compared with simple flat plate structures, the geometric shapes of curved plates and structures with holes commonly found in aircraft structures are more complex, which puts forward higher requirements for automatically constructing corresponding finite element models containing curved fibers. Not only should the changes in fiber orientation angles, fiber overlaps and gaps be considered, but also the reasonable division process of finite element meshes under complex geometric shapes should be considered. These factors hinder the application of variable stiffness laminates containing curved fiber plies in aircraft lightweight design and optimization. Summary of the invention

[0003] In order to solve the above problems, the purpose of the present invention is to provide a modeling method of variable stiffness composite laminates based on finite elements, which can realize automatic finite element modeling of variable stiffness composite laminates, solve the problem that variable angle plies need to be input one by one and the efficiency is low, and can realize automatic finite element modeling of porous structures and curved plates, and at the same time take the influence of fiber gaps and overlaps into modeling considerations, laying a foundation for using variable stiffness composite materials for structural design and mechanical property analysis.

[0004] The present invention is achieved through the following technical solutions:

[0005] A modeling method for variable stiffness composite laminate based on finite element, comprising:

[0006] S1: Define the model parameters of the composite laminate, including the laminate geometry parameters, composite fiber layup parameters and finite element unit size parameters;

[0007] S2: For the composite laminate model defined in S1, a script file is written using a programming language to divide the quadrilateral mesh, obtain the nodes corresponding to each mesh and their geometric coordinates in the O1XYZ rectangular coordinate system, and number all mesh nodes and finite element units;

[0008] S3: Calculate the geometric coordinates of the center point of the finite element unit according to the grid node coordinates obtained in S2;

[0009] S4: Establish a plane rectangular coordinate system ROS, so that the angle between the fiber curve at the origin O and the R axis is T0, and the angle between the R axis and the X axis is

[0010] S5: according to the fiber laying parameters defined in S1, a fiber reference path center curve is defined in the plane rectangular coordinate system ROS established in S4;

[0011] S6: For the geometric coordinates of the unit center point obtained in S3, the fiber direction angle under the coordinates of each unit center point is calculated layer by layer and fiber by fiber, and the angle is used as the fiber direction angle of the corresponding unit in the current layer, and the fiber layer number of the unit is updated at the same time, so as to obtain a finite element unit containing fiber layer information;

[0012] S7: defining material properties for the finite element unit generated in S6 to obtain a finite element model containing material information;

[0013] S8: define boundary conditions and loads for the finite element model generated in S7;

[0014] S9: Generate a calculation input file from the finite element model obtained in S8, and submit it to the corresponding finite element software for finite element analysis to obtain corresponding mechanical properties analysis results.

[0015] Preferably, in S1, the geometric parameters of the laminate include the length and width of the laminate, the curve equation of the curved plate section z=z(x), the coordinates of the center of the circular hole in the O1XYZ rectangular coordinate system, and the radius of the circular hole.

[0016] Preferably, in S1, the fiber laying parameters of the composite material include the fiber path width, the fiber path translation distance, the fiber direction periodic change starting angle T0, the fiber direction periodic change ending angle T1, the single-layer fiber thickness, the fiber laying number, the angle between the fiber curve coordinate system and the laminate geometric coordinate system

[0017] Preferably, in S1, the finite element unit size parameters include mesh size and the number of hole boundary transition unit layers.

[0018] Preferably, in S2, for the quadrilateral meshing process of the curved plate structure, it is projected into a plane along the normal direction. After the plane shape is meshed, coordinate transformation is performed according to the curved plate profile curve equation z=z(x) defined in S1 to calculate the three-dimensional coordinates of the quadrilateral mesh nodes; the transformation process from two-dimensional coordinates to three-dimensional coordinates is defined as follows:

[0019]

[0020] Where x1, y1 represent the X and Y coordinates of the curved plate projected on the XY plane, and z(x1) represents the Z coordinate corresponding to x1 calculated according to the profile curve equation.

[0021] Preferably, in S3, the geometric coordinates of the center point of the finite element unit are calculated by the following formula:

[0022]

[0023] In the formula, x c ,y c 、z c They represent the X, Y, and Z coordinates of the unit center point in the rectangular coordinate system O1XYZ, respectively. i ,y i 、z i (i=1, 2, 3, 4) respectively represent the X, Y, and Z coordinate values ​​of the grid node in the rectangular coordinate system O1XYZ.

[0024] Preferably, in S5, the direction angle of the center curve of the fiber reference path is defined by the following formula:

[0025]

[0026] Where θ represents the direction angle of the center curve of the fiber reference path, r represents the R coordinate value in the rectangular coordinate system ROS, a represents the length of the laminate, T0 represents the starting angle of the periodic change in the fiber direction, T1 represents the ending angle of the periodic change in the fiber direction, and m is the result of rounding up r / (a / 2).

[0027] Preferably, in S5, the fiber reference path center curve is defined by the following formula:

[0028]

[0029] Wherein, θ represents the direction angle of the center curve of the fiber reference path, r represents the R coordinate value in the rectangular coordinate system ROS, s represents the S coordinate value in the rectangular coordinate system ROS, a represents the length of the laminate, T0 represents the starting angle of the periodic change in the fiber direction, T1 represents the ending angle of the periodic change in the fiber direction, and m is the result of rounding up r / (a / 2).

[0030] Preferably, in S6, each fiber is obtained by translating the center curve of the fiber reference path of the layer defined in S5 once or multiple times according to the fiber path translation distance defined in S1.

[0031] Preferably, in S6, during the calculation of the fiber direction angle and the number of plies at the coordinates of each unit center point, it is necessary to perform coordinate transformation between XO1Y and ROS on the unit center point.

[0032] Preferably, in S6, during the calculation of the fiber direction angle and the number of plies at the coordinates of the center point of each unit, for the overlapping parts between the fibers, that is, the common units of multiple fibers, the fiber direction angles of the unit in different fibers are recorded respectively, and the number of plies is accumulated and calculated.

[0033] Preferably, in S6, during the calculation of the fiber direction angle and the number of plies at the coordinates of the center point of each unit, for the gap between the fibers, that is, the unit that does not belong to any fiber in this layer, it is not necessary to record the fiber direction angle of the unit in this layer, and the number of plies remains unchanged.

[0034] Preferably, in S7, the material properties include the elastic modulus E1 along the fiber direction, the elastic modulus E2 perpendicular to the fiber direction, the Poisson's ratio ν 12 , in-plane shear modulus G 12 , fiber density ρ, and the fiber direction angle and fiber ply number described in S6.

[0035] Compared with the prior art, the present invention has the following beneficial technical effects:

[0036] The modeling method of the variable stiffness composite laminate based on finite element disclosed by the present invention defines the model parameters of the composite laminate; uses a programming language to write a script file to divide the quadrilateral grid; calculates the geometric coordinates of the unit center point; calculates the fiber direction angle under the coordinates of each unit center point layer by layer and fiber by fiber, and uses it as the fiber direction angle of the corresponding unit in the current layer, and updates the fiber layer number of the unit; defines material properties; defines boundary conditions and loads; generates a calculation input file for the finite element model, and submits it to the corresponding finite element software for analysis, thereby realizing the finite element automatic modeling of the variable stiffness composite laminate, solving the problem that the variable angle layer needs to be input one by one and the efficiency is low. Compared with the prior art, the present invention can realize the finite element automatic modeling of the hole-containing structure and the curved plate, and at the same time takes the influence of the fiber gap and overlap into the modeling consideration, laying a foundation for the use of variable stiffness composite materials for structural design and mechanical property analysis. BRIEF DESCRIPTION OF THE DRAWINGS

[0037] In order to more clearly illustrate the technical solution of the present invention, the following is a description of the drawings used in the description of the embodiments:

[0038] Figure 1 The flowchart of the modeling method of variable stiffness composite laminate based on finite element;

[0039] Figure 2 This is a schematic diagram of a 1 / 4 model of a composite material curved plate with holes;

[0040] Figure 3 It is a schematic diagram of two-dimensional grid division;

[0041] Figure 4 It is a schematic diagram of three-dimensional mesh division;

[0042] Figure 5 A schematic diagram of the ROS definition of a rectangular coordinate system;

[0043] Figure 6 This is a schematic diagram of the center curve of the first layer fiber placement path;

[0044] Figure 7 Schematic diagram of the center curve of the second layer fiber placement path

[0045] Figure 8 Schematic diagram of the 1 / 4 model unit thickness of a composite material curved plate with holes. DETAILED DESCRIPTION

[0046] The present invention is further explained below with reference to the accompanying drawings and in combination with specific embodiments. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0047] Figure 1 A flow chart showing a modeling method of a variable stiffness composite laminate based on finite elements according to an exemplary embodiment of the present invention is shown. Figure 1 The process shown in the figure is as follows:

[0048] The first step is to define the composite laminate model parameters:

[0049] ① Geometric parameters of laminate: Figure 2 Given a 1 / 4 model of a curved plate with a hole, the side length is 15mm×15mm, the coordinates of the center of the circular hole in the O1XYZ rectangular coordinate system are (0,0,0), the radius is 3mm, and the equation of the profile curve is z=0.02x 2 .

[0050] ② Composite material laying parameters: The number of plies of the laminate is 4, symmetrically distributed, the starting angle T0 of the periodic change of the fiber direction of each layer is 0°, the ending angle T1 of the periodic change of the fiber direction of the first layer is 45°, and the T1 of the second layer is -45°; the thickness of the single fiber layer is 0.13mm, the fiber path width is 2mm, the fiber path translation distance is also 2mm, and the angle between the fiber curve coordinate system and the laminate geometric coordinate system is is 0°.

[0051] ③ Finite element unit size parameters: unit size 0.3mm, circular hole boundary transition 8 layers.

[0052] The second step is to use the programming language to write a script file to divide the quadrilateral mesh, in which the circular hole boundary mesh is radial to ensure the quality of mesh division, such as Figure 3 For the curved plate model, the two-dimensional grid is first divided, and then the curved plate profile curve equation is used. Perform coordinate transformation and calculate the three-dimensional coordinates of the quadrilateral mesh nodes, such as Figure 4 The transformation process from two-dimensional coordinates to three-dimensional coordinates is defined as follows:

[0053]

[0054] In the formula, x1, y1 represent the X and Y coordinates of the node of the curved plate projected on the XY plane, 0.02x1 2 Indicates the Z coordinate corresponding to x1 calculated according to the profile curve equation. So far, the nodes corresponding to each grid and their geometric coordinates in the O1XYZ rectangular coordinate system are obtained, and all grid nodes and units are numbered.

[0055] Step 3: Calculate the geometric coordinates of each unit center point based on the grid node coordinates using the following formula:

[0056]

[0057] In the formula, x c ,y c 、z c They represent the X, Y, and Z coordinates of the unit center point in the rectangular coordinate system O1XYZ, respectively. i ,y i 、z i (i=1, 2, 3, 4) respectively represent the X, Y, and Z coordinate values ​​of the grid node in the rectangular coordinate system O1XYZ.

[0058] Step 4: Establish a plane rectangular coordinate system ROS, whose origin O is the central symmetric point of the center curve of the fiber reference path. The angle between the fiber curve at the origin O and the R axis is T0, and point O coincides with point O1. The angle between the R axis and the X axis is like Figure 5 shown.

[0059] Step 5: First, according to the starting angle T0 and ending angle T1 of the fiber direction periodic change defined in the first step, the direction angle of the center curve of the fiber reference path is defined by the following formula in the plane rectangular coordinate system ROS:

[0060]

[0061] Where θ represents the direction angle of the center curve of the fiber reference path, r represents the R coordinate value in the rectangular coordinate system ROS, a represents the length of the laminate, T0 represents the starting angle of the periodic change in the fiber direction, T1 represents the ending angle of the periodic change in the fiber direction, and m is the result of rounding up r / (a / 2).

[0062] Then, according to the direction angle of the fiber reference path center curve defined above, the fiber reference path center curve is defined by the following formula in the plane rectangular coordinate system ROS:

[0063]

[0064] In the formula, θ represents the direction angle of the fiber reference path center curve, r represents the R coordinate value in the rectangular coordinate system ROS, s represents the S coordinate value in the rectangular coordinate system ROS, a represents the length of the laminate, T0 represents the starting angle of the fiber direction periodic change, T1 represents the ending angle of the fiber direction periodic change, and m is the result of r / (a / 2) rounded up. The first and second layer fiber placement path center curves obtained by translating the fiber reference path center curve are as follows: Figure 6 and 7 shown.

[0065] Step 6: Calculate the fiber direction angle of each unit layer by layer and fiber by fiber, and use it as the fiber direction angle of the corresponding unit in the current layer, and update the fiber layer number of the unit. The golden search method is preferably used here, and any method known to professional technicians to find the maximum value of a function within a given range can also be used. The specific steps are as follows:

[0066] ① Calculate the maximum number of times the center curve of the fiber reference path of this layer is translated in the positive direction of the Y axis.

[0067] The corresponding fiber path A is obtained.

[0068] ② Starting from the above fiber path A, the golden search method is used to determine whether the unit is in the fiber, that is, first transform the unit center coordinates from the XO1Y coordinate system to the ROS coordinate system, and then determine whether the distance between the unit center coordinates and the center curve of the fiber path A is less than or equal to 1 / 2 of the fiber path width. If so, the fiber direction angle under the unit center point coordinates is used as the fiber direction angle of the unit in the current fiber layer, and the fiber layer number of the unit is increased by 1; otherwise, the next unit is calculated.

[0069] ③ After traversing all units, reduce the number of fiber path translations by 1 and repeat ② until all units are no longer in the fiber, completing the calculation of this layer.

[0070] ④ Repeat ① to ③ for each layer to obtain the fiber placement model of the laminate, which includes the effect of fiber overlap, such as Figure 8 As shown in the figure, the darker the grid color in the figure, the more fiber layers there are in the unit, indicating that there is fiber overlap in the unit.

[0071] Step 7: Define material properties, including elastic modulus E1 along the fiber direction, elastic modulus E2 perpendicular to the fiber direction, Poisson's ratio ν 12 , in-plane shear modulus G 12 , fiber density ρ, and the fiber direction angle and number of fiber plies described in step 6.

[0072] Step 8. Define boundary conditions and loads. The boundary conditions include symmetric boundary conditions and fixed constraints, and the loads are forced displacements.

[0073] Step 9: Generate a calculation input file from the obtained finite element model and submit it to the corresponding finite element software for finite element analysis to obtain the corresponding mechanical properties analysis results.

[0074] Through the above discussion and examples, the modeling of variable stiffness composite laminates based on finite element is realized.

[0075] It should be noted that the above is only a part of the implementation of the present invention, and equivalent changes made to the system described in the present invention are all included in the protection scope of the present invention. Those skilled in the art of the present invention can replace the specific examples described in a similar manner, as long as they do not deviate from the structure of the present invention or exceed the scope defined by the claims, they all belong to the protection scope of the present invention.

Claims

1. A modeling method for variable stiffness composite laminates based on finite element, characterized in that: include: S1: Define the model parameters of the composite laminate, including the laminate geometry parameters, composite fiber layup parameters and finite element unit size parameters; S2: For the composite laminate model defined in S1, a script file is written using a programming language to divide the quadrilateral mesh, obtain the nodes corresponding to each mesh and their geometric coordinates in the O1XYZ rectangular coordinate system, and number all mesh nodes and finite element units; S3: Calculate the geometric coordinates of the center point of the finite element unit according to the grid node coordinates obtained in S2; S4: Establish a plane rectangular coordinate system ROS, so that the angle between the fiber curve at the origin O and the R axis is T0, and the angle between the R axis and the X axis is φ; S5: according to the fiber laying parameters defined in S1, a fiber reference path center curve is defined in the plane rectangular coordinate system ROS established in S4; S6: For the geometric coordinates of the unit center point obtained in S3, the fiber direction angle under the coordinates of each unit center point is calculated layer by layer and fiber by fiber, and the angle is used as the fiber direction angle of the corresponding unit in the current layer, and the fiber layer number of the unit is updated at the same time, so as to obtain a finite element unit containing fiber layer information; S7: defining material properties for the finite element unit generated in S6 to obtain a finite element model containing material information; S8: define boundary conditions and loads for the finite element model generated in S7; S9: Generate a calculation input file from the finite element model obtained in S8, and submit it to the corresponding finite element software for finite element analysis to obtain corresponding mechanical properties analysis results.

2. The method according to claim 1, characterized in that: The S1 laminate geometric parameters include the length and width of the laminate, the curve equation of the curved plate section z=z(x), the coordinates of the center of the circular hole in the O1XYZ rectangular coordinate system and the radius of the circular hole.

3. The method according to claim 1, characterized in that: The fiber laying parameters of the S1 composite material include fiber path width, fiber path translation distance, fiber direction periodic change starting angle T0, fiber direction periodic change ending angle T1, single-layer fiber thickness, fiber laying number, and the angle φ between the fiber curve coordinate system and the laminate geometric coordinate system.

4. The method according to claim 1, characterized in that: The S1 finite element unit size parameters include mesh size and the number of hole boundary transition unit layers.

5. The method according to claim 2, characterized in that: In the process of dividing the quadrilateral mesh in S2, for the curved plate structure, it is projected into a plane along the normal direction. After the plane shape is meshed, coordinate transformation is performed according to the curved plate profile curve equation z=z(x) defined in S1 to calculate the three-dimensional coordinates of the quadrilateral mesh nodes; the transformation process from two-dimensional coordinates to three-dimensional coordinates is defined as follows: Where x1, y1 represent the X and Y coordinates of the curved plate projected on the XY plane, and z(x1) represents the Z coordinate corresponding to x1 calculated according to the profile curve equation.

6. The method according to claim 1, characterized in that: The geometric coordinates of the center point of the S3 finite element unit are calculated by the following formula: In the formula, x c ,y c 、z c They represent the X, Y, and Z coordinates of the unit center point in the rectangular coordinate system O1XYZ, respectively. i ,y i 、z i (i=1, 2, 3, 4) respectively represent the X, Y, and Z coordinate values ​​of the grid node in the rectangular coordinate system O1XYZ.

7. The method according to claim 1 or 3, characterized in that: The direction angle of the center curve of the S5 fiber reference path is defined by the following formula: Where θ represents the direction angle of the center curve of the fiber reference path, r represents the R coordinate value in the rectangular coordinate system ROS, a represents the length of the laminate, T0 represents the starting angle of the periodic change in the fiber direction, T1 represents the ending angle of the periodic change in the fiber direction, and m is the result of rounding up r / (a / 2).

8. The method according to claim 1, 3 or 7, characterized in that: The S5 fiber reference path center curve is defined by the following formula: Wherein, θ represents the direction angle of the center curve of the fiber reference path, r represents the R coordinate value in the rectangular coordinate system ROS, s represents the S coordinate value in the rectangular coordinate system ROS, a represents the length of the laminate, T0 represents the starting angle of the periodic change in the fiber direction, T1 represents the ending angle of the periodic change in the fiber direction, and m is the result of rounding up r / (a / 2).

9. The method according to claim 1 or 3 or 7 or 8, characterized in that: Each fiber in S6 is obtained by translating the center curve of the fiber reference path of this layer defined in S5 once or multiple times according to the fiber path translation distance defined in S1.

10. The method according to claim 1, characterized in that: In the calculation process of the fiber direction angle and the number of plies at the coordinates of each unit center point of S6, it is necessary to perform coordinate transformation between XO1Y and ROS on the unit center point.

11. The method according to claim 1, characterized in that: In the process of calculating the fiber direction angle and the number of plies at the coordinates of the center point of each unit of S6, for the overlapping parts between fibers, that is, the common units of multiple fibers, the fiber direction angles of the unit in different fibers are recorded respectively, and the number of plies is accumulated and calculated.

12. The method according to claim 1, characterized in that: In the process of calculating the fiber direction angle and the number of plies at the coordinates of the center point of each unit in S6, for the gap between fibers, that is, the unit that does not belong to any fiber in this layer, there is no need to record the fiber direction angle of the unit in this layer, and the number of plies remains unchanged.

13. The method according to claim 1, characterized in that: The S7 material properties include the elastic modulus E1 along the fiber direction, the elastic modulus E2 perpendicular to the fiber direction, the Poisson's ratio ν 12 , in-plane shear modulus G 12 , fiber density ρ, and the fiber direction angle and fiber ply number described in S6.

Images