Mesoscale calculation method for ABD matrix of braided composite materials and computer-readable storage medium

Abstract

The present invention provides a mesoscale calculation method of the ABD matrix of a woven composite material and a computer-readable storage medium, comprising the following steps: generating a mesoscale RVE for a woven composite material through composite material modeling software; and importing the generated mesoscale RVE into ABAQUS finite element analysis software; running an ABD matrix plug-in in the ABAQUS finite element analysis software, inputting material parameters in the ABD matrix plug-in, and specifying material directions for yarns; the material parameters include yarn parameters and matrix parameters; setting periodic boundary conditions through the ABD matrix plug-in to apply 6 periodic boundary conditions to the mesoscale RVE, thereby realizing 6 variations of the mesoscale RVE; analyzing the variation data of the 6 variations of the mesoscale RVE through the ABAQUS finite element analysis software, thereby calculating the ABD matrix of the thin woven composite material.

Description

Technical Field

[0001] The present invention relates to the technical field of composite material performance analysis, and in particular to a mesoscale calculation method of an ABD matrix of a woven composite material and a computer-readable storage medium. Background Art

[0002] Composite materials are becoming increasingly popular in the design of various high-performance structures, such as aerospace, automotive, and civil engineering components, due to their remarkable mechanical properties, such as high stiffness and strength. Unidirectional laminates and woven composites are two commonly used fiber-reinforced composites, each with its own unique characteristics and advantages. Unidirectional laminates consist of multiple layers of fibers, with all fibers in each layer oriented in the same direction and embedded in a matrix material. Woven composites consist of fiber yarns that are crisscrossed to form a multi-directional, braided structure. Woven composites offer several advantages over unidirectional laminates, including enhanced mechanical isotropy, improved impact resistance, and superior damage tolerance. However, these benefits come with trade-offs, such as increased manufacturing complexity and higher cost. Modeling the stiffness behavior of woven composites is also a considerably more complex task compared to unidirectional laminates.

[0003] The ABD matrix is ​​a stiffness matrix used to describe the elastic properties of composite laminates and is the basis for the mechanical analysis of laminate structures. The classical laminate theory has been proven to accurately predict the ABD matrix of unidirectionally reinforced carbon fiber composites, but there are significant errors in woven composites. A multiscale homogenization method based on mesoscale RVE and FEA has gradually become the mainstream method for ABD matrix prediction. After experimental verification, this method can accurately predict the ABD matrix of woven composites. However, the process of implementing this method involves many complex links: data transmission, FEA, post-processing, etc., making it difficult to directly apply to any form of woven composite materials. Summary of the Invention

[0004] In order to overcome the above technical defects, the purpose of the present invention is to provide a mesoscale calculation method of the ABD matrix of a woven composite material and a computer-readable storage medium, which are used to simplify the calculation process of the ABD matrix of the woven composite material.

[0005] The present invention discloses a mesoscale calculation method for predicting the ABD matrix of a thin woven composite material, comprising the following steps: generating a mesoscale RVE for a woven composite material through composite material modeling software; and importing the generated mesoscale RVE into ABAQUS finite element analysis software; running an ABD matrix plug-in in the ABAQUS finite element analysis software, inputting material parameters in the ABD matrix plug-in, and specifying material directions for yarns; the material parameters include yarn parameters and matrix parameters; setting periodic boundary conditions through the ABD matrix plug-in to apply 6 periodic boundary conditions to the mesoscale RVE, thereby realizing 6 variations of the mesoscale RVE; analyzing the variation data of the 6 variations of the mesoscale RVE through the ABAQUS finite element analysis software, thereby calculating the ABD matrix of the thin woven composite material.

[0006] Preferably, the generation of a mesoscale RVE for woven composite materials by composite material modeling software includes: obtaining the gap size S between adjacent warp yarns and adjacent weft yarns, the length L, width W and height H of the mesoscale RVE, thereby obtaining the weaving form of the woven composite material; selecting a periodic mesoscale RVE according to the weaving form; taking a cross-sectional micrograph of the woven composite material, and determining the cross-sectional shape, size and yarn undulation path of each yarn according to the cross-sectional micrograph, thereby obtaining the true geometric characteristics of the yarn; adjusting the mesoscale RVE according to the true geometric characteristics of the yarn to obtain a mesoscale RVE model close to the actual situation.

[0007] Preferably, running the ABD matrix plug-in in the ABAQUS finite element analysis software also includes: the ABD matrix plug-in deletes the default settings of the mesoscale RVE, and the default settings include material parameters of the yarn and matrix, solver settings, constraint equations, and boundary conditions.

[0008] Preferably, running the ABD matrix plug-in in the ABAQUS finite element analysis software, inputting material parameters in the ABD matrix plug-in, and specifying material direction for the yarn also includes: selecting a finite element analysis calculation solver corresponding to the grid division form of the mesoscale RVE in the ABD matrix plug-in; the grid division form includes volume grid, voxel grid and dry grid.

[0009] Preferably, when the dry grid form is used, the adhesive behavior between the yarn surfaces is approximated by setting the viscous contact properties of the yarn surfaces, including: defining the relationship between the stress and separation displacement between the two adhesive surfaces as:

[0010] Where t and δ represent stress and separation displacement; subscript n represents the normal direction, s and t represent the two transverse shear directions; knn represents the normal stiffness, k ss and k tt represents the tangential stiffness, k nn 、k ss and k tt Collectively referred to as uncoupled traction stiffness; the stiffness value is set by the uncoupled traction stiffness; the stiffness value between the nodes in contact with each other on the surfaces of two yarns is used to create an association relationship between the nodes in contact with each other, thereby approximating the adhesive behavior between the yarn surfaces.

[0011] Preferably, the six periodic boundary conditions are applied to the mesoscale RVE to realize the six variations of the mesoscale RVE, including: establishing a number of reference points RP on the surface around the mesoscale RVE, the several reference points RP forming the mid-surface edge of the mesoscale RVE; the unit node of the upper surface edge in the periodic grid with the reference point RP as the base point is defined as the upper node, and the unit node of the lower surface edge in the periodic grid with the reference point RP as the base point is defined as the lower node; the reference point RP, the upper node and the lower node in the same periodic grid are rigidly linked using MPCbeam units to realize coupling of the nodes and all degrees of freedom of RP in all directions; six periodic boundary conditions are applied to the several reference points RP respectively to perform deformation loading on the mid-surface of the mesoscale RVE to realize deformation control of the entire mesoscale RVE.

[0012] Preferably, the six periodic boundary conditions are applied to the mesoscale RVE to achieve the six variations of the mesoscale RVE, and further comprises: if the upper node / the lower node is offset from the surrounding surfaces of the mesoscale RVE, then: the coordinates of the two vertices of the mid-body diagonal of the mesoscale RVE (x min ,y min ,z min ) and (x max ,y max ,z max ) (thereby obtaining the coordinates of the remaining six vertices of the mesoscale RVE); a rectangular parallelepiped region is established by the ABAQUS finite element analysis software, and the coordinates of the two vertices of the mid-body diagonal of the rectangular parallelepiped region are (x min -t x ,y min -t y ,z min -t z ) and (x max +t x ,y max +t y ,z max +t z), determining the boundary of the rectangular parallelepiped region according to the coordinates of two vertices of the mid-body diagonal line of the rectangular parallelepiped region; all grid nodes in the rectangular parallelepiped region form a first set; Where n is the tolerance value; all mesh nodes in the first set are sorted according to the size of the x-coordinate / y-coordinate, the upper and lower endpoints are found among the nodes with the same x-coordinate / y-coordinate, the midpoint of the upper and lower endpoints is set as the reference point RP, and the nodes with the same x-coordinate / y-coordinate are rigidly linked to the reference point RP using MPC beam units.

[0013] Preferably, the applying six periodic boundary conditions to the mesoscale RVE to achieve six variations of the mesoscale RVE includes applying displacement boundary conditions as a percentage of unit strain.

[0014] Preferably, the composite material modeling software includes Texge and Digmat.

[0015] The present invention also discloses a computer-readable storage medium having a computer program stored thereon. When the computer program is executed by a processor, the steps of the mesoscale calculation method for predicting the ABD matrix of a thin woven composite material are implemented.

[0016] Compared with the existing technology, the above technical solution has the following beneficial effects:

[0017] 1. The present invention installs a customized ABD matrix plug-in in the ABAQUS finite element analysis software. This allows the plug-in to automatically set periodic boundary conditions, sequentially apply six predefined loading conditions, and calculate the ABD matrix of the laminate as the final output. The present invention is simple to operate and the overall method is highly integrated, eliminating the need for further programming and analysis by the user, thereby saving time and effort. It also avoids the user from having to repeat tedious operations in multi-scale modeling, numerical calculations, and data processing. BRIEF DESCRIPTION OF THE DRAWINGS

[0018] Figure 1 It is a plain woven composite laminate and homogenization model;

[0019] Figure 2 are the cross-sectional shapes before deformation (ABCD) and after deformation (A'B'C'D') at any position in the homogenized model of the laminate;

[0020] Figure 3 Create a process for mesoscale RVE of double-layer plain woven composite laminates;

[0021] Figure 4 Three forms of meshing for mesoscale RVE finite element models;

[0022] Figure 5 Schematic diagram of periodic grid;

[0023] Figure 6 Schematic diagram of applying PBC to mesoscale RVE of woven composite laminates;

[0024] Figure 7 Schematic diagram of six deformations of the mid-surface of a mesoscale RVE under periodic boundary conditions;

[0025] Figure 8 Modeling six deformations of mesoscale RVE under periodic boundary conditions;

[0026] Figure 9 is the displacement-force relationship in the x direction of the CP-X reference point;

[0027] Figure 10 Numbering conventions for mesoscale RVE yarns in different types of woven composites;

[0028] Figure 11 Execute the flow chart for AMWC;

[0029] Figure 12 It is the lateral surface phenomenon of the mesoscale RVE in the vertex offset;

[0030] Figure 13 This is the method for displaying vertices that deviate from the plane. DETAILED DESCRIPTION

[0031] The advantages of the present invention are further described below with reference to the accompanying drawings and specific embodiments.

[0032] Exemplary embodiments will be described in detail herein, with examples illustrated in the accompanying drawings. In the following description, when referring to the drawings, identical numerals in different figures represent identical or similar elements, unless otherwise indicated. The embodiments described in the following exemplary embodiments are not intended to represent all possible embodiments consistent with the present disclosure. Rather, they are merely examples of apparatus and methods consistent with certain aspects of the present disclosure, as detailed in the appended claims.

[0033] The terms used in this disclosure are for the purpose of describing specific embodiments only and are not intended to limit the disclosure. As used in this disclosure and the appended claims, the singular forms "a," "an," "the," and "the" are intended to include the plural forms as well, unless the context clearly indicates otherwise. It should also be understood that the term "and / or" as used herein refers to and encompasses any and all possible combinations of one or more of the associated listed items.

[0034] It should be understood that although the terms first, second, third, etc. may be used in this disclosure to describe various information, such information should not be limited to these terms. These terms are only used to distinguish information of the same type from each other. For example, without departing from the scope of this disclosure, first information may also be referred to as second information, and similarly, second information may also be referred to as first information. Depending on the context, the word "if" as used herein may be interpreted as "at the time of" or "when" or "in response to determining."

[0035] In the description of the present invention, it should be understood that the terms "longitudinal", "transverse", "up", "down", "front", "back", "left", "right", "vertical", "horizontal", "top", "bottom", "inside", "outside", etc., indicating the orientation or position relationship, are based on the orientation or position relationship shown in the accompanying drawings, and are only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the device or element referred to must have a specific orientation, be constructed and operated in a specific orientation, and therefore cannot be understood as limiting the present invention.

[0036] In the description of the present invention, unless otherwise specified and limited, it should be noted that the terms "installed", "connected" and "connected" should be understood in a broad sense. For example, it can be a mechanical connection or an electrical connection, or it can be the internal communication between two components. It can be a direct connection or an indirect connection through an intermediate medium. For ordinary technicians in this field, the specific meanings of the above terms can be understood according to the specific circumstances.

[0037] In the following description, the suffixes such as "module", "component" or "unit" used to represent elements are only used to facilitate the description of the present invention and have no specific meaning. Therefore, "module" and "component" can be used interchangeably.

[0038] First, the ABD matrix is ​​explained. The ABD matrix is ​​a basic concept in the classical laminate theory (CLT) to characterize the effective stiffness characteristics of composite laminates. The ABD matrix is ​​a 6×6 matrix derived based on the CLPT. It characterizes the relationship between the load, strain and curvature of the composite laminate and is used to define the elastic properties of the entire laminate.

[0039] CLPT assumes that the laminate material is uniformly distributed and the interlaminar displacement is continuous. The laminate homogenization model is as follows: Figure 1 The deformation of the section at any position in the homogenized model is shown as Figure 2 This figure shows the cross-sectional shape of a section perpendicular to the y-axis of the laminate before and after deformation, where the xy plane is equidistant from the upper and lower surfaces of the laminate and is called the mid-plane.

[0040] Select a point C on the mid-surface (point C' after deformation), and the displacement function in the x, y, and z directions can be expressed as a function of the coordinates x and y. The specific expression is as follows:

[0041] u o =u o (x,y),v o =v o (x,y),w=f(x,y)(1-1)

[0042] The deflection angles of the mid-plane around the y-axis and around the x-axis are:

[0043]

[0044] The distance between point B and point C in the cross section is zB. The displacement of this point (B') after deformation is:

[0045]

[0046]

[0047] The displacement expression after deformation of any point in the cross section is:

[0048]

[0049] Where z is the distance between any point in the cross section and the mid-surface before deformation.

[0050] According to the linear strain-displacement relationship under the small displacement assumption, the strain expression at any point in the cross section is as follows:

[0051]

[0052]

[0053] If the point is on the mid-plane, the strain expression is:

[0054]

[0055] The curvature expressions of the laminate in three directions are:

[0056]

[0057] In most cases, a single-layer plate is considered as the basic building block of a laminate structure and is made of orthotropic materials. The mechanical properties are expressed using engineering constants (E x ,E y ,E z ,G xy ,G xz ,G yz ,ν xy ,νxz ,ν yz ) is used to describe the thickness of the single-layer plate (z direction) is very small compared to the dimensions in other plane directions (x, y), so it is approximately considered that σ z =τ xz =τ yz =0, which defines the stress-strain relationship of the single-layer plate as shown in formula (1-8).

[0058]

[0059] where the matrix is the transformation matrix of the two-dimensional stiffness matrix Q, which is obtained from the stress rotation axis formula. The specific transformation relationship is shown in formula (1-9):

[0060]

[0061] Where θ is the angle between the material principal direction 1 and the x-axis in the xy coordinate system. The relationship between each item in the Q matrix and the engineering constants of the single-layer board is shown in formula (1-10):

[0062]

[0063] Assumptions Figure 1 The laminate shown is composed of n layers of single-layer plates. The strain-strain relationship of the kth layer can be obtained by combining formulas (1-5), (1-6), (1-7) and (1-8):

[0064]

[0065] Figure 2 Shows the force diagram of the laminate, where N x , N y The unit length force in the x and y directions, N xy is the shear force per unit length, and the dimensions of the three in mm are N / mm. x , M y are the unit length bending moments of the laminate around the y and x axes, M xy is the torque per unit length, and the dimension of the three in mm is N. The above terms can be obtained by integrating the stress on each single layer along the thickness of the laminate:

[0066]

[0067] Substituting formula (1-11) into formula (1-12) yields the final form of the ABD matrix, as shown in formula (1-13):

[0068]

[0069] In the formula It can be seen that formula (1-13) is a symmetric matrix, where the sub-matrices [A] and [D] represent the in-plane tensile and out-of-plane bending stiffness characteristics of the laminate, respectively, and [B] represents the coupling between in-plane and out-of-plane loads and deformations. ij =A ji ,B ij =B ji ,D ij =D ji (i,j=1,2,6). A ij It is called in-plane (tensile or shear) stiffness, B ij It is called coupling (tensile and bending coupling) stiffness, D ij It is called the out-of-plane (bending or torsional) stiffness.

[0070]

[0071] The present invention installs a plug-in for the ABD matrix of woven composites (ABD Matrix of Woven Composites, AMWC) in the ABAQUS finite element analysis software to simplify the calculation process of the ABD matrix of woven composites. In the calculation of the ABD matrix of the composite material of the present invention, the deformation of each entry of the ABD matrix is ​​calculated by considering the deformation of the mesoscale RVE under six different loading conditions - three in-plane and three out-of-plane conditions, thereby obtaining the individual values ​​of the ABD matrix of the composite material. Compared with the existing calculation method in the ABAQUS finite element analysis software, the calculation process of the present invention is simpler and more automated.

[0072] The calculation process of the ABD matrix of the present invention is a calculation method based on the mesoscale RVE. The process of obtaining the ABD matrix in the mesoscale link is also called the homogenization method based on the periodic mesoscale RVE. This homogenization method first requires creating a mesoscale RVE at the mesoscale of the woven composite laminate, then applying the PBC on the mesoscale RVE, and finally obtaining the ABD matrix of the laminate through six deformation loadings. The reasonable selection of the mesoscale RVE, accurate modeling and correct application of the PBC become the key factors determining whether the homogenization is successful.

[0073] The present invention discloses a mesoscale calculation method for predicting the ABD matrix of thin woven composite materials (see Appendix Figure 11 ), including the above-mentioned mesoscale RVE generation and selection process, modeling process, PBC application process and ABD matrix calculation process, as follows:

[0074] S100, selection and modeling process of mesoscale RVE: generating a mesoscale RVE for weaving a composite material using composite material modeling software; and importing the generated mesoscale RVE into ABAQUS finite element analysis software;

[0075] S200, a process of importing a mesoscale RVE into the ABAQUS finite element analysis software: running an ABD matrix plug-in in the ABAQUS finite element analysis software, inputting material parameters in the ABD matrix plug-in, and specifying a material direction for the yarn; the material parameters include yarn parameters and matrix parameters;

[0076] S300, PBC application process: setting periodic boundary conditions through the ABD matrix plug-in to apply six periodic boundary conditions to the mesoscale RVE, thereby realizing six variations of the mesoscale RVE;

[0077] S400, ABD matrix calculation process: using the ABAQUS finite element analysis software to analyze the deformation data of the six variations of the mesoscale RVE, thereby calculating and obtaining the ABD matrix of the thin woven composite material.

[0078] For the mesoscale RVE in step S100, in composite material mechanics, the mesoscale RVE is also called REV (Representative elementary volume) or UC (Unit cell), which represents the overall structure by selecting a minimum volume. Measuring or analyzing this minimum volume structure can accurately obtain the response of the overall structure, which requires that the volume of the mesoscale RVE must be large enough to contain as many microstructural features as possible to obtain the correct macroscopic response, but also small enough to reduce the analysis cost. Figure 3 Taking a double-layer woven composite material as an example, the process of selecting and establishing a mesoscale RVE and the important geometric parameters that need to be considered are demonstrated.

[0079] In order to establish a correct mesoscale RVE finite element model, it is generally necessary to use a scanning electron microscope to obtain the real microstructural characteristics of the woven composite material. Figure 3As shown, the first step is to select a suitable periodic mesoscale RVE according to the weaving form. The geometric parameters that need to be determined in this process are: the size S of the gap between adjacent warp yarns and adjacent weft yarns, the length L, width W and height H of the mesoscale RVE. The second step is to take a cross-sectional micrograph of the woven composite laminate (the cross section of the yarn is usually assumed to be an ideal ellipse, with a and b representing its major axis and minor axis, respectively. The fluctuation of each yarn is determined by changing the position of the control point (the red dot in the figure)). Based on the micrograph, the cross-sectional shape, size and yarn undulation path of each yarn are determined, and the mesoscale RVE is adjusted again according to the actual geometric characteristics of the yarn, and finally a mesoscale RVE model close to the actual situation is obtained.

[0080] Preferably, the present invention uses Texgen, a woven composite material modeling software, to model the mesoscale RVE geometry. The mesoscale RVE geometry model created by this software is subsequently transferred to ABAQUS as a .inp file for analysis. However, it should be noted that the modeling software of the present invention is not limited to Texgen; Digmat or user-defined models are also acceptable. It can be understood that the mesoscale RVE geometry model adapted by the plug-in of the present invention is applicable not only to woven composite materials, but also to any thin plate model that satisfies the Kirchochoff assumption and mesh periodicity conditions.

[0081] Regarding the process of importing the mesoscale RVE into the ABAQUS finite element analysis software in step S200, when importing the mesoscale RVE into the ABAQUS finite element analysis software, the directly imported model file will be accompanied by some default settings, including the material parameters of the yarn and the matrix, solver settings, constraint equations, boundary conditions, etc. The ABD matrix plug-in AMWC will first delete all these settings.

[0082] And it is necessary to select the corresponding ABAQUS analysis and calculation method (solution) according to the grid form of the mesoscale RVE built, see the attached Figure 4 , usually including three commonly used grid division forms: volume grid (Volume mesh), voxel grid (Volxe mesh) and dry grid (Drymesh).

[0083] like Figure 4As shown in (a), the first mesh division form is volume mesh. When using this division method, the mesoscale RVE finite element model generated by Texgen is composed of yarn and matrix, and the unit type is C3D4 / C3D10 tetrahedral unit. The matrix and the adjacent surface mesh of the yarn share the same nodes. The mesoscale RVE finite element model generated by this division method has fewer units and high computational efficiency. The premise for the correct application of PBC to the mesoscale RVE is that the mesoscale RVE grid is periodic, that is, the number of nodes on the relative edges / surfaces in the mesoscale RVE is consistent and the nodes completely overlap after being projected to the opposite edges / surfaces, such as Figure 5 shown.

[0084] As attached Figure 4 As shown in (b), the second meshing method is the voxel mesh. Similar to the volume mesh, the mesoscale RVE finite element model generated using this method consists of yarns and a matrix, using C3D8 / C3D10 hexahedral elements. Adjacent matrix and yarn meshes share common nodes. Unlike the volume mesh, all elements generated by the volxe mesh are cube meshes of identical size, or voxel meshes. Using this mesh, the volxe mesh can divide the mesoscale RVE into a periodic grid, regardless of the complexity of the weave.

[0085] As attached Figure 4 As shown in (c), the third mesh division form is dry mesh. There are only yarns in the mesoscale RVE, and no resin matrix is ​​included. The unit type is C3D8 / C3D10 hexahedral unit. In actual situations, there is relative sliding between the yarn surfaces, and this sliding is affected by the viscous behavior of the matrix. In order to correctly capture this behavior, it is crucial to accurately simulate the structural response of the mesoscale RVE. Volume mesh and Volxe mesh solve this problem well by performing mesh common node processing between the yarn and the matrix, but in Drymesh, it is necessary to set the yarn surface viscous contact properties (*CONTACT—Cohesive Behavior) to approximate this adhesive behavior. Specific implementation method: The stiffness value between the nodes in contact with each other on the two yarn surfaces is used to create an association relationship between the nodes in contact with each other, thereby approximating the adhesive behavior between the yarn surfaces. The setting of this stiffness value needs to be determined by the cohesive traction separation law, which defines the stress and separation displacement relationship between the two adhesive surfaces: Where t and δ represent stress and separation displacement; subscript n represents the normal direction, s and t represent the two transverse shear directions; k nn represents the normal stiffness, kss and k tt represents the tangential stiffness, k nn 、k ss and k tt They are collectively referred to as uncoupled traction stiffness; that is, the stiffness value is set by the uncoupled traction stiffness.

[0086] You also need to enter material parameters (including yarn and matrix parameters) into the ABD matrix plug-in AMWC and specify the material orientation for the yarns. When generating the mesoscale RVE, Texgen assigns a name and number to each yarn according to a specific pattern, and this pattern is retained after importing into ABAQUS. Each yarn is placed in a differently numbered Set, and the ABD matrix plug-in AMWC assigns the correct material orientation to the yarns according to the numbering pattern.

[0087] Attachment Figure 10 The rules for numbering the RVE yarns in different types of braided composites are shown. Figure 10 (a) The left image shows a three-layer 2D braided composite mesoscale RVE. The right image shows the numbering pattern for the first layer of yarns: the six yarns in the y-direction (the yarn axis is perpendicular to the y-axis) are first numbered counterclockwise, yarn 0 to yarn 5. The four yarns in the x-direction are then numbered counterclockwise, yarn 6 to yarn 9. The second and third layers are then numbered, starting with yarn 10, following the same numbering pattern as the first layer. The GUI input parameters for this process are: Number of yarns (x): 4; Number of yarns (y): 6; Number of layers: 3. The ABD matrix plugin AMWC creates two collections named Allyarn_x and Allyarn_y, placing the x- and y-direction yarns in each collection: Allyarn_y = [yarn0-yarn5 yarn10-yarn15 yarn20-yarn25], and Allyarn_x = [yarn6-yarn9yarn16-yarn19yarn26-yarn29]. Then, select the Allyarn_y collection to align the material's principal direction with the y-axis, and select the Allyarn_x collection to align the material's principal direction with the x-axis.

[0088] Attachment Figure 10(b) The mesoscale RVE and yarn numbering scheme for a 3D woven composite. In AMWC, a 3D woven composite is considered one layer. The order is to first number all yarns in the y direction, yarn0-yarn6, and then all yarns in the x direction, yarn7-yarn15. The corresponding GUI input parameters for this process are: Number of yarns (x): 9; Number of yarns (y): 7; Number of layers: 1. The corresponding total set of yarns in the y and x directions is Allyarn_y = [yarn0-yarn6], and Allyarn_x = [yarn7-yarn15].

[0089] Attachment Figure 10 (c) is a mesoscale RVE of a 3D braided composite material containing half a yarn. Two half yarns are considered as one yarn when numbering. The rest of the rules are the same as those in the attached Figure 10 (b) Consistent.

[0090] Regarding the PBC application process in step S300, according to homogenization theory, the entire woven composite single-layer plate / laminate can be formed by repeatedly forming multiple mesoscale RVEs. Regardless of the position of the mesoscale RVE within the plate, the strain and curvature of each mesoscale RVE should exhibit a consistent response. This requirement is achieved by applying a specific PBC to the mesoscale RVEs. PBC, also known as an equal displacement boundary condition, aims to equalize the displacements of symmetrical nodes on two opposing surfaces of the mesoscale RVE, thereby ensuring that the opposing surfaces of the mesoscale RVE remain parallel before and after deformation.

[0091] Attachment Figure 6 Taking a mesoscale RVE of a plain woven composite single-layer plate structure as an example, the process of applying PBC is demonstrated. First, several reference points RP are established on the surfaces around the mesoscale RVE. The reference points RP form the mid-surface edges of the mesoscale RVE. The number of RPs is consistent with the number of unit nodes on the upper / lower surface edges of the mesoscale RVE (periodic grid). The unit nodes of the upper surface edge in the periodic grid with the reference point RP as the base point are defined as upper nodes, and the unit nodes of the lower surface edge in the periodic grid with the reference point RP as the base point are defined as lower nodes. The reference point RP, the upper node, and the lower node in the same periodic grid are rigidly linked using MPC beam units to achieve coupling of all degrees of freedom of the nodes and RPs in all directions.

[0092] The top view of the mesoscale RVE mid-surface and the distribution of reference point RP are shown in the attached figure. Figure 6(b) Six periodic boundary conditions are then applied to the reference points RP to deform the mid-surface of the mesoscale RVE and control the deformation of the entire mesoscale RVE. The mid-surface reference points have the following relative displacement relationships in different directions:

[0093]

[0094] The superscript and They represent the jth pair of reference points on the x-direction opposite sides (Left side and Right side) of the mesoscale RVE midplane, see Figure 6 (b) purple dot; and They represent the jth pair of reference points on the y-direction opposite sides (Bottomside and Topside) of the mesoscale RVE midplane, see Figure 6 Yellow point in (b). u, v, w are the displacements of a reference point in the x, y, and z directions, respectively. x θ is the angle between the mid-plane at a reference point and the horizontal plane where the y-axis is located after the reference point rotates around the x-axis. y It is the angle between the mid-plane at a reference point and the horizontal plane where the x-axis is located after the reference point is rotated about the y-axis.

[0095] Formula (1-15) defines the loading of the PBC and also needs to be implemented with the degree of freedom control equation (such as the *EQUATION command) and boundary conditions (*BOUNDARY command) in ABAQUS. For example, the degree of freedom association equation of the first line of the PBC formula (1-15) in ABAQUS is: in and Same as above, CP-X represents a reference point established at any location outside the mesoscale RVE model (see Figure 6 (a) red dot). CP-X represents the x-direction translational freedom at point CP-X. a, b and c are displacement amplitude coefficients. Given the coefficient values ​​and u CP-X Specific values ​​can be achieved The displacement relationship between the three reference points CP-X. For example, let a=1, b=-1, c=-1, u CP-X =ε x L, means reference point and The relative displacement between them in the x direction is equal to ε xWhen full PBC is applied, a pair of mid-surface reference points has six constraint equations (each reference point has six degrees of freedom). All reference points on the edge of the RL face are associated with CP-X, and all reference points on the edge of the TB face are associated with CP-Y. The constraint equations for each degree of freedom are modified accordingly according to formula (1-15). If the RL face has N pairs of reference points on the edge and the TB face has M pairs of reference points on the edge, the total number of equations is 6 (N + M).

[0096] After the PBC is successfully applied, the deformation loading is performed on the mid-surface of the mesoscale RVE to achieve the following Figure 7 The six deformations are shown. The black surface represents the mid-surface before deformation, and the blue surface represents the mid-surface after deformation. The deformation of the entire mesoscale RVE driven by the mid-surface can be seen. Figure 8 shown.

[0097] In the ABD matrix calculation process of step S400, the ABAQUS finite element analysis software is mainly used to analyze the variant data of the six variants of the mesoscale RVE, thereby calculating the ABD matrix of the thin woven composite material. Specifically, the mesoscale RVE is firstly subjected to the ABAQUS finite element analysis software (e.g., using a linear analysis method). Figure 8 The finite element analysis corresponding to the six deformations shown above can obtain the ABD matrix by extracting the data of the two reference points CP-X and CP-Y and performing post-processing analysis. Figure 8 The tensile deformation in (a) is used as an example to illustrate this process. Figure 8 As shown in (a), when the mesoscale RVE is in the x-direction stretching deformation, that is, ε x ≠0,ε y =ε xy =κ x =κ y =κ xy =0, from formula (1-8):

[0098] N x =A 11 ε x ,N y =A 21 ε y ,N xy =A 61 ε xy ,M x =B 11 κ x ,M y =B 12 κ y ,M xy =B 16 κ xy (1-17)

[0099] According to the homogenization theory:

[0100]

[0101]

[0102] in is the resultant force in the x direction of all nodes (n) on the surface of the mesoscale RVE x=L, is the resultant force in the y direction of all nodes (m) on the surface y=W, is the resultant force in the y direction of all nodes (n) on the surface x=L, is the resultant moment of all nodes (n) on the surface of the mesoscale RVE x=L in the x direction (bending around the y axis), is the resultant moment in the y direction (bending around the x axis) of all nodes (m) on the mesoscale RVE y=W surface, is the resultant moment in the y direction (bending around the x-axis) of all nodes (n) on the surface of the mesoscale RVE x=L.

[0103] Since all nodes on the mesoscale RVE surface have been coupled to the mid-surface reference point, and the mid-surface reference point is also coupled to CP-X or CP-Y, after performing finite element analysis on the mesoscale RVE, the support reaction force / support reaction moment on CP-X or CP-Y is extracted to obtain the surface resultant force / resultant moment in formula (1-18). In order to obtain the first term A in (1-17), 11 For example, the tensile strain ε in the mesoscale RVEx direction x =1( Figure 8 (a) deformation), and substituting it into the first term of formula (1-15), we can know that the stretching distance between the left and right surfaces of the mesoscale RVE in the x direction is L, and the displacement-force relationship in the reference point CP-Xx direction is extracted, as follows: Figure 9 As shown, where F x (L,y,z) is the surface force x=L Substitute into formula (1-18) to get N x . Finally, N x and ε x =1 into (1-17) to determine A 11 . According to this method, the remaining stiffness terms in (1-17) are solved one by one, and all the elements in the first column of the ABD matrix are determined. The remaining columns in the ABD matrix can be obtained from the corresponding deformation results (the second column corresponds to Figure 8 (b), the third column corresponds to Figure 8 (c), the fourth column corresponds to Figure 8 (d), the fifth column corresponds to Figure 8 (e), the sixth column corresponds to Figure 8(f)), the post-processing method is consistent with the stiffness terms in the first column, and the ABD matrix of the composite material represented by the REV is calculated.

[0104] As mentioned above, in the process of applying PBC, it is necessary to couple the nodes with the same x, y coordinates on the side of the mesoscale RVE with the intermediate reference point by MPC beam (see Appendix Figure 6 ), but Texgen will inevitably fail in the modeling process. Figure 12 The mesh error shown is that the nodes are deviated from the plane ( Figure 12 In the case of XZ plane), although the offset distance is so small that it is difficult to detect with the naked eye, it will still cause the node to be out of alignment with the mid-surface reference point ( Figure 12 The phenomenon of missing MPC beam between the yellow dots in the middle.

[0105] Therefore, the present invention introduces the Tolerance variable (tolerance value), see Figure 13 Taking the left surface in the y direction as an example, the coordinates (x min ,y min ,z min ) and (x max ,y max ,z max ), so that the coordinates of the remaining 6 vertices of the mesoscale RVE can be obtained.

[0106] Then, a rectangular area (rectangular black frame) is established by the ABAQUS finite element analysis software. The coordinates of the two vertices (red triangle and blue triangle) of the mid-body diagonal of the rectangular area are (x min -t x ,y min -t y ,z min -t z ) and (x max +t x ,y max +t y ,z max +t z ), determine the boundary of the rectangular area according to the coordinates of the two vertices of the mid-body diagonal line of the rectangular area; all grid nodes in the rectangular area form a first set, named left_plane_nodes.

[0107] Theoretically iIf the (i=x,y,z) value is greater than the maximum value of the node offset plane, all mesoscale RVE side surface nodes including the offset nodes can be selected. i It should be smaller than the minimum side length of the unit in the mesoscale RVE to prevent the rectangular box selection area from including points near the surface ( Figure 13 The yellow node in the enlarged box in the figure). Therefore, t i The specific value of is determined by the following formula: Where n is the tolerance value.

[0108] To avoid non-convergence issues that may occur during the finite element analysis, it is preferable to apply displacement boundary conditions as a percentage of unit strain. Because all finite element analyses in AMWC are geometrically linear, the magnitude of the applied deformation strain does not affect the final ABD matrix results.

[0109] The present invention also discloses a computer-readable storage medium having a computer program stored thereon. When the computer program is executed by a processor, the steps of the mesoscale calculation method for predicting the ABD matrix of a thin woven composite material are implemented.

[0110] It should be noted that the embodiments of the present invention have better practicability and do not impose any form of limitation on the present invention. Any technician familiar with the field may use the technical content disclosed above to change or modify it into an equivalent effective embodiment. However, any modification or equivalent changes and modifications made to the above embodiments based on the technical essence of the present invention without departing from the content of the technical solution of the present invention are still within the scope of the technical solution of the present invention.

Claims

1. A mesoscale calculation method for predicting the ABD matrix of thin woven composite materials, characterized in that: The steps include: Generating a mesoscale RVE for braiding a composite material using composite material modeling software; and importing the generated mesoscale RVE into ABAQUS finite element analysis software; Running an ABD matrix plug-in in the ABAQUS finite element analysis software, inputting material parameters in the ABD matrix plug-in, and specifying material directions for the yarns; the material parameters include yarn parameters and matrix parameters; Setting periodic boundary conditions through the ABD matrix plug-in to apply six periodic boundary conditions to the mesoscale RVE, thereby realizing six variations of the mesoscale RVE; Analyzing the deformation data of the six variations of the mesoscale RVE using the ABAQUS finite element analysis software, thereby calculating the ABD matrix of the thin woven composite material; The six periodic boundary conditions are applied to the mesoscale RVE to achieve the six variations of the mesoscale RVE, including: Establishing a plurality of reference points RP on the surface surrounding the mesoscale RVE, wherein the plurality of reference points RP form the mid-surface edge of the mesoscale RVE; The unit nodes of the upper surface edge in the periodic grid with the reference point RP as the base point are defined as upper nodes, and the unit nodes of the lower surface edge in the periodic grid with the reference point RP as the base point are defined as lower nodes; the reference point RP, the upper node, and the lower node in the same periodic grid are rigidly linked using MPC beam units to achieve coupling of all degrees of freedom of the node and RP in all directions; Applying six periodic boundary conditions to the plurality of reference points RP to perform deformation loading on the mid-surface of the mesoscale RVE, thereby achieving deformation control of the entire mesoscale RVE; Applying six periodic boundary conditions to the mesoscale RVE to achieve six variations of the mesoscale RVE also includes: if the upper node / the lower node is offset from the surrounding surfaces of the mesoscale RVE, then: The coordinates of the two vertices of the mid-body diagonal of the mesoscale RVE are identified by the ABAQUS finite element analysis software. and ; A rectangular parallelepiped region is established by the ABAQUS finite element analysis software. The coordinates of the two vertices of the diagonal line of the rectangular parallelepiped region are: and , determining the boundary of the rectangular parallelepiped region according to the coordinates of two vertices of the mid-body diagonal line of the rectangular parallelepiped region; all grid nodes in the rectangular parallelepiped region form a first set; ; Where n is the tolerance value; All mesh nodes in the first set are sorted according to the size of their x-coordinates and y-coordinates. The upper and lower endpoints are found among the nodes with the same x-coordinates and y-coordinates. The midpoint of the upper and lower endpoints is set as the reference point RP. The nodes with the same x-coordinates and y-coordinates are rigidly linked to the reference point RP using MPC beam units.

2. The mesoscale calculation method for predicting the ABD matrix of thin woven composite materials according to claim 1 is characterized in that: The method of generating a mesoscale RVE for a woven composite material using composite material modeling software includes: Obtaining the gap size S between adjacent warp yarns and adjacent weft yarns, the length L, width W, and height H of the mesoscale RVE, thereby obtaining a weaving form of the woven composite material; selecting a periodic mesoscale RVE according to the weaving form; A cross-sectional micrograph of the woven composite material is taken, and the cross-sectional shape, size, and yarn undulation path of each yarn are determined based on the cross-sectional micrograph, thereby obtaining the true geometric characteristics of the yarn; the mesoscale RVE is adjusted according to the true geometric characteristics of the yarn to obtain a mesoscale RVE model close to the actual situation.

3. The mesoscale calculation method for predicting the ABD matrix of thin woven composite materials according to claim 1 is characterized in that: The ABD matrix plug-in running in the ABAQUS finite element analysis software also includes: The ABD matrix plug-in deletes the default settings of the mesoscale RVE, which include material parameters of the yarn and matrix, solver settings, constraint equations, and boundary conditions.

4. The mesoscale calculation method for predicting the ABD matrix of thin woven composite materials according to claim 1 is characterized in that: The step of running an ABD matrix plug-in in the ABAQUS finite element analysis software, inputting material parameters in the ABD matrix plug-in, and specifying material directions for the yarns further comprises: A finite element analysis solver corresponding to the mesh division form of the mesoscale RVE is selected in the ABD matrix plug-in; the mesh division form includes volume mesh, voxel mesh and dry mesh.

5. The mesoscale calculation method for predicting the ABD matrix of thin woven composite materials according to claim 4 is characterized in that: When in the dry mesh form, the adhesive behavior between the yarn surfaces is approximated by setting the yarn surface adhesive contact properties, including: The relationship between stress and separation displacement between two adhesive surfaces is defined as: ； Where t and δ represent stress and separation displacement; the subscript n represents the normal direction, and s and t represent the two transverse shear directions; represents the normal stiffness, and represents the tangential stiffness, 、 and Collectively referred to as uncoupled traction stiffness; Setting a stiffness value by means of the uncoupled traction stiffness; The stiffness values ​​between the nodes in contact with each other on the surfaces of two yarns are used to generate an association relationship between the nodes in contact with each other, thereby approximating the adhesive behavior between the yarn surfaces.

6. The mesoscale calculation method for predicting the ABD matrix of thin woven composite materials according to claim 1 is characterized in that: The six periodic boundary conditions are applied to the mesoscale RVE to achieve the six variations of the mesoscale RVE, including: Apply displacement boundary conditions as a percentage of unit strain.

7. The mesoscale calculation method for predicting the ABD matrix of thin woven composite materials according to claim 1 is characterized in that: The composite material modeling software includes Texge and Digmat.

8. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the steps of the mesoscale calculation method for predicting the ABD matrix of a thin woven composite material described in any one of claims 1 to 7 are implemented.

Images