Multi-scale variable stiffness optimization method for fiber reinforced composites considering stress constraints

Abstract

The present invention discloses a multi-scale variable stiffness optimization method for fiber-reinforced composite materials considering stress constraints, and relates to the field of multi-scale variable stiffness optimization design for fiber-reinforced composite materials. The present invention adopts a normal distribution fiber optimization interpolation format to interpolate the fiber laying angle, and applies density filtering to the angle design variable, thereby realizing the multi-scale variable stiffness optimization design for fiber-reinforced composite materials, and obtaining the optimization design result with continuous fibers; then based on the Tsai-Wu failure criterion of the composite material, the P-norm method is used to condense the stress constraint function of the large-scale unit, and the explicit sensitivity formula of the stress constraint function to the design variable is derived, thereby realizing the multi-scale variable stiffness optimization design for fiber-reinforced composite materials considering stress constraints. The present invention solves the problem of difficulty in optimization solution of the existing method.

Description

Technical Field

[0001] The present invention relates to the field of multi-scale variable stiffness optimization design of fiber reinforced composite materials, and in particular to a multi-scale variable stiffness optimization method for fiber reinforced composite materials considering stress constraints. Background Art

[0002] Aerospace structures often face extremely harsh working environments during their service. In the past few decades, fiber-reinforced composite materials have been increasingly widely used in the fields of aerospace, automobiles, and new energy equipment due to their high specific strength, high specific stiffness, fatigue resistance, high temperature resistance, and designability. There are more and more optimization design methods for fiber-reinforced composite materials, but stress issues are often not taken into consideration during the design stage. It is often necessary to repeatedly perform finite element analysis and mechanical experiments on composite parts for verification, which is time-consuming and costly. Incorporating structural stress into the optimization design can effectively solve the above problems. The main difficulties in carrying out multi-scale variable stiffness optimization design of fiber-reinforced composites considering stress constraints are as follows: 1) The number of stress constraints will increase with the increase of grid density, and a large number of constraints will make it difficult to solve the optimization algorithm; 2) Discrete structure optimization methods are prone to cause fiber discontinuity, which is not conducive to subsequent production and manufacturing; 3) The mechanical theory of composite materials and its failure criteria are complex, and it is difficult to derive the sensitivity of objective functions and constraint functions; 4) There is a nonlinear relationship between the stress value of the structure and the topological optimization configuration. The stress size of any unit is closely related to the density change of its adjacent area, which will lead to local stress oscillation or non-convergence of the optimization iteration process. Summary of the invention

[0003] In view of the above-mentioned deficiencies in the prior art, the present invention provides a multi-scale variable stiffness optimization method for fiber reinforced composite materials taking into account stress constraints to solve the above-mentioned technical problems.

[0004] In order to achieve the above-mentioned invention object, the technical solution adopted by the present invention is: a multi-scale variable stiffness optimization method of fiber reinforced composite materials considering stress constraints, comprising the following steps:

[0005] S1: Determine the optimization problem, including the design domain of the structure to be optimized, the constraints of the structure, and the load conditions;

[0006] S2: Perform finite element modeling based on the optimization problem and initialize various parameters;

[0007] S3: Perform finite element analysis based on the initialized parameters, including calculation of the overall displacement vector, unit stress and Tsai-Wu failure factor of the structure;

[0008] S4: Calculate the objective function and constraint function based on the finite element analysis results, including the overall compliance, volume and Tsai-Wu failure factor of the structure after P-norm method condensation;

[0009] S5: performing sensitivity analysis based on the objective function and the constraint function, including calculating the sensitivity information of the objective function to the topological design variables and the angle design variables, and the sensitivity information of the constraint function to the topological design variables and the angle design variables;

[0010] S6: input the current topological design variables, angle design variables, objective function, constraint function and sensitivity information into the MMA optimization solver for optimization and solution to obtain new design variables;

[0011] S7: Repeat the above steps S3-S6 until the new design variables reach the convergence standard or the maximum number of iterations, output the optimized design results and the Tsai-Wu failure factor cloud diagram, and complete the multi-scale variable stiffness optimization of fiber reinforced composite materials considering stress constraints.

[0012] Furthermore, S3 includes the following sub-steps:

[0013] S31: According to the first-order shear deformation beam theory, the stiffness matrix of the unit is obtained, and the stiffness matrix is ​​grouped to obtain the overall stiffness matrix of the composite laminate, and the overall displacement vector of the structure is obtained by solving the equilibrium equation;

[0014] S32: Calculate the element stress at the Gaussian integration point based on the overall displacement vector of the structure;

[0015] S33: Based on the unit stress at the Gaussian integration point, the unit stress in the main direction of the material is calculated through the coordinate transformation tensor, the unit stress in the main direction of the material is interpolated, and the Tsai-Wu failure factor is obtained.

[0016] Further, in said S31, The stiffness matrix of the element for:

[0017]

[0018] in, For the The in-plane part of the element stiffness matrix, and For the The in-plane and bending coupled parts of the element stiffness matrix, For the The bending part of the element stiffness matrix, For the The shear part of the element stiffness matrix;

[0019]

[0020]

[0021]

[0022]

[0023]

[0024] in, is the number of laminate layers considered, is the area of ​​the finite element, For the The in-plane strain-displacement matrix of the element, For the The bending strain-displacement matrix of the unit, For the The shear strain-displacement matrix of the unit, To use SIMP method to interpolate Unit No. The result of the layer plane elasticity tensor is, To use SIMP method to interpolate Unit No. The result for the layer shear elastic tensor is, For laminate The coordinates of the top of the layer, For laminate The coordinates of the bottom of the layer;

[0025]

[0026]

[0027]

[0028] in, and are coordinates, is the shape function;

[0029] Perform stiffness matrix grouping to obtain the overall stiffness matrix of the composite laminate :

[0030]

[0031] in, is the number of the element stiffness matrix, is the total number of element stiffness matrices;

[0032] The overall displacement vector of the structure is obtained by solving the equilibrium equation :

[0033]

[0034] in, is the structural load vector.

[0035] Furthermore, the unit stress at the Gaussian integration point in S32 is:

[0036]

[0037] in, For the Unit No. Layer plane stress vector, For the Unit No. Layer shear stress vector, For the Unit No. Layer plane strain vector, For the Unit No. Layer plane strain vector;

[0038]

[0039]

[0040] in, For the The displacement vector of the unit, For the The thickness of the layer.

[0041] Furthermore, the unit stress in the main direction of the material is calculated by the coordinate transformation tensor in S33 as follows:

[0042]

[0043] in, For the Unit No. plane stress vectors of the layer elements in the principal directions of the material, For the Unit No. Shear stress vector of the layer in the principal direction of the material, For the Unit No. The transposed matrix of the intra-plane part, For the Unit No. The transposed matrix of the layer clipping part;

[0044] Interpolate element stresses in the principal directions of the material:

[0045]

[0046] in, For the Unit No. The principal plane stress vector after layer interpolation, After interpolation Unit No. Layer principal shear stress vector, For the The actual topological design variables of the unit, is the stress interpolation coefficient;

[0047] Tsai-Wu failure factor satisfy:

[0048]

[0049]

[0050]

[0051]

[0052]

[0053]

[0054]

[0055] in, , , , , and is the tensor coefficient, For the Unit No. The magnitude of the stress along the fiber direction of the layer, For the Unit No. The magnitude of the vertical stress in the fiber layer, For the Unit No. The magnitude of the shear stress in the interlaminar plane is is the longitudinal tensile strength of the material, is the longitudinal compressive strength of the material, is the transverse tensile strength of the material, is the transverse compressive strength of the material, is the shear strength of the material.

[0056] Furthermore, the overall compliance of the structure in S4 for:

[0057]

[0058] Among them, the superscript is the transpose symbol;

[0059] Volume constraint function for:

[0060]

[0061] in, is the area of ​​the unit;

[0062] Tsai-Wu failure factor after P-norm method aggregation for:

[0063]

[0064] in, is the P-norm aggregation coefficient, For the Unit No. Tsai-Wu failure factor of the layer.

[0065] Furthermore, the sensitivity information of the objective function in S5 to the topological design variables and the angle design variables is:

[0066] Overall flexibility of the structure Design variables for real topologies The sensitivity information is:

[0067]

[0068] in, For the The unit displacement extraction vector of unit number, For the The element stiffness matrix of the element number, For the The face-to-face stiffness matrix of the element, For the The surface-bending stiffness matrix of the element, For the The bending-surface stiffness matrix of the element, For the The bending-bending stiffness matrix of the element, For the The shear stiffness matrix of the element;

[0069]

[0070]

[0071]

[0072]

[0073]

[0074] in, For the The in-plane displacement-strain matrix of the element, For the The bending displacement-strain matrix of the element, For the The shear displacement-strain matrix of the element, To use SIMP method to interpolate Unit No. The result for the plane elasticity tensor of the layer is, To use SIMP method to interpolate Unit No. Results for the layer shear elastic tensor;

[0075] Overall flexibility of the structure Design variables for angles The sensitivity information is:

[0076]

[0077] in, After filtering Unit No. Layer angle design variables, is the filter radius, express Unit No. The distance between the units, express Unit No. The distance between the units, For the The unit number set within the unit filtering radius.

[0078] Furthermore, the sensitivity information of the constraint function in S5 to the topological design variables and the angle design variables is:

[0079] Volume constraint function Topology Design Variables The sensitivity information is:

[0080]

[0081] Volume constraint function Design variables for angles The sensitivity information is 0;

[0082] Tsai-Wu failure factor Topology Design Variables The sensitivity information is:

[0083]

[0084] in, For the Unit No. Tsai-Wu failure factor of the layer, For the Unit No. The stress vector after layer relaxation, , and is the sensitivity coefficient;

[0085]

[0086]

[0087]

[0088] in, For the The actual topological design variables of the unit, For the Unit No. The coordinate transposed matrix of the layer, For the Number unit displacement vector;

[0089] Tsai-Wu failure factor Design variables for angles The sensitivity information is:

[0090]

[0091] in, , and is the sensitivity coefficient;

[0092]

[0093]

[0094] .

[0095] The beneficial effects of the present invention are:

[0096] (1) The present invention uses the first-order shear deformation theory, which can be extended to realize the multi-scale variable stiffness optimization design of fiber-reinforced composite laminates;

[0097] (2) The present invention uses the Tsai-Wu failure criterion, which is more accurate in failure assessment of fiber-reinforced composite materials.

[0098] (3) The present invention applies density filtering to angle design variables to obtain a more continuous fiber angle design;

[0099] (4) The present invention provides the derivation of explicit sensitivity formulas for structural flexibility, volume constraints, and stress constraints obtained by P-norm condensation to topological design variables and angle design variables. BRIEF DESCRIPTION OF THE DRAWINGS

[0100] Figure 1 Flowchart of the multi-scale variable stiffness optimization method for fiber reinforced composites considering stress constraints.

[0101] Figure 2 Schematic diagram of the optimization design domain, loads, and boundary constraints for an L-beam.

[0102] Figure 3 Schematic diagram of quadrilateral unit in natural coordinate system.

[0103] Figure 4 Topological results and fiber orientation diagram for the optimal design of an L-beam without considering stress constraints.

[0104] Figure 5 Tsai-Wu failure factor cloud diagram for the optimal design of L-beam without considering stress constraints.

[0105] Figure 6 Optimize the design topology and fiber orientation for an L-beam considering stress constraints.

[0106] Figure 7 Tsai-Wu failure factor cloud diagram for optimal design of L-beam considering stress constraints. DETAILED DESCRIPTION

[0107] The present invention will be further described below in conjunction with the accompanying drawings and specific embodiments.

[0108] The present invention is based on the first order shear deformation theory, uses the Normal Distribution Fiber Optimization (NDFO) interpolation format to interpolate the fiber laying angle, and applies the density filter to the angle design variable, realizing the multi-scale variable stiffness optimization design of fiber reinforced composite materials, and obtaining the optimization design result with continuous fibers. Subsequently, based on the Tsai-Wu failure criterion of composite materials, the P-norm method is used to condense the large-scale unit stress constraint function, and the explicit sensitivity formula of the stress constraint function to the design variable is derived, realizing the multi-scale variable stiffness optimization design of fiber reinforced composite materials considering stress constraints.

[0109] like Figure 1 As shown, a multi-scale variable stiffness optimization method for fiber reinforced composite materials considering stress constraints includes the following steps:

[0110] S1: Determine the optimization problem, including the design domain of the structure to be optimized, the constraints of the structure, and the load conditions;

[0111] S2: Perform finite element modeling based on the optimization problem and initialize various parameters;

[0112] S3: Perform finite element analysis based on the initialized parameters, including calculation of the overall displacement vector, unit stress and Tsai-Wu failure factor of the structure;

[0113] S4: Calculate the objective function and constraint function based on the finite element analysis results, including the overall compliance, volume and Tsai-Wu failure factor of the structure after P-norm method condensation;

[0114] S5: performing sensitivity analysis based on the objective function and the constraint function, including calculating the sensitivity information of the objective function to the topological design variables and the angle design variables, and the sensitivity information of the constraint function to the topological design variables and the angle design variables;

[0115] S6: input the current topological design variables, angle design variables, objective function, constraint function and sensitivity information into the MMA optimization solver for optimization and solution to obtain new design variables;

[0116] S7: Repeat the above steps S3-S6 until the new design variables reach the convergence standard or the maximum number of iterations, output the optimized design results and the Tsai-Wu failure factor cloud diagram, and complete the multi-scale variable stiffness optimization of fiber reinforced composite materials considering stress constraints.

[0117] In one embodiment of the present invention, an optimization problem is determined: taking an L-shaped beam as an example, the effectiveness of the present invention is demonstrated. Figure 2 The geometry, loads and boundary constraints of the L-beam are given, and the entire area occupied by the L-beam is a designable domain. The geometric dimensions are mainly based on the dimensions that may appear in reality and can be manufactured by general printers. In this example, the length and width of 100mm*100mm are taken as an example. In order to distinguish whether the stress constraint is effective, the load size cannot be too small (so that the structure will not fail without applying stress constraints) or too large (so that the structure that meets the conditions does not exist and the optimization result cannot be obtained). In this example, the 6 nodes on the right are taken, and a downward force of 100N is applied to each node, with a total load of 600N.

[0118] Finite element modeling is performed based on the optimization problem, and various parameters are initialized: first, by establishing constant parameters in MATLAB, the length and width of the rectangular space occupied by the L-shaped beam are set, and the geometric dimensions of the structure are stored; then, by establishing a series of matrices in MATLAB, the rectangular space is meshed, and the horizontal and vertical coordinates of the center of the unit and the node numbers of the four nodes corresponding to the unit and the coordinate information of all nodes are stored according to the unit number and node number, respectively, to realize the finite element meshing of the rectangular space; then, according to the geometric dimensions of the L-shaped beam, the units outside the space occupied by the L-shaped beam are selected in the rectangular space, the unit numbers of these units are stored in a new vector and marked as non-design domains, and the rest are marked as design domains.

[0119] Then, the constant parameters are set, and the initial value of all topological design variables is 0.5; the initial value of the angle design variable is the value when the angle value is 0°; the constant parameter in the SIMP method is set to 3; the parameter in the NDFO method is set to 4 for all units initially, and then it is judged every two steps whether the weight corresponding to the angle design variable is close to 1. If the weight corresponding to the angle design variable is greater than 0.95, the NDFO parameter of the unit is reduced to 97% of the original; the filter radius in the density filtering is set to the length of 3 units; eta in the Heaviside projection is set to 0.5, beta is set to 1, and the increasing rule of beta in the iteration process is set to double every 50 steps; the number of alternative fiber angles is set to 16; the small value to prevent the overall stiffness matrix from being singular is set to 1E-4; the limit value of the volume constraint is set to 0.4; and the constant parameter of stress interpolation is set to 0.8.

[0120] Next, we start with the assembly of the stiffness matrix by interpolating the angles of the underlying composite finite element calculation. Tier The fiber laying angles in the design domain occupied by the finite element units are consistent. The NDFO interpolation format uses coefficients right Alternative fiber placement angles Perform weighted summation to obtain the equivalent fiber laying angle :

[0121]

[0122] Among them, the equivalent fiber laying angle The subscript of , and the alternative fiber placement angle The subscript of This is because the fiber laying angles of different units and layers are different. Each unit and each layer has its own equivalent fiber laying angle. , and this angle is used as the actual (Physical) fiber laying angle in the finite element analysis. represents the discrete alternative fiber placement angles considered in this study, which are independent of the number of units and layers, so the subscript only uses In the above formula is the normalized weight coefficient:

[0123]

[0124] in, It is about the angle design variable The normal distribution function of Indicates the number of alternative fiber placement angles:

[0125]

[0126] Among them, the superscript Indicates that the coefficient is used to interpolate the angle; Indicates layer, Indicates Unit No. Indicates No. Alternative fiber layup angle; The value range is ; yes The angle value corresponding to the alternative fiber laying angle; Indicates Tier The NDFO angle penalty coefficient of the unit, which changes as the optimization progresses.

[0127] The multi-scale variable stiffness topology optimization of fiber-reinforced composite materials has the problems of discontinuous fibers in the optimization results and difficulty in manufacturing. To solve the above problems, the density filtering method proposed by Bruns and Tortorelli et al. and mathematically proved by Bourdin is adopted. The formula is as follows:

[0128]

[0129] According to the first-order shear deformation theory, the stress-strain relationship of an orthotropic material in the main direction of the material under plane stress state can be expressed as:

[0130]

[0131] in, is the stress tensor, is the strain tensor, and are the bending stiffness tensor and the shear stiffness tensor, respectively. For simplicity, the subscripts in the stiffness tensor are contracted by the following index correspondence rule:

[0132] 1-11

[0133] 2-22

[0134] 6-12

[0135] 4-23

[0136] 5-13

[0137] So the first line can be expressed as:

[0138]

[0139] The second line is represented as:

[0140]

[0141] Subscript 11 indicates the direction along the fiber (called axial, longitudinal or first principal direction), subscript 22 indicates the direction within the plane of the single-layer board and perpendicular to the fiber (called transverse or second principal direction), and subscript 12 indicates the 12th principal direction plane, which is the plane where the single-layer board is located.

[0142] Engineering elastic constants are often used in engineering to represent the elastic properties of materials. These engineering elastic constants include the generalized elastic modulus , Poisson's ratio and shear modulus , these constants can be determined by simple tensile and shear tests. The elements in the stiffness tensor can be expressed in terms of engineering elastic constants as follows:

[0143]

[0144] The above calculations can be used to obtain the stiffness tensor in the material coordinate system. However, the fiber laying direction of fiber reinforced composite materials is usually different from the global coordinate system in finite element analysis. The fiber laying direction is at an angle to the Axis Angle (that is, the fiber laying angle is ), by transforming the coordinates of the tensor The fiber laying angle can be calculated as The stiffness tensor at and for:

[0145]

[0146] The coordinate transformation tensor and They are:

[0147]

[0148] Because different units and layers have different fiber laying angles , so the coordinate transformation tensor , and the transformed stiffness tensor Subscripts are added to indicate differences between units and layers.

[0149] The material interpolation topology optimization based on artificial density has mesh dependence and checkerboard phenomenon. To solve the above problems, the density filtering method proposed by Bruns and Tortorelli et al. and mathematically proved by Bourdin is adopted. The formula is as follows:

[0150]

[0151] In order to obtain clear (0-1, discretized) topology optimization results, the improved Heaviside projection is used to project the filtered artificial density to obtain the real artificial density. , as shown below:

[0152]

[0153] SIMP interpolation format through real artificial density Optimize the topology. When Fill the design domain with material; When The design domain where the number element is located is not filled with material. The SIMP interpolation format is as follows:

[0154]

[0155] in, To use SIMP method to interpolate Unit No. The result of the layer plane elasticity tensor is, To use SIMP method to interpolate Unit No. The result for the layer shear elastic tensor is, Small values ​​introduced to prevent the stiffness matrix from becoming singular.

[0156] The S3 includes the following sub-steps:

[0157] S31: According to the first-order shear deformation beam theory, the stiffness matrix of the unit is obtained, and the stiffness matrix is ​​grouped to obtain the overall stiffness matrix of the composite laminate, and the overall displacement vector of the structure is obtained by solving the equilibrium equation;

[0158] S32: Calculate the element stress at the Gaussian integration point based on the overall displacement vector of the structure;

[0159] S33: Based on the unit stress at the Gaussian integration point, the unit stress in the main direction of the material is calculated through the coordinate transformation tensor, the unit stress in the main direction of the material is interpolated, and the Tsai-Wu failure factor is obtained.

[0160] In the S31 The stiffness matrix of the element for:

[0161]

[0162] in, For the The in-plane part of the element stiffness matrix, and For the The in-plane and bending coupled parts of the element stiffness matrix, For the The bending part of the element stiffness matrix, For the The shear part of the element stiffness matrix;

[0163]

[0164]

[0165]

[0166]

[0167]

[0168] in, is the number of laminate layers considered, is the area of ​​the finite element, For the The in-plane strain-displacement matrix of the element, For the The bending strain-displacement matrix of the unit, For the The shear strain-displacement matrix of the unit, To use SIMP method to interpolate Unit No. The result of the layer plane elasticity tensor is, To use SIMP method to interpolate Unit No. The result for the layer shear elastic tensor is, To prevent small values ​​introduced by the stiffness matrix singularity, For laminate The coordinates of the top of the layer, For laminate The coordinates of the bottom of the layer;

[0169]

[0170]

[0171]

[0172] in, and are coordinates, is the shape function. This embodiment adopts a second-order 4-node rectangular element, such as Figure 3 The shape function is expressed as:

[0173]

[0174] in, and is the unit coordinate;

[0175] Perform stiffness matrix grouping to obtain the overall stiffness matrix of the composite laminate :

[0176]

[0177] in, is the number of the element stiffness matrix, is the total number of element stiffness matrices;

[0178] The overall displacement vector of the structure is obtained by solving the equilibrium equation :

[0179]

[0180] in, is the structural load vector.

[0181] After the node displacement is obtained by finite element analysis, the strain of the unit at the Gaussian integration point can be calculated by the node displacement, and then the stress of the unit at the Gaussian integration point can be calculated, as shown in the following formula:

[0182]

[0183] in, is the strain in the x direction, is the strain in the y direction, is the shear strain in the xy plane, is the shear strain in the xz plane, is the shear strain in the yz plane, , , is the displacement in the x, y, and z directions;

[0184] The element stress at the Gaussian integration point in S32 is:

[0185]

[0186] in, For the Unit No. Layer plane stress vector, For the Unit No. Layer shear stress vector, For the Unit No. Layer plane strain vector, For the Unit No. Layer plane strain vector;

[0187]

[0188]

[0189] in, For the The displacement vector of the unit, For the The thickness of the layer.

[0190] use and The finite element analysis obtains stress and strain in the global coordinate system. When solving the Tsai-Wu failure coefficient, the coordinate transformation tensor needs to be converted into stress in the main direction of the material. The stress state at the Gaussian integration point of each layer of the unit in the global coordinate system is obtained. Finally, in order to calculate the Tsai-Wu failure coefficient, the coordinate transformation tensor is used Compute the stresses in the principal directions of the element material:

[0191] The unit stress in the main direction of the material is calculated by the coordinate transformation tensor in S33 as follows:

[0192]

[0193] in, For the Unit No. plane stress vectors of the layer elements in the principal directions of the material, For the Unit No. Shear stress vector of the layer in the principal direction of the material, For the Unit No. The transposed matrix of the intra-plane part, For the Unit No. The transposed matrix of the layer clipping part;

[0194] Interpolate element stresses in the principal directions of the material:

[0195]

[0196] in, For the Unit No. The principal plane stress vector after layer interpolation, After interpolation Unit No. Layer principal shear stress vector, For the The actual topological design variables of the unit, is the stress interpolation coefficient, usually taken as 0.5;

[0197] Fiber reinforced composite materials are orthotropic materials. When the external load acts along the main direction of the material, it becomes the main direction load, and the corresponding stress is called the main direction stress. If the load direction is inconsistent with the main direction of the material, the stress in the load direction can be converted into the stress in the main direction of the material through coordinate transformation.

[0198] For isotropic materials, the maximum stress and maximum strain referred to in various strength theories refer to the principal stress and principal strain of the material; but for anisotropic materials, since the maximum applied stress does not necessarily correspond to the dangerous state of the material (for example, when the stress exceeding the maximum transverse tensile stress acts in the longitudinal direction, the material is not in a dangerous state at this time because the main direction of the material is reinforced by fibers), the maximum principal stress that is independent of the material direction is meaningless.

[0199] There are five basic strengths of orthotropic materials:

[0200] , ——Longitudinal tensile and compressive strength (along the main direction of the material 1)

[0201] , ——Transverse tensile and compressive strength (along the main direction of the material 2)

[0202] ——Shear strength (along 1-2 plane)

[0203] Most experimental determinations of material strength are based on unidirectional stress states, but actual structural problems often involve plane stress states or spatial stress states. Here we discuss the strength of orthotropic single-layer materials, so we involve plane strength theory. It is assumed that the material is macroscopically uniform and some microscopic failure mechanisms are not considered.

[0204] Tsai-Wu strength theory states that the following inequality holds true when the structure fails:

[0205] Tsai-Wu failure factor satisfy:

[0206]

[0207]

[0208]

[0209]

[0210]

[0211]

[0212]

[0213] in, , , , , and is the tensor coefficient, For the Unit No. The magnitude of the stress along the fiber direction of the layer, For the Unit No. The magnitude of the vertical stress in the fiber layer, For the Unit No. The magnitude of the shear stress in the interlaminar plane is is the longitudinal tensile strength of the material, is the longitudinal compressive strength of the material, is the transverse tensile strength of the material, is the transverse compressive strength of the material, is the shear strength of the material.

[0214] This embodiment calculates the Tsai-Wu failure factors at the midpoints of different layers and then averages them to evaluate whether structural failure occurs in each layer of the unit.

[0215] After calculating the Tsai-Wu failure factors at the four Gaussian integration points of each unit, because establishing constraints on the Tsai-Wu failure factors of each Gaussian integration point will bring a huge amount of stress constraints, making the problem difficult to solve. In order to reduce the number of constraints, this scheme uses the P-norm method to condense stress constraints.

[0216] The overall compliance of the S4 structure for:

[0217]

[0218] Among them, the superscript is the transpose symbol;

[0219] Volume constraint function for:

[0220]

[0221] in, is the area of ​​the unit;

[0222] Tsai-Wu failure factor after P-norm method aggregation for:

[0223]

[0224] in, is the P-norm aggregation coefficient, For the Unit No. Tsai-Wu failure factor of the layer.

[0225] The sensitivity information of the objective function in S5 to the topological design variables and the angle design variables is:

[0226] Overall flexibility of the structure Design variables for real topologies The sensitivity information is:

[0227] The derivatives of both sides of the equation with respect to the filtered topology design variables are:

[0228]

[0229] For the structural static equilibrium equation ( )The derivatives of both sides are:

[0230]

[0231] After transformation, we get:

[0232]

[0233] Substituting into the above formula, we can get:

[0234]

[0235] because It is obtained by expanding the stiffness matrix of each unit, so the above formula is expanded as follows:

[0236]

[0237] in, For the The unit displacement extraction vector of unit number, For the The element stiffness matrix of the element number, For the The face-to-face stiffness matrix of the element, For the The surface-bending stiffness matrix of the element, For the The bending-surface stiffness matrix of the element, For the The bending-bending stiffness matrix of the element, For the The shear stiffness matrix of the element;

[0238]

[0239]

[0240]

[0241]

[0242]

[0243] in, For the The in-plane displacement-strain matrix of the element, For the The bending displacement-strain matrix of the element, For the The shear displacement-strain matrix of the element, To use SIMP method to interpolate Unit No. The result for the plane elasticity tensor of the layer is, To use SIMP method to interpolate Unit No. Results for the layer shear elastic tensor;

[0244] and Topology Design Variables The partial derivative of can be calculated by the following formula:

[0245]

[0246]

[0247]

[0248] in, and For the Unit No. The plane elasticity tensor and shear elasticity tensor of the layer, is the Heaviside projection coefficient, The density filtration method is used to process the Number of unit topology design variables.

[0249] Overall flexibility of the structure Design variables for angles The sensitivity information is:

[0250] Design variables for angles The sensitivity solution of the first half is the same as . and right The partial derivative of can be calculated by the following formula:

[0251]

[0252] The following is simplified for brevity:

[0253]

[0254]

[0255]

[0256]

[0257] in yes right The derivative of :

[0258]

[0259] is a micro-design variable The function of It can be expressed as follows:

[0260]

[0261]

[0262] in, For the Unit No. Tier The normalized NDFO weight coefficients corresponding to the alternative fiber placement angles are: For the Unit No. Tier NDFO weight coefficients corresponding to the alternative fiber laying angles;

[0263]

[0264]

[0265] in, For the Unit No. Angular design variables after layer density filtering, is the filter radius, express Unit No. The distance between the units, express Unit No. The distance between the units, For the The unit number set within the unit filtering radius.

[0266] The sensitivity information of the constraint function in S5 to the topological design variables and the angle design variables is:

[0267] Volume constraint function Topology Design Variables The sensitivity information is:

[0268]

[0269] Volume constraint function Design variables for angles The sensitivity information is 0;

[0270] Tsai-Wu failure factor Topology Design Variables The sensitivity information is:

[0271]

[0272] in, For the Unit No. Tsai-Wu failure factor of the layer, For the Unit No. The principal direction stress vectors of the layer, , and is the sensitivity coefficient;

[0273]

[0274]

[0275]

[0276] in, For the The actual unit topology design variables, For the Unit No. The coordinate transposed matrix of the layer, For the The unit displacement vector of the unit number;

[0277] Direct solution It will greatly increase the computational consumption. To reduce the amount of calculation, the adjoint vector method is used as follows:

[0278]

[0279] make:

[0280]

[0281] in It can be obtained by solving the following linear equations:

[0282]

[0283]

[0284] in, is the adjoint vector parameter, For the The coordinate transformation matrix of the unit number;

[0285] Tsai-Wu failure factor Design variables for angles The sensitivity information is:

[0286]

[0287] in, , and is the sensitivity coefficient;

[0288]

[0289]

[0290] .

[0291] right The calculation of is also done using the adjoint vector method as follows:

[0292]

[0293] make:

[0294]

[0295] Observation can reveal .

[0296] After obtaining the sensitivity information, it is introduced into the MMA algorithm to obtain a new round of design variables, and then steps S3-S6 are repeated until the change of the design variable is less than 1% of the value range or when the number of optimization iterations reaches 500. The optimization algorithm stops iterating and performs post-processing of the results, outputting the optimized design results and the Tsai-Wu failure factor cloud diagram.

[0297] In one embodiment of the present invention, when the stress constraint mentioned in the above method is not introduced, the optimization result obtained is as follows: Figure 4 As shown, we can see that the topological structure of the optimization result has right angles. Figure 5It can be found that the structure has obvious stress concentration at the right angle, where the Tsai-Wu failure factor exceeds 1, and the structure fails. After introducing stress constraints, the optimization results are as follows: Figure 6 As shown, we can see that the right angle is replaced in the topological configuration. Figure 7 It can be found that the stress concentration phenomenon has disappeared. From the data in the color scale bar on the left, it can be seen that the maximum Tsai-Wu failure factor of the unit does not exceed 1, and the structure has not failed.

[0298] In summary, the multi-scale variable stiffness optimization design method for fiber reinforced composite materials considering stress constraints proposed in the present invention is effective. It is worth noting that the method is not limited to the given L-beam example, but is also applicable to the optimization design of other situations.

[0299] Those skilled in the art will appreciate that the embodiments described herein are intended to help readers understand the principles of the present invention, and should be understood that the protection scope of the present invention is not limited to such specific statements and embodiments. Those skilled in the art can make various other specific variations and combinations that do not deviate from the essence of the present invention based on the technical revelations disclosed in the present invention, and these variations and combinations are still within the protection scope of the invention.

Claims

1. A multi-scale variable stiffness optimization method for fiber reinforced composite materials considering stress constraints, characterized in that: The following steps are involved: S1: Determine the optimization problem, including the design domain of the structure to be optimized, the constraints of the structure, and the load conditions; S2: Perform finite element modeling based on the optimization problem and initialize various parameters; S3: Perform finite element analysis based on the initialized parameters, including calculation of the overall displacement vector, unit stress and Tsai-Wu failure factor of the structure; S4: Calculate the objective function and constraint function based on the finite element analysis results, including the overall compliance, volume and Tsai-Wu failure factor of the structure after P-norm method condensation; S5: performing sensitivity analysis based on the objective function and the constraint function, including calculating the sensitivity information of the objective function to the topological design variables and the angle design variables, and the sensitivity information of the constraint function to the topological design variables and the angle design variables; S6: input the current topological design variables, angle design variables, objective function, constraint function and sensitivity information into the MMA optimization solver for optimization and solution to obtain new design variables; S7: Repeat the above steps S3-S6 until the new design variable reaches the convergence standard or the maximum number of iterations, output the optimization design results and the Tsai-Wu failure factor cloud map, and complete the multi-scale variable stiffness optimization of fiber reinforced composite materials considering stress constraints; The S3 includes: S31: According to the first-order shear deformation beam theory, the stiffness matrix of the unit is obtained, and the stiffness matrix is ​​grouped to obtain the overall stiffness matrix of the composite laminate, and the overall displacement vector of the structure is obtained by solving the equilibrium equation; In the S31 The stiffness matrix of the element for: in, For the The in-plane part of the element stiffness matrix, and For the The in-plane and bending coupled parts of the element stiffness matrix, For the The bending part of the element stiffness matrix, For the The shear part of the element stiffness matrix; in, is the number of laminate layers considered, is the area of ​​the finite element, For the The in-plane strain-displacement matrix of the element, For the The bending strain-displacement matrix of the unit, For the The shear strain-displacement matrix of the unit, To use SIMP method to interpolate Unit No. The result of the layer plane elasticity tensor is, To use SIMP method to interpolate Unit No. The result for the layer shear elastic tensor is, For laminate The coordinates of the top of the layer, For laminate Coordinates of the bottom of the layer, superscript is the transpose symbol; Perform stiffness matrix grouping to obtain the overall stiffness matrix of the composite laminate : in, is the number of the element stiffness matrix, is the total number of element stiffness matrices; The overall displacement vector of the structure is obtained by solving the equilibrium equation : in, is the structural load vector.

2. The multi-scale variable stiffness optimization method for fiber reinforced composite materials considering stress constraints according to claim 1, characterized in that: The S3 also includes: S32: Calculate the element stress at the Gaussian integration point based on the overall displacement vector of the structure; S33: Based on the unit stress at the Gaussian integration point, the unit stress in the main direction of the material is calculated through the coordinate transformation tensor, the unit stress in the main direction of the material is interpolated, and the Tsai-Wu failure factor is obtained.

3. The multi-scale variable stiffness optimization method for fiber reinforced composite materials considering stress constraints according to claim 2, characterized in that: In S31 , and for: in, and is the coordinate, is the shape function; The element stress at the Gaussian integration point in S32 is: in, For the Unit No. Layer plane stress vector, For the Unit No. Layer shear stress vector, For the Unit No. Layer plane strain vector, For the Unit No. Layer plane strain vector; in, For the The displacement vector of the unit, For the The thickness of the layer.

4. The multi-scale variable stiffness optimization method for fiber reinforced composite materials considering stress constraints according to claim 3, characterized in that: The unit stress in the main direction of the material is calculated by the coordinate transformation tensor in S33 as follows: in, For the Unit No. plane stress vectors of the layer elements in the principal directions of the material, For the Unit No. Shear stress vector of the layer in the principal direction of the material, For the Unit No. The transposed matrix of the intra-plane part, For the Unit No. The transposed matrix of the layer clipping part; Interpolate element stresses in the principal directions of the material: in, For the Unit No. The principal plane stress vector after layer interpolation, After interpolation Unit No. Layer principal shear stress vector, For the The actual topological design variables of the unit, is the stress interpolation coefficient; Tsai-Wu failure factor satisfy: in, , , , , and is the tensor coefficient, For the Unit No. The magnitude of the stress along the fiber direction of the layer, For the Unit No. The magnitude of the vertical stress in the fiber layer, For the Unit No. The magnitude of the shear stress in the interlaminar plane is is the longitudinal tensile strength of the material, is the longitudinal compressive strength of the material, is the transverse tensile strength of the material, is the transverse compressive strength of the material, is the shear strength of the material.

5. The multi-scale variable stiffness optimization method for fiber reinforced composite materials considering stress constraints according to claim 4, characterized in that: The overall compliance of the S4 structure for: Among them, the superscript is the transpose symbol; Volume constraint function for: in, is the area of ​​the unit; Tsai-Wu failure factor after P-norm method aggregation for: in, is the P-norm aggregation coefficient, For the Unit No. Tsai-Wu failure factor of the layer.

6. The multi-scale variable stiffness optimization method for fiber reinforced composite materials considering stress constraints according to claim 5, characterized in that: The sensitivity information of the objective function in S5 to the topological design variables and the angle design variables is: Overall flexibility of the structure Design variables for real topologies The sensitivity information is: in, For the The unit displacement extraction vector of unit number, For the The element stiffness matrix of the element number, For the The face-to-face stiffness matrix of the element, For the The surface-bending stiffness matrix of the element, For the The bending-surface stiffness matrix of the element, For the The bending-bending stiffness matrix of the element, For the The shear stiffness matrix of the element; in, For the The in-plane displacement-strain matrix of the element, For the The bending displacement-strain matrix of the element, For the The shear displacement-strain matrix of the element, To use SIMP method to interpolate Unit No. The result for the plane elasticity tensor of the layer is, To use SIMP method to interpolate Unit No. Results for the layer shear elastic tensor; Overall flexibility of the structure Design variables for angles The sensitivity information is: in, After filtering Unit No. Layer angle design variables, is the filter radius, express Unit No. The distance between the units, express Unit No. The distance between the units, For the The unit number set within the unit filtering radius.

7. The multi-scale variable stiffness optimization method for fiber reinforced composite materials considering stress constraints according to claim 6, characterized in that: The sensitivity information of the constraint function in S5 to the topological design variables and the angle design variables is: Volume constraint function Topology Design Variables The sensitivity information is: Volume constraint function Design variables for angles The sensitivity information is 0; Tsai-Wu failure factor Topology Design Variables The sensitivity information is: in, For the Unit No. Tsai-Wu failure factor of the layer, For the Unit No. The stress vector after layer relaxation, , and is the sensitivity coefficient; in, For the The actual topological design variables of the unit, For the Unit No. The coordinate transposed matrix of the layer, For the Number unit displacement vector; Tsai-Wu failure factor Design variables for angles The sensitivity information is: in, , and is the sensitivity coefficient; 。

Images