Method for establishing microscopic finite element model of reinforced fiber concrete structure

Abstract

The present invention discloses a method for establishing a microscopic finite element model of a reinforced fiber concrete structure, belonging to the technical field of concrete microscopic modeling. By using the MATLAB language in conjunction with the ANSYS / LS-DYNA finite element software to perform simulation, a four-phase finite element model of a reinforced concrete structure is obtained. The method arranges steel bars that meet the requirements in a preset area, places aggregates and fibers in the placement area, and ensures that the aggregates, fibers, and steel bars do not invade each other. A large number of fibers are modeled, and directionally arranged steel bars are accurately modeled. Aggregates of different particle sizes are also placed therein, thereby solving the problem of mutual invasion between the fibers, steel bars, and aggregates. The simulation model is very close to the actual concrete structure, the model generation time is short, and the model accuracy is high.

Description

Technical Field

[0001] The present 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 often simplified as a homogeneous material at the macroscale for mechanical analysis, but in reality it is a typical heterogeneous material. Currently, many concrete models rely on macromechanical methods, ignoring the complex internal structural properties of concrete. This is applicable to a certain extent for the overall study of material properties. However, this simplified approach cannot fully reveal the internal complexity, deformation patterns, and failure mechanisms of concrete. Experimental studies have 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 mesoscale models of concrete, it is possible to deeply study the material's response under dynamic and static loads. This approach enables a more comprehensive understanding of concrete's actual properties and provides theoretical support for engineering design and architectural structure optimization. To ensure that mesoscale numerical simulations are more realistic, it is crucial to establish accurate mesoscale models. In this process, the numerical characterization of aggregate and the random distribution of fibers are important components of concrete mesoscale simulations. Therefore, the development of mesoscale models that accurately reflect the internal components and their distribution is crucial.

[0003] Currently, there are two main methods for modeling concrete mesoscale models. One is based on real concrete CT scan results, generating a corresponding concrete mesoscale model based on the CT image, such as the CT image recognition method used in invention patent applications (CN 105699258 A, CN105427382A); the other is through parametric modeling, such as the invention patent applications (CN118133493A, CN113987750A), which use programming languages ​​to generate three-dimensional spatial information of aggregate, mortar, and interface transition zones in concrete, and then import it into finite element simulation software to establish a mesoscale model.

[0004] Most current concrete mesoscopic models have the problem of a single component and fail to accurately characterize the various components in concrete. These models usually face several key problems: first, the aggregate gradation control 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, existing concrete mesoscopic 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 properties of concrete materials. Therefore, improving the accuracy of mesoscopic models in component characterization, gradation control, spatial randomness processing, and interactions between components is crucial to further promote concrete mesoscopic research and application. Summary of the Invention

[0005] The purpose of the present invention is to provide a method for establishing a microscopic finite element model of a reinforced fiber concrete structure that 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 by comprising 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] Calculate the length of the steel bars in the X and Y directions based on the thickness of the concrete cover in the X, Y, and Z directions and the length, width, and height of the cube. Calculate the spacing between the transverse and longitudinal bars based on the number of transverse and longitudinal bars arranged. Arrange the first transverse and longitudinal bars in the placement area and determine their coordinates. Complete the steel bar arrangement by translating them based on the spacing and save the 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. 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 step-by-step process of increasing in size first and decreasing in size later to simulate the aggregate particle size distribution in real concrete. This ensures that the aggregate spheres do not exceed the boundaries and that there is no intrusion or overlap between aggregates.

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

[0014] After the aggregate is pre-placed, the intrusion between the aggregate and the steel bar is determined, and the intruded aggregate is processed. The processed aggregate needs to be re-determined whether it has intrusion with other aggregates until all aggregates and all steel bars do not intrude or intersect, and the aggregate coordinates are saved. The processing methods include translation and diameter change.

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

[0016] Determine the minimum and maximum placement coordinate points of the X, Y, and Z axes, and then determine the fiber placement range. Determine the placement quantity based on the selected fiber size parameters and perform random placement. After placement, 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 fibers are placed and processed, the intrusion between fibers and aggregates is determined. When intrusion occurs between fibers and aggregates, the intruded fibers are handled by moving them out of the aggregates and ensuring that the moved fibers are located within the boundary.

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

[0020] After completing the intrusion determination and adjustment between all fibers and aggregates, perform the intrusion determination between fibers and steel bars, delete the fibers that have intruded into the steel bars, ensure that all fibers and all steel bars do not intrude or intersect, and save the fiber coordinates;

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

[0022] After all required fiber, aggregate, and reinforcement arrangements are completed, the coordinate information of the aggregate is converted to a finite element mesh node. This involves traversing each cell in the mesh through a triple loop, calculating the difference between the distance from the sphere center coordinate to the cell center point and the aggregate radius R. If the difference is less than 0, the cell is classified as this aggregate, marked, and renamed.

[0023] Convert the fiber and steel bar endpoint coordinates into finite element simulation nodes, renumber and arrange the node information so that it becomes unit information that can be recognized by the finite element software. After completing the arrangement of the steel bar unit information, define it 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] Use the fopen function in 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, 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 saved in end_node1 in matrix form, and the fiber ending point coordinates are saved 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 based on the relationship between R-1 and the aggregate radius R. When R-1 is less than the radius R, calculate the index n1 of this unit in the matrix Element and classify this unit as an aggregate. After completion, proceed to the unit replacement of the second aggregate until all aggregate units are replaced. Classify all units except the aggregate into a matrix, and save all updated unit node information in the Element in matrix form.

[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, and 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, and 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 sequence, 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 element information: Create a column vector BH1 from 1 to A+O, a column vector BH2 ​​from A+O+1 to A+B+2·O, and a column vector BH3 from A+B+2·O+1 to A+B+3·O for element numbering. Create a column vector PART of length A+B and define its values ​​from 1 to A as 1 for naming steel finite element components. Define its values ​​from A+1 to A+O as 2 for naming fiber finite element components. Combine PART, BH1, BH2, and BH3 to form the overall element 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]. 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 n1 is calculated as follows: n1=(i+1)+j·Nx+k·Nx·Ny,

[0043] Calculation formula for unit node 1: Element1 = (i+1) + j·(Nx+1) + k·(Nx+1)·(Ny+1),

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

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

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

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

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

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

[0050] Calculation formula for unit node 8: 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 the aggregates invade and overlap with each other. Every time a new aggregate is randomly generated, it needs to be judged for invasion 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 invade and overlap. If there is overlap, it needs to be regenerated and judged again;​​​​​​​​​​​​​​​​​​​​​​​​

[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 coordinates PointP (Cx, Cy, Cz) are extracted, and the starting 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 using the following formula:

[0075]

[0076]

[0077]

[0078]

[0079] Dist=|PointP-PointC|,

[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 and steel bar, we decide 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 steel bars have invaded, and a secondary judgment is made.

[0082] When Dist < (R + g - dia) · 3 / 5, the spatial position of the aggregate is translated. When the aggregate is translated, its radius is also regenerated to ensure that the aggregate after translation 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 input 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] To determine the starting point and direction vector of the fiber, use the rand function rand(1), which generates 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] Randomly generate the Y coordinate of the starting point of each fiber: 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 will be 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 will be 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 fibers and aggregates in S6, if intrusion is considered, the intruded fibers are processed by moving them 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 end point coordinates exceed the minimum placement boundary, the coordinates are automatically corrected to be offset 1 mm inward from the minimum placement boundary. When the starting point coordinates and end point coordinates exceed the maximum placement boundary, the coordinates are automatically corrected to be offset 1 mm inward from the maximum placement boundary.

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

[0117] Currently, most concrete micro-modeling technologies only include a single random fiber component or random aggregate component, and rarely describe all the 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 This is a finite element model diagram of the steel bars of the reinforced fiber concrete structure of the present invention;

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

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

[0121] Figure 4 This is a 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 steel bar aggregate in 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 clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. 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 making any creative efforts shall fall 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 is a model constructed by taking a single-layer reinforced fiber concrete structure as an example to describe the embodiment 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 within the designated area. Since the geometry of steel bars approximates line segments, the directional line segments can be considered steel bars. Based on this, random fibers are generated. Since the geometry of fibers approximates line segments, the random line segments can be considered fibers. The specific process is as follows:

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

[0130] Take a cube as an example to determine the size of the 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.

[0131] Rebar related parameter settings, X-axis concrete cover thickness px, Y-axis concrete cover thickness py, Z-axis concrete cover thickness pz, steel bar cross-section diameter G-dia, then 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 the 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 end point coordinates 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 reinforcements, first calculate the starting coordinates of the first longitudinal reinforcement in the placement area, calculate the end coordinates of the first longitudinal reinforcement using the reinforcement-related parameters, save the three coordinates of the reinforcement in the Cartesian coordinate system in matrix form, translate the starting coordinates and end coordinates of the first longitudinal reinforcement, and obtain the starting coordinates and end coordinates of the second longitudinal reinforcement until the arrangement of the longitudinal reinforcement 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 layout of all steel bars, save the starting coordinates of the horizontal and vertical steel bars in matrix form in g_node1, and save the ending coordinates of the horizontal and vertical steel bars in matrix form in g_node2.

[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 elements 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 divided evenly according to the unit size in each direction. i, j, k represent the index 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 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 index of the x, y, z direction 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 of n1 is as follows: n1=(i+1)+j·Nx+k·Nx·Ny,

[0173] Calculation formula for unit node 1: Element1 = (i+1) + j·(Nx+1) + k·(Nx+1)·(Ny+1),

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

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

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

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

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

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

[0180] Calculation formula for unit node 8: 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 reinforcement layout and unit grid division, the aggregate model is placed. The aggregate is placed in a step-by-step manner from large to small. By controlling the different particle sizes of aggregate, the aggregate particle size distribution 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. Each time a new aggregate is randomly generated, it needs to be judged for 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 needs to be regenerated and judged again.

[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 steel bars

[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 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-PointC|,

[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 and steel bar, we decide how to translate the aggregate. The judgment conditions are as follows:

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

[0212] When Dist < (R + g - dia) · 3 / 5, the spatial position of the aggregate is translated. When the aggregate is translated, its radius is also regenerated to ensure that the aggregate after translation 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) need to be re-judged to see if they intrude into other aggregates, until all aggregates and all steel bars do not intrude 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 on the X axis of the delivery range, y0 is the minimum delivery coordinate point on the Y axis of the delivery range, z0 is the minimum delivery coordinate point on the Z axis of the delivery range, Lx1 is the maximum delivery coordinate point on the X axis of the delivery range, Ly1 is the maximum delivery coordinate point on the Y axis of the delivery range, Lz1 is the maximum delivery coordinate point on the Z axis of the delivery range, 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 input quantity m is as follows:

[0230] To determine the starting point and direction vector of the fiber, use the rand function rand(1), which generates 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] Randomly generate the Y coordinate of the starting point of each fiber: 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 fibers are placed, their starting and ending coordinates are determined. When the starting and ending coordinates exceed the minimum placement boundary of the placement area, their coordinates are automatically corrected to be offset 1 mm inward 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 fibers are placed and processed, the fiber-aggregate intrusion determination is performed. The following is one specific method, but other methods or nomenclatures are also possible. The fiber-rebar intrusion determination is also a specific method, but other methods are also possible.

[0259] Extract the center coordinates PointP (Cx, Cy, Cz) of the first spherical aggregate and the corresponding sphere radius R. Extract the starting coordinates C (X11, Y11, Z11) and end coordinates D (X12, Y12, Z12) of the first fiber. Calculate the shortest distance from the center of the spherical aggregate to the fiber and determine its relationship with the radius R of the spherical aggregate 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 relationship between Distance and the sphere radius R, it is determined whether the fiber and aggregate have invaded. When the fiber and aggregate have invaded, the invading fiber is moved out of the aggregate to solve the problem. The specific calculation formula is as follows:

[0266]

[0267]

[0268]

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

[0270]

[0271]

[0272] NewX21=PointGx-0.5·n·T2x,

[0273] NewY21=PointGy-0.5·n·T2y,

[0274] NewZ21=PointGz-0.5·n·T2z,

[0275] NewX22=PointGx+0.5·n·T2x,

[0276] NewY22=PointGy+0.5·n·T2y,

[0277] NewZ22=PointGz+0.5·n·T2z,

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

[0279] NewX21=x-min+0.001,

[0280] NewY21=y-min+0.001,

[0281] NewZ21=z-min+0.001,

[0282] NewX22=x-min+0.001,

[0283] NewY22=y-min+0.001,

[0284] NewZ22=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 offset 1 mm inward from the maximum delivery boundary. The calculation formula is as follows:

[0286] NewX21=x-max-0.001,

[0287] NewY21=y-max-0.001,

[0288] NewZ21=z-max-0.001,

[0289] NewX22=x-max-0.001,

[0290] NewY22=y-max-0.001,

[0291] NewZ22=z-max-0.001,

[0292] Step 3.3: Fiber and steel intrusion determination

[0293] After completing the invasion judgment between all fibers and aggregates, the invasion judgment between fibers and steel bars is performed, and the starting coordinates of the first fiber are extracted: NewPOINT1 (NewX21, NewY21, NewZ21).

[0294] With the end point coordinates: NewPOINT2 (NewX22, NewY22, NewZ22),

[0295] Extract the starting coordinates of the first steel bar POINT-G1(X G 1,Y G 1,Z G 1) and the end point coordinate POINT-G1(X G 2,Y G 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 endpoint 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 coordinates of the fiber that meets the requirements are saved in end_node1 in matrix form, and the ending 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 all required fiber, aggregate, and rebar arrangements are completed, the aggregate coordinate information is converted to finite element mesh nodes. This involves traversing each cell in the mesh using a triple loop, calculating the difference between the distance from the sphere's center coordinate to the cell's center and the radius R. If the difference is less than 0, the cell is classified as this aggregate, labeled, and renamed. The fiber and rebar endpoint coordinates are converted to finite element simulation nodes, and the node information is renumbered and reorganized to become unit information recognizable by the finite element software. After the rebar unit information is organized, it is defined 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 is used to save the information as a keyword K file recognizable 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.

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

[0319] Extract the spatial coordinates of the first aggregate and calculate the distance R-1 between each cell center and the center of the aggregate sphere. Based on the relationship between R-1 and the aggregate radius R, determine the cell type. If R-1 is less than the radius R, calculate the cell's index n1 in the Element matrix and classify it as an aggregate. Once completed, proceed to the cell replacement for the second aggregate, and continue until all aggregates have been replaced. All cells except the aggregate are classified as a matrix, and all updated cell node information is stored in the Element matrix as a matrix.

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

[0321] Node information arrangement: Create a column vector BH4 with a length equal to the number of rebars, O, and set its value to 0. Concatenate BH4 with g_node1 and g_node2 to form a longer node vector. Extract the row number A of end_node1 (i.e., the number of nodes in end_node1) and the row number B of end_node2 (i.e., the number of nodes in end_node2), and create a column vector BH ranging from 1 to A+B+O (the sum of the number of rebars and fibers) for subsequent node numbering. Arrange the node information in end_node1 and end_node2 in order, and concatenate 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 node coordinates and other information. Arrange the processed node information in the desired order to obtain the total node set.

[0322] Element information arrangement: Create a column vector BH1 from 1 to A+O, a column vector BH2 ​​from A+O+1 to A+B+2·O, and a column vector BH3 from A+B+2·O+1 to A+B+3·O for element numbering. Create a column vector PART of length A+B and define its values ​​from 1 to A as 1 for naming steel finite element components. Define its values ​​from A+1 to A+O as 2 for naming fiber finite element components. Combine PART, BH1, BH2, and BH3 to form the overall element 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 embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any technician familiar with the technical field, within the technical scope disclosed by the present invention, who makes equivalent replacements or changes based on the technical solution and inventive concept of the present invention, should be covered by the scope of protection 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. The characteristics are as follows: 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 layout by shifting the reinforcement 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. 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 step-by-step process of increasing in size first and decreasing in size later to simulate the aggregate particle size distribution in real concrete. This ensures that the aggregate spheres do not exceed the boundaries and that there is no intrusion or 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 they are within the boundaries. S6. Determine and adjust the intrusion between fiber and aggregate: When fibers and aggregates invade each other, the intruding fibers are handled by moving them out of the aggregate and ensuring that the moved fibers are located inside the boundary. S7. Determine and adjust the penetration between fiber and steel bar: Determine the intrusion between fibers and steel bars, delete fibers that intrude into steel bars, ensure that all fibers and all steel bars do not intrude or intersect, and save the fiber coordinates; S8. Constructing a microscopic finite element model of reinforced concrete structure: After all required fibers, aggregates, and steel bars are arranged, the coordinate information of the aggregate is converted into finite element mesh nodes. This involves traversing each cell in the mesh through a triple loop, calculating the difference between the distance from the sphere center coordinate to the cell center point and the aggregate radius. If the difference is less than 0, the cell is classified as this aggregate, marked, and renamed. Convert the fiber and steel bar endpoint coordinates into finite element simulation nodes, renumber and arrange the node information so that it becomes unit information that can be recognized by the finite element software. After completing the arrangement of the steel bar unit information, define it 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; Use the fopen function in 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, 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, wherein: 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 saved in end_node1 in matrix form, and the fiber ending point coordinates are saved 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, wherein: When the aggregate coordinate information in S8 is converted into 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 aggregate radius R. When R-1 is less than the radius R, the index n1 of this unit in the matrix Element is calculated, and this unit is divided into aggregate. After completion, the unit replacement of the second aggregate is entered until the unit replacement of all aggregates is completed. The units except the aggregate are divided into the matrix, and all updated unit node information is saved in the Element in matrix form. 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, and 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, and 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 sequence, 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 element 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 naming steel finite element components, define its values ​​from A+1 to A+O as 2, for naming fiber finite element components, and combine PART, BH1, BH2, and BH3 into the overall element matrix.

4. The method for establishing a microscopic finite element model of a reinforced concrete structure according to claim 3, wherein: 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]. 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 microscopic finite element model of a reinforced concrete structure according to claim 4, wherein: 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 n1 is calculated as follows: n1=(i+1)+j·Nx+k·Nx·Ny, Calculation formula for unit node 1: Element1 = (i+1) + j·(Nx+1) + k·(Nx+1)·(Ny+1), Calculation formula for unit node 2: Element2 = (i+2)+j·(Nx+1)+k·(Nx+1)·(Ny+1), Calculation formula for unit node 3: Element3 = (i+2)+(j+1)·(Nx+1)+k·(Nx+1)·(Ny+1), Calculation formula for unit node 4: Element4 = (i+1)+(j+1)·(Nx+1)+k·(Nx+1)·(Ny+1), Calculation formula for unit node 5: Element5 = (i+1)+j·(Nx+1)+(k+1)·(Nx+1)·(Ny+1), Calculation formula for unit node 6: Element6 = (i+2)+j·(Nx+1)+(k+1)·(Nx+1)·(Ny+1), Calculation formula for unit node 7: Element7 = (i+2)+(j+1)·(Nx+1)+(k+1)·(Nx+1)·(Ny+1), Calculation formula for unit node 8: 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 microscopic finite element model of a reinforced ultra-high performance concrete structure according to claim 5, wherein: 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 formula for calculating 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 formula for calculating the aggregate radius is as follows: R = 1 / 2·Maxsize - a + rand(1)·1 / 6·Maxsize - a, When Vpersen·11 / 12 < Asize, the formula for calculating the aggregate radius is as follows: R = Minsize - a + rand(1)·Minsize - a, 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, which 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 placing the coordinates of the aggregate center, it is necessary to judge whether the aggregates invade and overlap with each other. For each randomly generated aggregate, it needs to be judged for invasion 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 invade and overlap. If overlap occurs, it needs to be regenerated and judged again; 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 (Cx, Cy, Cz) of the aggregate center. The judgment conditions are as follows: R < A1, (Cx + R) ≤ Lx, (Cy + R) ≤ Ly, ​ ​ Vtotal=4 / 3·π·R 3 。 ​ In S4, when determining the intrusion of aggregate and steel bars, the radius R of the spherical aggregate and its corresponding spherical center coordinates PointP (Cx, Cy, Cz) are extracted, and the starting 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 using the following formula: ​ ​ ​ ​ ​ ​ ​ NewR=0.64·(g-dia+R).

8. The method for establishing a microscopic finite element model of a reinforced ultra-high performance concrete structure according to claim 7, wherein: 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 input quantity m is as follows:

9. The method for establishing a microscopic finite element model of a reinforced ultra-high performance concrete structure according to claim 8, wherein: The method for determining the starting point and direction vector of the fiber in S5 is as follows: To determine the starting point and direction vector of the fiber, use the rand function rand(1), which generates 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), Randomly generate the Y coordinate of the starting point of each fiber: 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 will be 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 will be automatically corrected to be offset 1 mm inward from the maximum delivery boundary.

10. The method for establishing a microscopic finite element model of a reinforced ultra-high performance concrete structure according to claim 9, wherein: When judging the intrusion between fibers and aggregates in S6, if intrusion is considered, the intruded fibers are processed by moving them 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 end point coordinates exceed the minimum placement boundary, the coordinates are automatically corrected to be offset 1 mm inward from the minimum placement boundary. When the starting point coordinates and end point coordinates exceed the maximum placement boundary, the coordinates are automatically corrected to be offset 1 mm inward from the maximum placement boundary.

Images