Establishment method of microcosmic finite element model of reinforced fiber concrete structure

Abstract

The invention discloses a method for establishing a microcosmic finite element model of a reinforced fiber concrete structure, and belongs to the technical field of concrete microcosmic model.The method comprises the steps that an MATLAB language is combined with ANSYS / LS-DYNA finite element software for analog simulation to obtain a four-phase finite element model of the reinforced concrete structure, and according to the method, reinforcing steel bars meeting requirements are arranged in a preset area; aggregate and fiber are put in the putting area, it is ensured that the aggregate, the fiber and the steel bars do not invade one another, modeling is conducted on a large number of fibers, accurate modeling is conducted on the directionally-arranged steel bars, the aggregate with different particle sizes is put into the modeling system, and the defect that the aggregate, the steel bars and the aggregate invade one another is overcome; the simulation model is very close to a real concrete structure, the model generation time is short, and the model precision is high.

Description

Technical Field

[0001] The invention belongs to the technical field of concrete micro-modeling, and in particular relates to a method for establishing a micro-finite element model of a reinforced fiber concrete structure. Background Art

[0002] Concrete is usually simplified as a homogeneous material at the macroscopic scale for mechanical analysis, but in fact it is a typical heterogeneous material. At present, many concrete models rely on macroscopic mechanical methods and ignore the complex structural characteristics inside concrete, which is suitable for the overall study of material properties to a certain extent. However, this simplified method cannot fully reveal the complexity, deformation mode and failure mechanism inside concrete. Experimental research has certain limitations in observing the internal deformation and failure of concrete. With the rapid development of computer technology, finite element simulation has become a widely used tool. By establishing a mesoscale model of concrete, the response of the material under dynamic and static loads can be deeply studied. This method can more comprehensively understand the actual performance of concrete and provide theoretical support for engineering design and building structure optimization. In order to make the mesoscale numerical simulation closer to the actual situation, it is crucial to establish an accurate mesoscale model. In this process, the numerical characterization of aggregates and the random distribution of fibers are important links in the mesoscale simulation of concrete. Therefore, it is crucial to develop a mesoscale model that can accurately reflect the internal components and their distribution of concrete.

[0003] At present, there are two main methods for modeling concrete mesoscopic models. One is based on the real concrete CT scanning results, and the corresponding concrete mesoscopic model is generated on the basis of the CT image, such as the CT image recognition method adopted in the invention patent applications (CN 105699258 A, CN105427382A); the other is through parametric modeling, such as the invention patent applications (CN118133493A, CN113987750A) using programming language to generate three-dimensional spatial information of aggregate, mortar and interface transition zone in concrete, and import it into finite element simulation software to establish the mesoscopic model.

[0004] Most of the current concrete meso models have the problem of single components and fail to accurately characterize the various components in concrete. These models usually face several key problems: first, the grading control of aggregates is insufficient, resulting in inaccurate distribution of aggregates in the model; second, the model still needs to be optimized when dealing with spatial randomness. In addition, when considering the random distribution of internal components, the existing concrete meso models often do not effectively avoid spatial overlap or intersection between components. These problems limit the accuracy and reliability of the model and affect the comprehensive understanding of the performance of concrete materials. Therefore, improving the accuracy of meso models in component characterization, grading control, spatial randomness processing, and interaction between components is crucial to further promote concrete meso research and application. Summary of the invention

[0005] The purpose of the present invention is to provide a method for establishing a mesoscopic finite element model of a reinforced fiber concrete structure which is easy to operate and has high accuracy, so as to solve the problems reflected in the background technology.

[0006] In order to achieve the above object, the present invention adopts the following technical solutions:

[0007] The method for establishing a microscopic finite element model of a reinforced fiber concrete structure is characterized in that the finite element model of the reinforced concrete structure is obtained by simulating using MATLAB language in conjunction with ANSYS / LS-DYNA finite element software, and the method is characterized in that the method comprises the following steps:

[0008] S1. Determine the size and volume of the concrete cube: the length of the cube is LX, the width of the cube is LY, and the height of the cube is LZ, then the volume of the concrete is: V = LX·LY·LZ;

[0009] S2. Complete the setting and calculation of steel bar related parameters:

[0010] According to the thickness of the concrete cover in the X-axis, Y-axis and Z-axis directions, and the length, width and height of the cube, calculate the length of the steel bars in the X-axis and Y-axis directions; according to the number of transverse and longitudinal steel bars, calculate the spacing of the transverse and longitudinal steel bars; arrange the first transverse and longitudinal steel bars in the delivery area and determine the coordinates, complete the steel bar arrangement according to the spacing, and save the steel bar coordinates;

[0011] S3. Complete the setting of aggregate-related parameters and aggregate placement:

[0012] The concrete placement area is evenly divided into multiple unit grids, and the placement space is noded according to the unit grid size of the finite element model. The model uses eight-node hexahedral units. Aggregates are placed in a large-to-small order to simulate the aggregate particle gradation in real concrete and ensure that the aggregate sphere does not exceed the boundary and that there is no intrusion and overlap between aggregates.

[0013] S4. Determine and adjust the intrusion between aggregate and steel bar:

[0014] After the pre-placement of aggregates is completed, the intrusion between aggregates and steel bars is judged, and the intruded aggregates are processed. The processed aggregates need to be re-judged whether they have intrusion with other aggregates, until all aggregates and all steel bars do not intrude and intersect, and the aggregate coordinates are saved. The processing methods include translation and diameter change, etc.

[0015] S5. Complete fiber-related parameter settings and fiber placement:

[0016] Determine the minimum and maximum delivery coordinate points of the X, Y and Z axes, and then determine the fiber delivery range. Determine the delivery quantity according to the selected fiber size parameters, and perform random delivery. After delivery, adjust the fiber coordinates to ensure that they are within the boundaries.

[0017] S6. Determine and adjust the intrusion between fiber and aggregate:

[0018] After all the fibers are placed and processed, the intrusion between the fibers and the aggregate is determined. When intrusion occurs between the fibers and the aggregate, the intruded fibers are handled by moving them out of the aggregate and ensuring that the moved fibers are located inside the boundary.

[0019] S7. Determine and adjust the penetration between fiber and steel bar:

[0020] After completing the intrusion judgment and adjustment between all fibers and aggregates, the intrusion judgment between fibers and steel bars is performed, and the fibers intruding the steel bars are deleted to ensure that all fibers and all steel bars do not intrude and intersect, and the fiber coordinates are saved;

[0021] S8. Construct a microscopic finite element model of reinforced concrete structure:

[0022] After all the required fiber, aggregate and steel bar arrangements are completed, the coordinate information of the aggregate is converted into finite element mesh nodes, that is, each cell in the mesh is traversed through a triple loop, and the difference between the distance from the sphere center coordinate to the cell center point and the aggregate radius R is calculated. When the value is less than 0, the cell is divided into this aggregate, and it is marked and renamed by element update;

[0023] Convert the fiber and steel bar endpoint coordinates into finite element simulation nodes, renumber and arrange the node information to make it unit information recognizable by the finite element software, and define the steel bar unit information as part-1, the fiber unit information as part-2, the aggregate unit information as part-3, and the matrix unit information as part-4.

[0024] The fopen function in MATLAB software is used to save it in the form of a keyword K file that can be recognized by the finite element analysis software ANSYS / LS-DYNA, thus completing the establishment of the microscopic finite element model of the reinforced fiber concrete structure.

[0025] As a further description of the technical solution of the present invention:

[0026] After all the steel bars are arranged in S2, the starting coordinates of the horizontal and vertical steel bars are saved in g_node1 in matrix form, and the ending coordinates of the horizontal and vertical steel bars are saved in g_node2 in matrix form;

[0027] In S7, the fiber starting point coordinates are stored in end_node1 in matrix form, and the fiber ending point coordinates are stored in end_node2 in matrix form.

[0028] As a further description of the technical solution of the present invention:

[0029] When the aggregate coordinate information in S8 is converted into finite element simulation model information,

[0030] Extract the spatial coordinate information of the first aggregate, calculate the distance R-1 from the center point of each unit to the center of the aggregate sphere, and determine the unit type division based on the relationship between R-1 and the radius R of the aggregate; when R-1 is smaller than the radius R, calculate the index n1 of this unit in the matrix Element, and divide this unit into aggregate; after completion, enter the unit replacement of the second aggregate until the unit replacement of all aggregates is completed; divide the units except the aggregate into the matrix, and save all the updated unit node information in the Element in the form of a matrix;

[0031] When arranging node information: create a column vector BH4 with a length equal to the number of steel bars O, define its value as 0, connect BH4 with g_node1 and g_node2 to form a longer node vector; extract the row number A of end_node1 and the row number B of end_node2, create a column vector BH from 1 to A+B+O for subsequent node numbering; arrange the node information in end_node1 and end_node2 in order, and connect the arranged node matrix with BH to form a longer node vector; initialize a zero matrix NODE of size (3·A+3·O,5) to store the coordinates and other information of the node, and arrange the processed node information in the required order to obtain the total node set;

[0032] When arranging unit information: create a column vector BH1 from 1 to A+O, create a column vector BH2 ​​from A+O+1 to A+B+2·O, create a column vector BH3 from A+B+2·O+1 to A+B+3·O for element numbering, create a column vector PART with a length of A+B, define its row values ​​from 1 to A as 1, for steel finite element component naming, define its A+1 to A+O values ​​as 2, for fiber finite element component naming. Combine PART, BH1, BH2, and BH3 into an overall unit matrix.

[0033] As a further description of the technical solution of the present invention:

[0034] The total number of units Nmesh in S3 is calculated as follows: Nmesh = Nx·Ny·Nz,

[0035] The total number of nodes Nnoda is calculated as follows: Nnoda = (Nx+1)·(Ny+1)·(Nz+1),

[0036] Where Nx is the number of units in the X-axis direction of the delivery range, Ny is the number of units in the Y-axis direction of the delivery range, and Nz is the number of units in the Z-axis direction of the delivery range;

[0037] In S3, the concrete placement area is evenly divided according to the unit size in each direction. i, j, and k represent the indexes in the x, y, and z directions, respectively, where i∈[1, Nx], j∈[1, Ny], and k∈[1, Nz], and n represents the global number of the node. The matrix Node0 is used to record the number of each node and its position in the x, y, and z directions. The calculation formula for n is as follows:

[0038] n=(i+1)+j·(Nx+1)+k·(Nx+1)·(Ny+1),

[0039] Perform a linear transformation on the matrix Node0 to adjust the node coordinates:

[0040]

[0041] As a further description of the technical solution of the present invention:

[0042] After completing the adjustment of the node coordinates, the corresponding unit information is generated. The model uses an eight-node hexahedral unit, each unit corresponds to eight nodes, i, j, k represent the indexes in the x, y, and z directions, respectively, where i∈[1, Nx-1], j∈[1, Ny-1], k∈[1, Nz-1], n1 represents the global number of the unit, and the matrix Element is used to save the unit node information, where the calculation formula of n1 is as follows: n1=(i+1)+j·Nx+k·Nx·Ny,

[0043] Element node 1 calculation formula: Element1 = (i+1) + j·(Nx+1) + k·(Nx+1)·(Ny+1),

[0044] Element node 2 calculation formula: Element2 = (i+2) + j·(Nx+1) + k·(Nx+1)·(Ny+1),

[0045] Element node 3 calculation formula: Element3 = (i+2)+(j+1)·(Nx+1)+k·(Nx+1)·(Ny+1),

[0046] Element node 4 calculation formula: Element4 = (i+1)+(j+1)·(Nx+1)+k·(Nx+1)·(Ny+1),

[0047] Element node 5 calculation formula: Element5 = (i+1) + j·(Nx+1) + (k+1)·(Nx+1)·(Ny+1),

[0048] Element node 6 calculation formula: Element6 = (i+2) + j·(Nx+1) + (k+1)·(Nx+1)·(Ny+1),

[0049] Element node 7 calculation formula: Element7 = (i+2) + (j+1)·(Nx+1) + (k+1)·(Nx+1)·(Ny+1),

[0050] Element node 8 calculation formula: Element8 = (i+1) + (j+1)·(Nx+1) + (k+1)·(Nx+1)·(Ny+1),

[0051] After completing the calculation of the unit nodes, the unit information is combined into the overall unit matrix Element, and the unit's global index number n1 is added to it to complete the unit modeling.

[0052] As a further description of the technical solution of the present invention:

[0053] When placing aggregate in S3, the aggregate radius is calculated first. The calculation steps are as follows:

[0054] The initial value of Asize is 0, and the rand function rand(1) is used, which represents generating a random number between 0 and 1.

[0055] When Asize<(Vpersen·3 / 4), the aggregate radius calculation formula is as follows:

[0056] R=2 / 3·Maxsize-a+rand(1)·1 / 3·Maxsize-a,

[0057] When Vpersen·3 / 4 < Asize ≤ Vpersen·11 / 12, the calculation formula for the aggregate radius is as follows:

[0058] R = 1 / 2·Maxsize - a + rand(1)·1 / 6·Maxsize - a,

[0059] When Vpersen·11 / 12 < Asize, the calculation formula for the aggregate radius is as follows:

[0060] R = Minsize - a + rand(1)·Minsize - a,

[0061] After calculating the aggregate radius, the coordinates of the aggregate center (Cx, Cy, Cz) are placed. These coordinates are randomly generated within the restricted range to ensure that the aggregate spheres do not exceed the boundary, and the boundary is defined by Lx, Ly, Lz:

[0062] Cx = Minsize - a + rand(1)·(Lx - 2·Minsize - a),

[0063] Cy = Minsize - a + rand(1)·(Ly - 2·Minsize - a),

[0064] Cz = Minsize - a + rand(1)·(Lz - 2·Minsize - a),

[0065] After placing the coordinates of the aggregate center, it is necessary to judge whether there is intrusion and overlap between the aggregates. For each randomly generated aggregate, it is necessary to judge whether there is intrusion and overlap with all the previously generated aggregates. By the corresponding relationship between the distance of the center coordinates and the sum of the radii, it is judged whether there is intrusion and overlap between the aggregates. If there is overlap, it is necessary to regenerate and re-judge;

[0066] After placing the coordinates of the aggregate center, it is also necessary to perform boundary processing on the aggregates to ensure that all aggregates meet the boundary conditions. Take A1 as the minimum value among the coordinates of the aggregate center (Cx, Cy, Cz), and the judgment conditions are as follows:

[0067] R < A1,

[0068] (Cx + R) ≤ Lx,

[0069] (Cy + R) ≤ Ly,

[0070] (Cz + R) ≤ Lz,

[0071] The volume formula of the aggregate is:

[0072] Vtotal = 4 / 3·π·R 3 。

[0073] As a further description of the technical solution of the present invention:

[0074] In S4, when determining the intrusion of aggregate and steel bars, the radius R of the spherical aggregate and its corresponding spherical center coordinate PointP (Cx, Cy, Cz) are extracted, and the starting point coordinate Point-A (X G 1. Y G 1. Z G 1) and the end point coordinates Point-B (X G 2. Y G 2. Z G 2) Calculate the shortest distance from a point on the steel bar to the center of the aggregate sphere using the following formula:

[0075]

[0076]

[0077]

[0078]

[0079] Dist=|PointP-Point C|,

[0080] Finally, Dist is taken as the shortest distance from the aggregate to the steel bar. By judging the relationship between Dist and the cross-section of the aggregate steel bar, it is decided how to translate the aggregate. The Dist judgment conditions are as follows:

[0081] When Dist<(R+g-dia), it is considered that the aggregate and the steel bar have invaded, and a secondary judgment is made.

[0082] When Dist<(R+g-dia)·3 / 5, the spatial position of the aggregate is translated, and the radius is regenerated during the translation process to ensure that the translated aggregate does not invade the steel bars or other aggregates, forming NewR and NewPointA, where:

[0083] NewR=(R+g-dia) / 2,

[0084]

[0085]

[0086] NewPointA=PointA+direction·0.2·(g-dia+R),

[0087] When (g-dia+0.01)·3 / 5≤Dist<(R+0.01), a new radius is generated:

[0088] NewR = 0.64·(g-dia+R).

[0089] As a further description of the technical solution of the present invention:

[0090] The calculation formula for the fiber delivery range in S5 is as follows:

[0091] x-min=x0,

[0092] x-max=x0+Lx1,

[0093] y-min=y0,

[0094] y-max=y0+Ly1,

[0095] z-min=z0,

[0096] z-max=z0+Lz1,

[0097] Among them, x0 is the minimum delivery coordinate point of the delivery range on the X axis, y0 is the minimum delivery coordinate point of the delivery range on the Y axis, z0 is the minimum delivery coordinate point of the delivery range on the Z axis, Lx1 is the maximum delivery coordinate point of the delivery range on the X axis, Ly1 is the maximum delivery coordinate point of the delivery range on the Y axis, and Lz1 is the maximum delivery coordinate point of the delivery range on the Z axis;

[0098] The calculation formula for the fiber delivery quantity m is as follows:

[0099] As a further description of the technical solution of the present invention:

[0100] The method for determining the starting point and direction vector of the fiber in S5 is as follows:

[0101] Determine the starting point and direction vector of the fiber, and use the rand function rand(1), which represents the generation of a random number between 0 and 1. The formula for calculating the starting point coordinates is as follows:

[0102] Randomly generate the X coordinate of the starting point of each fiber: x1 = LX rand (1),

[0103] The Y coordinate of the starting point of each fiber is randomly generated: y1 = LY rand (1),

[0104] Randomly generate the Z coordinate of the starting point of each fiber: z1 = LZ rand (1),

[0105] Use angle1 to represent the angle of the direction vector in the xy plane, and angle2 to represent the inclination angle of the direction vector in the zy plane. The calculation formula is as follows

[0106] angle1=2·π·rand(1),

[0107] angle2=2·π·rand(1),

[0108] angle1 and angle2 represent random angles in the range of 0 to 2π;

[0109] The fiber end point coordinates are calculated using the starting point coordinates and random angles angle1 and angle2. The calculation formula is as follows:

[0110] The X coordinate of the end point of each fiber is randomly generated as: x2 = x1 + Ln·cos(angle2)·cos(angle1),

[0111] The Y coordinate of the end point of each fiber is randomly generated as: y2 = y1 + Ln·cos(angle2)·sin(angle1),

[0112] The Z coordinate of the end point of each fiber is randomly generated as: z2=z1+Ln·sin(angle2),

[0113] When adjusting the fiber coordinates, when the starting point coordinates and the end point coordinates of the fiber exceed the minimum delivery boundary of the delivery area, the coordinates are automatically corrected to be offset 1 mm inward from the minimum delivery boundary; when the starting point coordinates and the end point coordinates exceed the maximum delivery boundary of the delivery area, the coordinates are automatically corrected to be offset 1 mm inward from the maximum delivery boundary.

[0114] As a further description of the technical solution of the present invention:

[0115] When judging the intrusion between fiber and aggregate in S6, if intrusion is considered, the intruded fiber will be processed by translating it out of the aggregate. After all fibers are placed and invaded, the fiber placement range is judged to ensure that all fibers are within the placement range. When the starting point coordinates and the end point coordinates exceed the minimum placement boundary, the coordinates are automatically corrected to the minimum placement boundary with an inward offset of 1 mm. When the starting point coordinates and the end point coordinates exceed the maximum placement boundary, the coordinates are automatically corrected to the maximum placement boundary with an inward offset of 1 mm.

[0116] In summary, due to the adoption of the above technical solution, the beneficial effects of the present invention are:

[0117] Currently, most concrete mesoscopic modeling technologies only include a single random fiber component or random aggregate component, and rarely describe all components of the concrete structure. This method models a large number of fibers, accurately models the directionally arranged steel bars, and also incorporates aggregates of different particle sizes, solving the mutual intrusion defects among fibers, steel bars, and aggregates. The simulation model is very close to the actual concrete structure, with a short model generation time and high model accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0118] Figure 1 It is a finite element model diagram of the steel bar of the reinforced fiber concrete structure of the present invention;

[0119] Figure 2 This is a finite element model diagram of the reinforced fiber concrete structure aggregate of the present invention;

[0120] Figure 3 This is a finite element model diagram of the reinforced fiber concrete structure of the present invention;

[0121] Figure 4 This is the finite element model diagram of the reinforced fiber concrete structure section of the present invention:

[0122] Figure 5 This is the spatial arrangement diagram of the steel bar aggregate of the reinforced fiber concrete structure of the present invention:

[0123] Figure 6 It is a modeling flow chart of the reinforced fiber concrete structure of the present invention. DETAILED DESCRIPTION

[0124] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0125] The present invention obtains a four-phase finite element model of a concrete structure by performing simulation using MATLAB language in conjunction with ANSYS / LS-DYNA finite element software.

[0126] The following takes a single-layer reinforced fiber concrete structure as an example to construct a model and describe the implementation mode of the present invention in detail. The established model is shown in FIG. Figure 1-6 The detailed steps are as follows:

[0127] Step 1: Generate basic information of ultra-high performance concrete reinforced composite structure in MATLAB software

[0128] Generate two rows of horizontal and vertical line segments in the pointed area. Since the geometric shape of the steel bar is similar to the line segment, the line segment arranged in a directional manner can be considered as the steel bar. On this basis, generate random fibers. Since the geometric shape of the fiber is similar to the line segment, the random line segment can be considered as the fiber. The specific process is as follows:

[0129] Step 1.1: Determine the concrete cube dimensions and volume

[0130] Take a cube as an example to determine the size of the cube concrete. The length of the cube is LX, the width of the cube is LY, and the height of the cube is LZ. Then the volume of the concrete is: V = LX · LY · LZ.

[0131] The relevant parameters of the steel bar are set as follows: the concrete cover thickness px in the X-axis direction, the concrete cover thickness py in the Y-axis direction, the concrete cover thickness pz in the Z-axis direction, and the steel bar section diameter G-dia. The steel bar length is:

[0132] G-lx=LX-2 px,

[0133] G-ly=LY-2·py,

[0134] Among them, G-lx is the length of the steel bar in the X-axis direction, and G-ly is the length of the steel bar in the Y-axis direction;

[0135] Determine the number of longitudinal reinforcement arrangements nx and the number of transverse reinforcement arrangements ny, and calculate the longitudinal reinforcement spacing JX and the transverse reinforcement spacing JY:

[0136] JX=G-lx / (nx-1),

[0137] JY = G-ly / (ny-1).

[0138] Step 1.2: Complete the setting and calculation of steel bar related parameters

[0139] First, specify the starting coordinates of the first transverse reinforcement in the placement area, calculate the end coordinates of the first reinforcement using reinforcement-related parameters, save the three coordinates of the reinforcement in the Cartesian coordinate system, translate the starting coordinates and end coordinates of the first transverse reinforcement, and obtain the starting coordinates and end coordinates of the second transverse reinforcement until the arrangement of the transverse reinforcement is completed.

[0140] The calculation formula for the starting point coordinates of the transverse reinforcement is as follows:

[0141] X1 g =px,

[0142] Y1 g =py+JY·(i-1)i∈[1,ny],

[0143] Z1 g =pz,

[0144] The calculation formula for the coordinates of the end point of the transverse reinforcement is as follows:

[0145] X2 g =px+G-lx,

[0146] Y2 g =py+JY·(i-1)i∈[1,ny],

[0147] Z2 g =pz,

[0148] After completing the arrangement of all transverse steel bars, first calculate the starting coordinates of the first longitudinal steel bar in the placement area, calculate the end coordinates of the first longitudinal steel bar using the steel bar related parameters, save the three coordinates of the steel bar in the Cartesian coordinate system in matrix form, translate the starting coordinates and end coordinates of the first longitudinal steel bar, and obtain the starting coordinates and end coordinates of the second longitudinal steel bar until the arrangement of the longitudinal steel bars is completed.

[0149] The calculation formula for the starting point coordinates of the longitudinal reinforcement is as follows:

[0150] X3 g =px+JX·(i-1)i∈[1,nx],

[0151] Y3 g =py,

[0152] Z3 g =pz+G-dia,

[0153] The calculation formula for the end point coordinates of the longitudinal reinforcement is as follows:

[0154] X4 g =px+JX·(i-1)i∈[1,nx],

[0155] Y4 g =py+G-ly,

[0156] Z4 g =pz+G-dia,

[0157] In the above coordinate naming, for easy distinction, the number 1 in the coordinate symbol represents the starting point of the transverse reinforcement, the number 2 represents the end point of the transverse reinforcement, the number 3 represents the starting point of the longitudinal reinforcement, and the number 4 represents the end point of the longitudinal reinforcement.

[0158] Step 1.3: Save the rebar coordinates

[0159] After completing the arrangement of all steel bars, the starting coordinates of the horizontal and vertical steel bars are saved in g_node1 in matrix form, and the ending coordinates of the horizontal and vertical steel bars are saved in g_node2 in matrix form.

[0160] Step 2: Aggregate placement area and related parameter settings

[0161] Step 2.1: Aggregate placement area setting

[0162] Nx is the number of units in the X-axis direction of the delivery range, Ny is the number of units in the Y-axis direction of the delivery range, Nz is the number of units in the Z-axis direction of the delivery range, Msx is the unit size in the X-axis direction of the delivery range, Msy is the unit size in the Y-axis direction of the delivery range, Msz is the unit size in the Z-axis direction of the delivery range, Maxsize-a is the maximum radius of the aggregate, Minsize-a is the minimum radius of the aggregate, Vpersen is the volume content of the aggregate, and the delivery range calculation formula is as follows:

[0163] Lx=Nx·Msx,

[0164] Ly=Ny·Msy,

[0165] Lz=Nz·Msz,

[0166] The total number of units is calculated as follows: Nmesh = Nx·Ny·Nz,

[0167] The total number of nodes is calculated as follows: Nnoda = (Nx+1)·(Ny+1)·(Nz+1),

[0168] The concrete placement area is evenly divided according to the unit size in each direction. i, j, k represent the indexes in the x, y, and z directions respectively, where i∈[1, Nx], j∈[1, Ny], k∈[1, Nz], and n represents the global number of the node. The matrix Node0 is used to record the number of each node and its position in the x, y, and z directions. The calculation formula for n is as follows:

[0169] n=(i+1)+j·(Nx+1)+k·(Nx+1)·(Ny+1),

[0170] Perform a linear transformation on the matrix Node0 to adjust the node coordinates:

[0171]

[0172] After the node coordinates are adjusted, the corresponding unit information is generated. The model uses eight-node hexahedral units, each unit corresponds to eight nodes, i, j, k represent the indexes in the x, y, and z directions, respectively, where i∈[1, Nx-1], j∈[1, Ny-1], k∈[1, Nz-1], n1 represents the global number of the unit, and the matrix Element is used to save the unit node information. The calculation formula for n1 is as follows: n1=(i+1)+j·Nx+k·Nx·Ny,

[0173] Element node 1 calculation formula: Element1 = (i+1) + j·(Nx+1) + k·(Nx+1)·(Ny+1),

[0174] Element node 2 calculation formula: Element2 = (i+2) + j·(Nx+1) + k·(Nx+1)·(Ny+1),

[0175] Element node 3 calculation formula: Element3 = (i+2)+(j+1)·(Nx+1)+k·(Nx+1)·(Ny+1),

[0176] Element node 4 calculation formula: Element4 = (i+1)+(j+1)·(Nx+1)+k·(Nx+1)·(Ny+1),

[0177] Element node 5 calculation formula: Element5 = (i+1) + j·(Nx+1) + (k+1)·(Nx+1)·(Ny+1),

[0178] Element node 6 calculation formula: Element6 = (i+2) + j·(Nx+1) + (k+1)·(Nx+1)·(Ny+1),

[0179] Element node 7 calculation formula: Element7 = (i+2) + (j+1)·(Nx+1) + (k+1)·(Nx+1)·(Ny+1),

[0180] Element node 8 calculation formula: Element8 = (i+1) + (j+1)·(Nx+1) + (k+1)·(Nx+1)·(Ny+1),

[0181] After completing the calculation of the unit nodes, the unit information is combined into the overall unit matrix Element, and the unit's global index number n1 is added to it to complete the unit modeling.

[0182] Step 2.2: Aggregate model placement

[0183] After completing the steel bar arrangement and unit grid division, the aggregate model is placed. The aggregate is placed in a large-to-small order. By controlling the different particle sizes of aggregates, the aggregate particle gradation in real concrete is simulated. The specific placement steps are as follows:

[0184] The initial value of Asize is 0, and the rand function rand(1) is used, which represents generating a random number between 0 and 1.

[0185] When Asize<(Vpersen·3 / 4), the aggregate radius calculation formula is as follows:

[0186] R = 2 / 3·Maxsize - a + rand(1)·1 / 3·Maxsize - a,

[0187] When Vpersen·3 / 4 < Asize ≤ Vpersen·11 / 12, the calculation formula for the aggregate radius is as follows:

[0188] R = 1 / 2·Maxsize - a + rand(1)·1 / 6·Maxsize - a,

[0189] When Vpersen·11 / 12 < Asize, the calculation formula for the aggregate radius is as follows:

[0190] R = Minsize - a + rand(1)·Minsize - a,

[0191] After calculating the aggregate radius, the coordinates of the aggregate center (Cx, Cy, Cz) are placed. These coordinates are randomly generated within the restricted range to ensure that the aggregate spheres do not exceed the boundary, where the boundary is defined by Lx, Ly, and Lz.

[0192] Cx = Minsize - a + rand(1)·(Lx - 2·Minsize - a),

[0193] Cy = Minsize - a + rand(1)·(Ly - 2·Minsize - a),

[0194] Cz = Minsize - a + rand(1)·(Lz - 2·Minsize - a),

[0195] After placing the coordinates of the aggregate center, it is necessary to judge whether there is intrusion and overlap between the aggregates. For each randomly generated aggregate, it is necessary to judge whether there is intrusion and overlap with all the previously generated aggregates. By the corresponding relationship between the distance of the center coordinates and the sum of the radii, it is judged whether the aggregates have intrusion and overlap. If there is overlap, it is necessary to regenerate and re - judge.

[0196] After placing the coordinates of the aggregate center, it is also necessary to perform boundary processing on the aggregates to ensure that all aggregates meet the boundary conditions. Take A1 as the minimum value among the coordinates of the aggregate center (Cx, Cy, Cz), and the judgment conditions are as follows:

[0197] R < A1,

[0198] (Cx + R) ≤ Lx,

[0199] (Cy + R) ≤ Ly,

[0200] (Cz + R) ≤ Lz,

[0201] The volume formula of the aggregate is:

[0202] Vtotal=4 / 3·π·R 3 ,

[0203] Step 2.3: Determine the penetration between aggregate and reinforcement

[0204] After all the aggregates are pre-placed, the intrusion between the aggregates and the steel bars is determined, the radius R of the spherical aggregate and its corresponding spherical center coordinates PointP (Cx, Cy, Cz) are extracted, and the starting point coordinates Point-A (X G 1. Y G 1. Z G 1) and the end point coordinates Point-B (X G 2. Y G 2. Z G 2) Calculate the shortest distance from a point on the steel bar to the center of the aggregate sphere. Point C is the coordinate of the point where the shortest distance from the steel bar to the aggregate is located. The calculation formula is as follows:

[0205]

[0206]

[0207]

[0208]

[0209] Dist=|PointP-Point C|,

[0210] Finally, Dist is taken as the shortest distance from the aggregate to the steel bar. By judging the relationship between Dist and the cross-section of the aggregate steel bar, it is decided how to translate the aggregate. The judgment conditions are as follows:

[0211] When Dist<(R+g-dia), it is considered that the aggregate and the steel bar have invaded, and a secondary judgment is made.

[0212] When Dist<(R+g-dia)·3 / 5, the spatial position of the aggregate is translated, and the radius is regenerated during the translation process to ensure that the translated aggregate does not invade the steel bars or other aggregates, forming NewR and NewPointA, where:

[0213] NewR=(R+g-dia) / 2,

[0214]

[0215]

[0216] NewPointA=PointA+direction·0.2·(g-dia+R),

[0217] When (g-dia+0.01)·3 / 5≤Dist<(R+0.01), a new radius is generated:

[0218] NewR=0.64·(g-dia+R),

[0219] Aggregates that have been processed by translation (including translation, diameter change and regeneration, etc.) need to be re-judged whether they have invaded other aggregates, until all aggregates and all steel bars do not invade or intersect, and the aggregate coordinates are saved.

[0220] Step 3: Fiber delivery area and related parameter settings

[0221] Step 3.1: Fiber calculation and placement

[0222] x0 is the minimum delivery coordinate point of the delivery range on the X axis, y0 is the minimum delivery coordinate point of the delivery range on the Y axis, z0 is the minimum delivery coordinate point of the delivery range on the Z axis, Lx1 is the maximum delivery coordinate point of the delivery range on the X axis, Ly1 is the maximum delivery coordinate point of the delivery range on the Y axis, Lz1 is the maximum delivery coordinate point of the delivery range on the Z axis, dia is the fiber cross-sectional diameter, ratio is the fiber volume content, and Ln is the fiber length. The delivery range calculation formula is as follows:

[0223] x-min=x0,

[0224] x-max=x0+Lx1,

[0225] y-min=y0,

[0226] y-max=y0+Ly1,

[0227] z-min=z0,

[0228] z-max=z0+Lz1,

[0229] The calculation formula for the fiber delivery quantity m is as follows:

[0230] Determine the starting point and direction vector of the fiber, and use the rand function rand(1), which represents the generation of a random number between 0 and 1. The formula for calculating the starting point coordinates is as follows:

[0231] Randomly generate the X coordinate of the starting point of each fiber: x1 = LX rand (1),

[0232] The Y coordinate of the starting point of each fiber is randomly generated: y1 = LY rand (1),

[0233] Randomly generate the Z coordinate of the starting point of each fiber: z1 = LZ rand (1),

[0234] Use angle1 to represent the angle of the direction vector in the xy plane, and angle2 to represent the inclination angle of the direction vector in the zy plane. The calculation formula is as follows

[0235] angle1=2·π·rand(1),

[0236] angle2=2·π·rand(1),

[0237] angle1 and angle2 represent random angles in the range of 0 to 2π;

[0238] The fiber end point coordinates are calculated using the starting point coordinates and random angles angle1 and angle2. The calculation formula is as follows:

[0239] The X coordinate of the end point of each fiber is randomly generated as: x2 = x1 + Ln·cos(angle2)·cos(angle1),

[0240] The Y coordinate of the end point of each fiber is randomly generated as: y2 = y1 + Ln·cos(angle2)·sin(angle1),

[0241] The Z coordinate of the end point of each fiber is randomly generated as: z2=z1+Ln·sin(angle2),

[0242] Step 3.2: Fiber coordinate adjustment

[0243] After all the fibers are placed, their starting point coordinates and end point coordinates are determined. When their starting point coordinates and end point coordinates exceed the minimum placement boundary of the placement area, their coordinates are automatically corrected to be 1 mm inwardly offset from the minimum placement boundary. The calculation formula is as follows:

[0244] x1=x-min+0.001,

[0245] y1=y-min+0.001,

[0246] z1=z-min+0.001,

[0247] x2=x-min+0.001,

[0248] y2=y-min+0.001,

[0249] z2=z-min+0.001,

[0250] When the starting point coordinates and the end point coordinates exceed the maximum delivery boundary of the delivery area, the coordinates are automatically corrected to be offset 1 mm inward from the maximum delivery boundary. The calculation formula is as follows:

[0251] x1=x-max-0.001,

[0252] y1=y-max-0.001,

[0253] z1=z-max-0.001,

[0254] x2=x-max-0.001,

[0255] y2=y-max-0.001,

[0256] z2=z-max-0.001,

[0257] Step 3.3: Fiber and aggregate intrusion determination

[0258] After all the fibers are placed and processed, the intrusion between the fibers and the aggregate is determined. The following is one of the specific methods, and of course other methods or other names can also be selected. The intrusion determination between the fibers and the steel bars is also a specific method, and other methods can also be selected.

[0259] Extract the center coordinates PointP (Cx, Cy, Cz) of the first spherical aggregate and the corresponding sphere radius R, and extract the starting coordinates C (X 1 1. Y 1 1. Z 1 1) and the end point coordinates D(X 1 2. Y 1 2. Z 1 2) By calculating the shortest distance from the center of the spherical aggregate to the fiber, the relationship between it and the radius R of the spherical aggregate is determined to determine whether the fiber and the sphere have invaded. The calculation formula is as follows:

[0260]

[0261]

[0262]

[0263]

[0264] Distance = |PointP-PointB|,

[0265] By judging the size relationship between Distance and the sphere radius R, it is determined whether the fiber and aggregate are intruded. When the fiber and aggregate are intruded, the intruded fiber is moved out of the aggregate to handle it. The specific calculation formula is as follows:

[0266]

[0267]

[0268]

[0269] PointG=PointH+0.15·R·T1,

[0270]

[0271]

[0272] NewX 2 1 = PointGx - 0.5 n T2x,

[0273] NewY 2 1 = PointGy - 0.5 n T2y,

[0274] NewZ 2 1 = PointGz - 0.5 n T2z,

[0275] NewX 2 2=PointGx+0.5·n·T2x,

[0276] NewY 2 2=PointGy+0.5·n·T2y,

[0277] NewZ 2 2=PointGz+0.5·n·T2z,

[0278] After all fibers are placed and invaded, the fiber placement range is determined to ensure that all fibers are within the placement range. When the starting point coordinates and the end point coordinates exceed the minimum placement boundary, the coordinates are automatically corrected to the minimum placement boundary offset by 1 mm. The calculation formula is as follows:

[0279] NewX 2 1 = x-min + 0.001,

[0280] NewY 2 1 = y-min + 0.001,

[0281] NewZ 2 1 = z-min + 0.001,

[0282] NewX 2 2 = x-min + 0.001,

[0283] NewY 2 2 = y-min + 0.001,

[0284] NewZ 2 2 = z-min + 0.001,

[0285] When the starting point coordinates and the end point coordinates exceed the maximum delivery boundary, the coordinates are automatically corrected to be 1 mm inwardly offset from the maximum delivery boundary. The calculation formula is as follows:

[0286] NewX 2 1 = x-max-0.001,

[0287] NewY 2 1 = y-max - 0.001,

[0288] NewZ 2 1 = z-max - 0.001,

[0289] NewX 2 2 = x-max-0.001,

[0290] NewY 2 2 = y-max - 0.001,

[0291] NewZ 2 2 = z-max - 0.001,

[0292] Step 3.3: Fiber and steel intrusion determination

[0293] After completing the intrusion judgment between all fibers and aggregates, the intrusion judgment between fibers and steel bars is performed, and the starting point coordinates of the first fiber are extracted: NewPOINT1 (NewX 2 1. NewY 2 1. NewZ 2 1),

[0294] and the end point coordinates: NewPOINT2(NewX 2 2. NewY 2 2. NewZ 2 2),

[0295] Extract the starting point coordinates of the first steel bar POINT-G1(X G 1,Y G 1,Z G 1) and the end point coordinates POINT-G1(X G 2,YG 2,Z G 2) When the fiber and the steel bar intersect in the matrix space, the shortest distance from the coordinates of the two end points of the fiber to the direction vector of the steel bar is calculated respectively. The calculation formula is as follows:

[0296]

[0297]

[0298]

[0299]

[0300]

[0301] Calculate the two intersection points of the vector with the shortest distance from the fiber end point to the steel bar and the steel bar direction vector. The calculation formula is as follows:

[0302]

[0303]

[0304] After finding the midpoint coordinates C of the two intersection points, use the midpoint coordinates C to calculate the shortest distance CD from the midpoint coordinates to the fiber. The calculation formula is as follows:

[0305]

[0306]

[0307]

[0308]

[0309]

[0310]

[0311]

[0312]

[0313] Finally, the minimum value among Dist, Dist1, and Dist2 is taken as the shortest distance from the fiber to the steel bar. By judging its relationship with the sum of the cross-sectional radius of the fiber and steel bar, it is decided whether to delete the fiber until all fibers and all steel bars do not intrude or intersect.

[0314] Step 3.4: Save fiber coordinates

[0315] After the fiber placement is completed, the starting point coordinates of the fiber that meets the requirements are saved in end_node1 in matrix form, and the ending point coordinates of the fiber are saved in end_node2 in matrix form.

[0316] Step 4: Construct a microscopic finite element model of the reinforced fiber reinforced concrete structure

[0317] After completing the arrangement of all fibers, aggregates, and steel bars that meet the requirements, the coordinate information of the aggregate is converted into finite element mesh nodes, that is, through a triple loop to traverse each cell in the grid, calculate the difference between the distance from the sphere center coordinate to the cell center point and the radius R. When its value is less than 0, its cell is divided into this aggregate, and it is marked and renamed. The fiber and steel bar end coordinates are converted into finite element simulation nodes, and the node information is renumbered and arranged to become unit information that can be recognized by the finite element software. After completing the arrangement of the steel bar unit information, it is defined as part-1, the fiber unit information is defined as part-2, the aggregate unit information is defined as part-3, and the matrix unit information is defined as part-4. Use the fopen function in the MATLAB software to save it in the form of a keyword K file that can be recognized by the finite element analysis software ANSYS / LS-DYNA, and the establishment of the microscopic finite element model of the reinforced fiber concrete structure is completed.

[0318] Step 4.1: Convert aggregate coordinate information into finite element simulation model information

[0319] Extract the spatial coordinate information of the first aggregate, calculate the distance R-1 from the center point of each unit to the center of the aggregate sphere, and determine the unit type division based on the relationship between R-1 and the radius R of the aggregate. When R-1 is less than the radius R, calculate the index n1 of this unit in the matrix Element, and divide this unit into aggregates. After completion, enter the unit replacement of the second aggregate until the unit replacement of all aggregates is completed. Divide the units except the aggregates into the matrix, and save all updated unit node information in the Element in the form of a matrix.

[0320] Step 4.2: Convert the steel bar coordinate and fiber coordinate information into finite element simulation model information

[0321] Node information arrangement: Create a column vector BH4 with a length equal to the number of steel bars O, define its value as 0, and connect BH4 with g_node1 and g_node2 to form a longer node vector. Extract the number of rows A of end_node1 (i.e. the number of nodes in end_node1) and the number of rows B of end_node2 (i.e. the number of nodes in end_node2), and create a column vector BH from 1 to A+B+O (the sum of the number of steel bars and the number of fibers) for subsequent node numbering. Arrange the node information in end_node1 and end_node2 in order, and connect the arranged node matrix with BH to form a longer node vector. Initialize a zero matrix NODE of size (3·A+3·O,5) to store the coordinates and other information of the nodes, and arrange the processed node information in the required order to obtain the total node set.

[0322] Unit information arrangement: create a column vector BH1 from 1 to A+O, create a column vector BH2 ​​from A+O+1 to A+B+2·O, create a column vector BH3 from A+B+2·O+1 to A+B+3·O for element numbering, create a column vector PART with a length of A+B, define its row values ​​from 1 to A as 1, used for steel finite element component naming, define its A+1 to A+O values ​​as 2, used for fiber finite element component naming. Combine PART, BH1, BH2, and BH3 into an overall unit matrix.

[0323] Step 4.3: Open the file model.k in write mode, write different model data into the file, add the file header information and format description, and obtain the finite element model file.

[0324] Step 4.4: Use the finite element analysis software ANSYS / LS-DYNA to directly open the obtained K file, in which the steel bar unit information is part-1, the fiber unit information is part-2, the aggregate unit information is part-3, and the matrix unit information is part-4.

[0325] The above description is only a preferred specific implementation manner of the present invention, but the protection scope of the present invention is not limited thereto. Any technician familiar with the technical field can make equivalent replacements or changes according to the technical scheme and inventive concept of the present invention within the technical scope disclosed by the present invention, which should be covered by the protection scope of the present invention.

Claims

1. The establishment method of the microscopic finite element model of reinforced fiber concrete structure is to obtain the finite element model of reinforced concrete structure by using MATLAB language combined with ANSYS / LS-DYNA finite element software for simulation, which is characterized by The steps include: S1. Determine the size and volume of concrete cube; S2. Complete the setting and calculation of steel bar related parameters: Arrange the first transverse and longitudinal reinforcement in the delivery area and determine the coordinates, complete the reinforcement arrangement according to the spacing, and save the reinforcement coordinates; S3. Complete the setting of aggregate-related parameters and aggregate placement: The concrete placement area is evenly divided into multiple unit grids, and the placement space is noded according to the unit grid size of the finite element model. The model uses eight-node hexahedral units. Aggregates are placed in a large-to-small order to simulate the aggregate particle gradation in real concrete and ensure that the aggregate sphere does not exceed the boundary and that there is no intrusion and overlap between aggregates. S4. Determine and adjust the intrusion between aggregate and steel bar: Determine the intrusion between aggregates and steel bars, process the intruded aggregates, and ultimately ensure that all aggregates and all steel bars do not intrude, and save the aggregate coordinates; S5. Complete fiber-related parameter settings and fiber placement: Determine the fiber placement range and quantity, and randomly place them. After placement, adjust the fiber coordinates to ensure that they are within the boundary; S6. Determine and adjust the intrusion between fiber and aggregate: When fiber and aggregate invade each other, the invaded fiber is handled by translating it out of the aggregate and ensuring that the translated fiber is located inside the boundary. S7. Determine and adjust the penetration between fiber and steel bar: Determine the intrusion between fibers and steel bars, delete the fibers that intrude into the steel bars, ensure that all fibers and all steel bars do not intrude or intersect, and save the fiber coordinates; S8. Construct a microscopic finite element model of reinforced concrete structure: After all the required fiber, aggregate and steel bar arrangements are completed, the coordinate information of the aggregate is converted into finite element mesh nodes, that is, each cell in the mesh is traversed through a triple loop, and the difference between the distance from the sphere center coordinate to the cell center point and the aggregate radius is calculated. When the value is less than 0, the cell is divided into this aggregate, and it is marked and renamed by element update; Convert the fiber and steel bar endpoint coordinates into finite element simulation nodes, renumber and arrange the node information to make it unit information recognizable by the finite element software, and define the steel bar unit information as part-1, the fiber unit information as part-2, the aggregate unit information as part-3, and the matrix unit information as part-4. The fopen function in MATLAB software is used to save it in the form of a keyword K file that can be recognized by the finite element analysis software ANSYS / LS-DYNA, thus completing the establishment of the microscopic finite element model of the reinforced fiber concrete structure.

2. The method for establishing a mesoscopic finite element model of a reinforced concrete structure according to claim 1, characterized in that: After all the steel bars are arranged in S2, the starting coordinates of the horizontal and vertical steel bars are saved in g_node1 in matrix form, and the ending coordinates of the horizontal and vertical steel bars are saved in g_node2 in matrix form; In S7, the fiber starting point coordinates are stored in end_node1 in matrix form, and the fiber ending point coordinates are stored in end_node2 in matrix form.

3. The method for establishing a mesoscopic finite element model of a reinforced concrete structure according to claim 2, characterized in that: When the aggregate coordinate information in S8 is converted into the finite element simulation model information, the spatial coordinate information of the first aggregate is extracted, and the distance R-1 from the center point of each unit to the center of the aggregate sphere is calculated. The unit type division is determined by the relationship between R-1 and the radius R of the aggregate; when R-1 is smaller than the radius R, the index n1 of this unit in the matrix Element is calculated, and this unit is divided into aggregates; after completion, the unit replacement of the second aggregate is entered until the unit replacement of all aggregates is completed; the units except the aggregates are divided into the matrix, and all updated unit node information is saved in the Element in the form of a matrix; When arranging node information: create a column vector BH4 with a length equal to the number of steel bars O, define its value as 0, connect BH4 with g_node1 and g_node2 to form a longer node vector; extract the row number A of end_node1 and the row number B of end_node2, create a column vector BH from 1 to A+B+O for subsequent node numbering; arrange the node information in end_node1 and end_node2 in order, and connect the arranged node matrix with BH to form a longer node vector; initialize a zero matrix NODE of size (3·A+3·O,5) to store the coordinates and other information of the node, and arrange the processed node information in the required order to obtain the total node set; When arranging unit information: create a column vector BH1 from 1 to A+O, create a column vector BH2 ​​from A+O+1 to A+B+2·O, create a column vector BH3 from A+B+2·O+1 to A+B+3·O for element numbering, create a column vector PART with a length of A+B, define its row values ​​from 1 to A as 1, which is used for naming steel finite element components, define its A+1 to A+O values ​​as 2, which is used for naming fiber finite element components, and combine PART, BH1, BH2, and BH3 into an overall unit matrix.

4. The method for establishing a mesoscopic finite element model of a reinforced concrete structure as claimed in claim 3, characterized in that: The total number of units Nmesh in S3 is calculated as follows: Nmesh = Nx·Ny·Nz, The total number of nodes Nnoda is calculated as follows: Nnoda = (Nx+1)·(Ny+1)·(Nz+1), Where Nx is the number of units in the X-axis direction of the delivery range, Ny is the number of units in the Y-axis direction of the delivery range, and Nz is the number of units in the Z-axis direction of the delivery range; In S3, the concrete placement area is evenly divided according to the unit size in each direction. i, j, and k represent the indexes in the x, y, and z directions, respectively, where i∈[1, Nx], j∈[1, Ny], and k∈[1, Nz], and n represents the global number of the node. The matrix Node0 is used to record the number of each node and its position in the x, y, and z directions. The calculation formula for n is as follows: n=(i+1)+j·(Nx+1)+k·(Nx+1)·(Ny+1), Perform a linear transformation on the matrix Node0 to adjust the node coordinates:

5. The method for establishing a mesoscopic finite element model of a reinforced concrete structure according to claim 4, characterized in that: After completing the adjustment of the node coordinates, the corresponding unit information is generated. The model uses an eight-node hexahedral unit, each unit corresponds to eight nodes, i, j, k represent the indexes in the x, y, and z directions, respectively, where i∈[1, Nx-1], j∈[1, Ny-1], k∈[1, Nz-1], n1 represents the global number of the unit, and the matrix Element is used to save the unit node information, where the calculation formula of n1 is as follows: n1=(i+1)+j·Nx+k·Nx·Ny, Element node 1 calculation formula: Element1 = (i+1) + j·(Nx+1) + k·(Nx+1)·(Ny+1), Element node 2 calculation formula: Element2 = (i+2) + j·(Nx+1) + k·(Nx+1)·(Ny+1), Element node 3 calculation formula: Element3 = (i+2)+(j+1)·(Nx+1)+k·(Nx+1)·(Ny+1), Element node 4 calculation formula: Element4 = (i+1)+(j+1)·(Nx+1)+k·(Nx+1)·(Ny+1), Element node 5 calculation formula: Element5 = (i+1) + j·(Nx+1) + (k+1)·(Nx+1)·(Ny+1), Element node 6 calculation formula: Element6 = (i+2) + j·(Nx+1) + (k+1)·(Nx+1)·(Ny+1), Element node 7 calculation formula: Element7 = (i+2) + (j+1)·(Nx+1) + (k+1)·(Nx+1)·(Ny+1), Element node 8 calculation formula: Element8 = (i+1) + (j+1)·(Nx+1) + (k+1)·(Nx+1)·(Ny+1), After completing the calculation of the unit nodes, the unit information is combined into the overall unit matrix Element, and the unit's global index number n1 is added to it to complete the unit modeling.

6. The method for establishing a mesoscopic finite element model of a reinforced ultra-high performance concrete structure according to claim 5, characterized in that: When placing aggregate in S3, the aggregate radius is calculated first. The calculation steps are as follows: The initial value of Asize is 0, and the rand function rand(1) is used, which represents generating a random number between 0 and 1. When Asize < (Vpersen·3 / 4), the calculation formula for the aggregate radius is as follows: R = 2 / 3·Maxsize - a + rand(1)·1 / 3·Maxsize - a, When Vpersen·3 / 4 < Asize ≤ Vpersen·11 / 12, the calculation formula for the aggregate radius is as follows: R = 1 / 2·Maxsize - a + rand(1)·1 / 6·Maxsize - a, When Vpersen·11 / 12 < Asize, the calculation formula for the aggregate radius is as follows: R = Minsize - a + rand(1)·Minsize - a, After completing the calculation of the aggregate radius, the coordinates of the aggregate sphere center (Cx, Cy, Cz) are placed. These coordinates are randomly generated within the restricted range to ensure that the aggregate spheres do not exceed the boundary, and the boundary is defined by Lx, Ly, Lz: Cx = Minsize - a + rand(1)·(Lx - 2·Minsize - a), Cy = Minsize - a + rand(1)·(Ly - 2·Minsize - a), Cz = Minsize - a + rand(1)·(Lz - 2·Minsize - a), After completing the placement of the aggregate sphere center coordinates, it is necessary to judge whether there is intrusion and overlap between the aggregates. For each randomly generated aggregate, it is necessary to judge the intrusion and overlap with all the previously generated aggregates. Through the corresponding relationship between the distance of the sphere center coordinates and the sum of the radii, it is judged whether the aggregates have intrusion and overlap. If there is overlap, it is necessary to regenerate and re-judge; After completing the placement of the aggregate sphere center coordinates, it is also necessary to perform boundary processing on the aggregates to ensure that all aggregates meet the boundary conditions. Take A1 as the minimum value among the aggregate sphere center coordinates (Cx, Cy, Cz), and the judgment conditions are as follows: R < A1, (Cx + R) ≤ Lx, (Cy + R) ≤ Ly, (Cz + R) ≤ Lz, The volume formula of the aggregate is: Vtotal=4 / 3·π·R 3 。 7. The method for establishing a mesoscopic finite element model of a reinforced ultra-high performance concrete structure according to claim 6, characterized in that: In S4, when determining the intrusion of aggregate and steel bars, the radius R of the spherical aggregate and its corresponding spherical center coordinate PointP (Cx, Cy, Cz) are extracted, and the starting point coordinate Point-A (X G 1. Y G 1. Z G 1) and the end point coordinates Point-B (X G 2. Y G 2. Z G 2) Calculate the shortest distance from a point on the steel bar to the center of the aggregate sphere using the following formula: Dist = |PointP - Point C|, Finally, take Dist as the shortest distance from the aggregate to the steel bar. By judging the size relationship between Dist and the cross-section of the aggregate steel bar, it is decided how to perform translation processing on the aggregate. The Dist judgment conditions are as follows: When Dist < (R + g - dia), it is considered that the aggregate and the steel bar have intrusion, and a secondary judgment is performed on it, When Dist < (R + g - dia)·3 / 5, the spatial position of the aggregate is translated, and the radius of the aggregate is regenerated when translating the aggregate, so as to ensure that the translated aggregate does not intrude into the steel bar and does not intrude into other aggregates, forming NewR and NewPointA, where: NewR = (R + g - dia) / 2, NewPointA = PointA + direction·0.2·(g - dia + R), When (g - dia + 0.01)·3 / 5 ≤ Dist < (R + 0.01), a new radius is regenerated: NewR = 0.64·(g-dia+R).

8. The method for establishing a mesoscopic finite element model of a reinforced ultra-high performance concrete structure according to claim 7, characterized in that: The calculation formula for the fiber delivery range in S5 is as follows: x-min=x0, x-max=x0+Lx1, y-min=y0, y-max=y0+Ly1, z-min=z0, z-max=z0+Lz1, Among them, x0 is the minimum delivery coordinate point of the delivery range on the X axis, y0 is the minimum delivery coordinate point of the delivery range on the Y axis, z0 is the minimum delivery coordinate point of the delivery range on the Z axis, Lx1 is the maximum delivery coordinate point of the delivery range on the X axis, Ly1 is the maximum delivery coordinate point of the delivery range on the Y axis, and Lz1 is the maximum delivery coordinate point of the delivery range on the Z axis; The calculation formula for the fiber delivery quantity m is as follows:

9. The method for establishing a mesoscopic finite element model of a reinforced ultra-high performance concrete structure according to claim 8, characterized in that: The method for determining the starting point and direction vector of the fiber in S5 is as follows: Determine the starting point and direction vector of the fiber, and use the rand function rand(1), which represents the generation of a random number between 0 and 1. The formula for calculating the starting point coordinates is as follows: Randomly generate the X coordinate of the starting point of each fiber: x1 = LX rand (1), The Y coordinate of the starting point of each fiber is randomly generated: y1 = LY rand (1), Randomly generate the Z coordinate of the starting point of each fiber: z1 = LZ rand (1), LX is the length of the cube, LY is the width of the cube, and LZ is the height of the cube; Use angle1 to represent the angle of the direction vector in the xy plane, and angle2 to represent the inclination angle of the direction vector in the zy plane. The calculation formula is as follows angle1=2·π·rand(1), angle2=2·π·rand(1), angle1 and angle2 represent random angles in the range of 0 to 2π; The fiber end point coordinates are calculated using the starting point coordinates and random angles angle1 and angle2. The calculation formula is as follows: The X coordinate of the end point of each fiber is randomly generated as: x2 = x1 + Ln·cos(angle2)·cos(angle1), The Y coordinate of the end point of each fiber is randomly generated as: y2 = y1 + Ln·cos(angle2)·sin(angle1), The Z coordinate of the end point of each fiber is randomly generated as: z2=z1+Ln·sin(angle2), When adjusting the fiber coordinates, when the starting point coordinates and the end point coordinates of the fiber exceed the minimum delivery boundary of the delivery area, the coordinates are automatically corrected to be offset 1 mm inward from the minimum delivery boundary; when the starting point coordinates and the end point coordinates exceed the maximum delivery boundary of the delivery area, the coordinates are automatically corrected to be offset 1 mm inward from the maximum delivery boundary.

10. The method for establishing a mesoscopic finite element model of a reinforced ultra-high performance concrete structure according to claim 9, characterized in that: When judging the intrusion between fiber and aggregate in S6, if intrusion is considered, the intruded fiber will be processed by translating it out of the aggregate. After all fibers are placed and invaded, the fiber placement range is judged to ensure that all fibers are within the placement range. When the starting point coordinates and the end point coordinates exceed the minimum placement boundary, the coordinates are automatically corrected to the minimum placement boundary with an inward offset of 1 mm. When the starting point coordinates and the end point coordinates exceed the maximum placement boundary, the coordinates are automatically corrected to the maximum placement boundary with an inward offset of 1 mm.

Images