Method for establishing mesoscopic finite element model of double-layer reinforced concrete structure

Abstract

The invention discloses a method for establishing a mesoscopic finite element model of a double-layer reinforced concrete structure, and belongs to the technical field of mesoscopic simulation of concrete structures. According to the method, steel bars meeting requirements are arranged in a preset area, and finite element grid node processing is conducted on an aggregate putting area. According to the method, grid division is directly completed in the aggregate generation process, so that an overall model is constructed, and the model adopts eight-node hexahedron units and has high similarity with a real reinforced concrete structure. Compared with other methods, the modeling time of the model is short, the convergence of a calculation result is good, and the method is suitable for various static and dynamic test simulations.

Description

Technical Field

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

[0002] As a widely used building material, concrete structure has the advantages of economy, easy availability and simple construction. At present, most concrete structure models rely on macroscopic mechanical methods. Although this method meets the research needs of overall material properties, it often ignores the complexity of the concrete structure. This simplified model cannot fully reveal the complex structure inside the concrete and its deformation and failure modes. Due to the limitations of experimental research, it is challenging to accurately observe the 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 structure, finite element simulation can conduct in-depth analysis of the material's response under dynamic and static loads, which is crucial for a comprehensive understanding of the performance of concrete structures and provides important theoretical support for engineering design and building structure optimization. Summary of the invention

[0003] The purpose of the present invention is to provide a method for establishing a mesoscopic finite element model of reinforced concrete structure which is easy to operate and has high accuracy. The method can accurately observe the deformation and damage inside the concrete through the established model, improve the shortcomings of traditional concrete detection, and provide important theoretical support for engineering design and building structure optimization.

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

[0005] The method for establishing a microscopic finite element model of a double-layer reinforced 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:

[0006] 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. The volume of the concrete is:

[0007] V = LX·LY·LZ;

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

[0009] 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, the length of the steel bars in the X-axis and Y-axis directions is calculated; according to the number of transverse and longitudinal steel bars arranged, the spacing of the transverse and longitudinal steel bars is calculated; the first transverse and longitudinal steel bars are arranged in the placement area and the coordinates are determined, the first layer of steel bars are arranged according to the spacing, the first layer of transverse and longitudinal steel bars are translated to obtain the second layer of steel bars and the coordinates are determined, and the vertical steel bars are placed after the transverse and longitudinal steel bars are arranged. By judging the odd-even relationship of the sum of the transverse and longitudinal steel bar numbers, it is decided whether to place the vertical steel bars at the intersection of the transverse and longitudinal steel bars, so as to ensure that the vertical steel bars are arranged alternately. The placement of the vertical steel bars is still carried out by coordinate translation, and the steel bar coordinates are saved;

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

[0011] 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.

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

[0013] After the pre-placement of aggregates is completed, the intrusion between the aggregates and the steel bars is judged, and the intruded aggregates are processed by translation, including translation and diameter change and regeneration. The processed aggregates need to be re-judged whether they have intrusion with other aggregates, until there is no intrusion and intersection with all aggregates and all steel bars, and the next aggregate information is extracted, and the judgment is continued until the intrusion judgment and processing of all aggregates and all steel bars are completed;

[0014] S5. Construct a microscopic finite element model of reinforced concrete structure:

[0015] After all the aggregates that meet the requirements are placed, each cell in the grid is traversed through a triple loop to calculate the difference between the distance from the sphere center coordinate to the cell center point and the aggregate radius. When the value is less than 0, the cell is divided into this aggregate, and it is marked and renamed with element update; the coordinates of the steel bar endpoints 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 the arrangement of the steel bar unit information is completed, it is defined as part-1, the aggregate unit information is defined as part-2, and the matrix unit information is defined as part-3; the fopen function in the 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, and the establishment of the microscopic finite element model is completed.

[0016] As a further description of the above technical solution:

[0017] The total number of reinforcement arrangements in S2 is O:

[0018] O=2·(nx+ny)+nz,

[0019] In the formula, nx is the number of longitudinal reinforcements, ny is the number of transverse reinforcements, and nz is the number of vertical reinforcements;

[0020] After all the steel bars are arranged in S2, the matrices storing the starting coordinates of the horizontal, longitudinal and vertical steel bars are concatenated into g_node1, and the matrices storing the ending coordinates of the horizontal, longitudinal and vertical steel bars are concatenated into g_node2.

[0021] As a further description of the above technical solution:

[0022] When arranging the node information in S5, 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;

[0023] Extract the row number D of node and create a column vector BH from 1 to D+O for subsequent node numbering;

[0024] Initialize a zero matrix Nodea of ​​size (2·O, 5) to store the coordinates and other information of the steel bar nodes. Arrange the processed node information in the required order to obtain the total node set, that is, fill the node information of g_node1 into rows 1 to O of the zero matrix Nodea, and fill the node information of g_node2 into rows 1+O to 2·O of the zero matrix Nodea.

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

[0026] When arranging the unit information in S5, create a column vector BH1 from D+1 to D+O, create a column vector BH2 ​​from D+O+1 to D+2·O, create a column vector BH3 from D+2·O+1 to D+3·O for element numbering, create a column vector PART with a length of O, define its row values ​​from 1 to O as 1, and combine PART, BH1, BH2, and BH3 into an overall unit matrix.

[0027] As a further description of the above technical solution:

[0028] The formula for calculating the delivery range in S3 is as follows:

[0029] Lx=Nx·Msx,

[0030] Ly=Ny·Msy,

[0031] Lz=Nz·Msz,

[0032] The total number of units Nmesh is calculated as follows:

[0033] Nmesh=Nx·Ny·Nz,

[0034] The calculation formula for the total number of nodes Nnoda is as follows:

[0035] Nnoda=(Nx+1)·(Ny+1)·(Nz+1),

[0036] Wherein, 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, and Msz is the unit size 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, 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:

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

[0039] Store the global number n of the current node in the first column of Node0, store the x coordinate of the current node in the second column, store the y coordinate of the current node in the third column, store the z coordinate of the current node in the fourth column, and perform a linear transformation on the matrix Node0 to adjust the node coordinates:

[0040]

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

[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 for n1 is as follows:

[0043] n1=(i+1)+j·Nx+k·Nx·Ny,

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

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

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

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

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

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

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

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

[0052] 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.

[0053] As a further description of the above technical solution:

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

[0055] The initial value of Asize is 0. The rand function rand(1) is used, which represents the generation of a random number between 0 and 1. 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 aggregate radius calculation formula is as follows:

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

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

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

[0061] In the formula, Maxsize - a is the maximum aggregate radius, Minsize - a is the minimum aggregate radius, and Vpersen is the aggregate volume content.

[0062] As a further description of the above technical solution:

[0063] After calculating the aggregate radius, the aggregate sphere center coordinates (Cx, Cy, Cz) are placed. The coordinates are randomly generated within the restricted range to ensure that the sphere does not exceed the boundary, and the boundary is defined by Lx, Ly, Lz.

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

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

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

[0067] After placing the sphere center coordinates, it is necessary to judge whether the aggregates invade and overlap with each other. Every time a random aggregate is 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 sphere 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 rejudged.

[0068] After placing 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:

[0069] R < A1,

[0070] (Cx + R) ≤ Lx,

[0071] (Cy + R) ≤ Ly,

[0072] (Cz + R) ≤ Lz,

[0073] The volume formula of aggregate is:

[0074] Vtotal=4 / 3·π·R 3 .

[0075] As a further description of the above technical solution:

[0076] In S4, when determining the intrusion of aggregate and steel bars, the center coordinates of the spherical aggregate (PointA, Cx, Cy, Cz) are extracted, and the starting coordinates of the first steel bar (Point-G1, X G 1,Y G 1,Z G 1) and the end point coordinates Point-G2(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:

[0077]

[0078]

[0079]

[0080]

[0081] Dist=|PointA-Point B|,

[0082] 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, it is decided how to translate the aggregate.

[0083] As a further description of the above technical solution:

[0084] Dist judgment conditions are as follows:

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

[0086] When Dist<(g-dia+0.01)·4 / 7, 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:

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

[0088]

[0089]

[0090] NewPointA=PointA+direction·0.4·(g-dia+R),

[0091] When (g-dia+0.01)·4 / 7≤Dist<(g-dia+0.01), a new radius is generated:

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

[0093] After the aggregate has been processed by translation, it is necessary to re-judge whether it has invaded other steel bars until there is no invasion and intersection with all steel bars. The next aggregate information is extracted and the judgment is continued until the intrusion judgment and processing of all aggregates and all steel bars are completed.

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

[0095] At present, most concrete mesoscopic modeling technologies only include a single aggregate component and a steel component, and rarely include all the components of reinforced concrete. This method models a large number of aggregates, and the aggregates also simulate the particle gradation of real concrete. At the same time, the steel bars arranged in two directions are accurately modeled, solving the intrusion defects between aggregates and steel bars. The simulation model is closer to the actual reinforced concrete structure, the model generation time is short, and the model accuracy is high. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0097] Figure 2 It is a finite element model diagram of the steel bar arrangement of the reinforced concrete structure of the present invention;

[0098] Figure 3 It is a finite element model diagram of the spatial arrangement of steel bar aggregates of the present invention;

[0099] Figure 4 It is a finite element model diagram of the reinforced concrete structure section of the present invention. DETAILED DESCRIPTION

[0100] 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.

[0101] The present invention obtains a reinforced concrete structure finite element model by simulating and emulating using MATLAB language in conjunction with ANSYS / LS-DYNA finite element software.

[0102] The following takes a double-layer reinforced concrete structure as an example to construct a model and describe the implementation of the present invention in detail. The constructed finite element model diagram can be found in Figure 1-4 The detailed steps are as follows:

[0103] Step 1: Generate basic information of reinforced concrete structure in MATLAB software

[0104] Generate two layers of horizontal and vertical line segments in the specified area. Since the geometric shape of the steel bar is similar to the line segment, the directional line segment can be considered as the steel bar. Then add the unit section attribute to it in the ANSYS / LS-DYNA finite element software to complete the modeling of the steel bar. On this basis, the specified delivery area is meshed, and spheres of different particle size ranges are randomly generated in the specified area. Since the geometric shape of the aggregate is similar to a sphere, the random sphere can be considered as the aggregate. Use the mapping grid algorithm to mark the grid occupied by the sphere as aggregate, and the rest as the matrix. The intrusion judgment of the aggregate and the steel bar is performed to ensure that the aggregate and the steel bar do not intersect, and the reinforced concrete micro-finite element model is obtained. The specific process is as follows:

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

[0106] 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. The volume of the concrete is:

[0107] V=LX·LY·LZ,

[0108] 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:

[0109] G-lx=LX-2·px,

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

[0111] 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;

[0112] Declare the spacing between the two layers of steel bars JZ, determine the number of longitudinal steel bars nx, the number of transverse steel bars ny, the number of vertical steel bars nz, and calculate the total number of steel bars O:

[0113] O=2·(nx+ny)+nz,

[0114] Calculate the longitudinal reinforcement spacing JX and the transverse reinforcement spacing JY:

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

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

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

[0118] First, specify the starting coordinates of the first transverse steel bar in the first layer in the placement area, calculate the end coordinates of the first 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 transverse steel bar to obtain the starting coordinates and end coordinates of the second transverse steel bar, translate the first layer of steel bars to obtain the corresponding coordinates of the second layer of steel bars, until the arrangement of the transverse steel bars is completed.

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

[0120] X11 g =px,

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

[0122] Z11 g =pz,

[0123] X11 g 、Y11 g and Z11 g Indicates the coordinates of the starting point of the first layer of transverse reinforcement (for easy distinction, the first number reflects the starting point or end point, and the second number reflects the number of layers. The starting point of the transverse reinforcement is represented by the number 1, and the end point is represented by the number 2. The starting point of the longitudinal reinforcement is represented by the number 3, and the end point is represented by the number 4. The starting point of the vertical reinforcement is represented by the number 5, and the end point is represented by the number 6. The vertical reinforcement is not layered, so there is no second number).

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

[0125] X21 g =px+G-lx,

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

[0127] Z21 g =pz,

[0128] X21 g 、Y21 gand Z21 g Indicates the coordinates of the end point of the first layer of transverse reinforcement;

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

[0130] X12 g =px,

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

[0132] Z12 g =pz+JZ,

[0133] X12 g 、Y12 g 、Z12 g Indicates the starting coordinates of the second layer of transverse reinforcement.

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

[0135] X22 g =px+G-lx,

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

[0137] Z22 g =pz+JZ,

[0138] X22 g 、Y22 g and Z22 g Indicates the coordinates of the end point of the second layer of transverse reinforcement;

[0139] After completing the arrangement of all transverse steel bars, first calculate the starting coordinates of the first longitudinal steel bar in the first layer 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 the form of a matrix, translate the starting coordinates and end coordinates of the first longitudinal steel bar to obtain the starting coordinates and end coordinates of the second longitudinal steel bar, translate the first layer of steel bars to obtain the corresponding coordinates of the second layer of steel bars, until the arrangement of the longitudinal steel bars is completed.

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

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

[0142] Y31 g =py,

[0143] Z31g =pz+G-dia,

[0144] X31 g 、Y31 g and Z31 g Indicates the starting coordinates of the first layer of longitudinal reinforcement;

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

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

[0147] Y41 g =py+G-ly,

[0148] Z41 g =pz+G-dia,

[0149] X41 g 、Y41 g and Z41 g Indicates the starting coordinates of the first layer of longitudinal reinforcement;

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

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

[0152] Y32 g =py,

[0153] Z32 g =pz+G-dia+JZ,

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

[0155] X42 g =px+G-lx,

[0156] Y42 g =py+JY·(i-1)i∈[1,nx],

[0157] Z42 g =pz+G-dia+JZ,

[0158] After completing the arrangement of horizontal and longitudinal reinforcements, vertical reinforcements are placed. First, calculate the starting coordinates of the first vertical reinforcement, and then obtain the coordinates of the remaining vertical reinforcements through translation.

[0159] Before placing vertical reinforcement, use the mod function to determine in advance whether (i+j) is an odd number. If it is an odd number, do not place vertical reinforcement at this coordinate point. For example, when i is 1 and j is 1, it means the intersection of the first longitudinal reinforcement and the first transverse reinforcement. At this time, (i+j) is 2, and the vertical reinforcement is placed here. When i is 1 and j is 2, it means the intersection of the first longitudinal reinforcement and the second transverse reinforcement. At this time, (i+j) is 3, which is an odd number, and the reinforcement is not placed here. This method is used to ensure the alternating arrangement of the vertical reinforcement placement. The placement of vertical reinforcement is still carried out by coordinate translation.

[0160] The calculation formula for the starting point of vertical reinforcement is as follows:

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

[0162] Y5 g =py+JY·(j-1)j∈[1,ny],

[0163] Z5 g =pz,

[0164] X5 g 、Y5 g and Z5 g Indicates the starting coordinates of the vertical reinforcement;

[0165] The calculation formula for the end point of vertical reinforcement is as follows:

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

[0167] Y6 g =py+JY·(i-1)j∈[1,ny],

[0168] Z6 g =pz+G-dia+JZ,

[0169] X6 g 、Y6 g and Z6 g Indicates the end point coordinates of the vertical reinforcement;

[0170] Step 1.3: Save the rebar coordinates

[0171] After completing the arrangement of all steel bars, the matrices storing the starting coordinates of the transverse, longitudinal and vertical steel bars are spliced ​​into g_node1, and the matrices storing the ending coordinates of the transverse, longitudinal and vertical steel bars are spliced ​​into g_node2.

[0172] Step 1.4: Aggregate placement area and related parameter settings

[0173] 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:

[0174] Lx=Nx·Msx,

[0175] Ly=Ny·Msy,

[0176] Lz=Nz·Msz,

[0177] The total number of units is calculated as follows:

[0178] Nmesh=Nx·Ny·Nz,

[0179] The total number of nodes is calculated as follows:

[0180] Nnoda=(Nx+1)·(Ny+1)·(Nz+1),

[0181] 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:

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

[0183] Store the global number n of the current node in the first column of Node0, store the x coordinate of the current node in the second column, store the y coordinate of the current node in the third column, store the z coordinate of the current node in the fourth column, and perform a linear transformation on the matrix Node0 to adjust the node coordinates:

[0184]

[0185] 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:

[0186] n1=(i+1)+j·Nx+k·Nx·Ny,

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

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

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

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

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

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

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

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

[0195] 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.

[0196] Step 1.5: Aggregate model placement

[0197] 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:

[0198] The initial value of Asize is 0. The rand function rand(1) is used, which represents generating 1 random number between 0 and 1. When Asize < (Vpersen·3 / 4), the calculation formula for the aggregate radius is as follows:

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

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

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

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

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

[0204] After calculating the aggregate radius, the aggregate sphere center coordinates (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, Lz.

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

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

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

[0208] After placing the aggregate sphere center coordinates, it is necessary to judge whether there is intrusion and overlap between the aggregates. Every time an 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 sphere 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 rejudged.

[0209] After placing 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). The judgment conditions are as follows:

[0210] R < A1,

[0211] (Cx + R) ≤ Lx,

[0212] (Cy+R)≤Ly,

[0213] (Cz+R)≤Lz,

[0214] The volume formula of aggregate is:

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

[0216] Step 1.6: Determine the penetration between aggregate and reinforcement

[0217] 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 PointA (Cx, Cy, Cz) are extracted, and the starting point coordinates Point-G1 (X G 1,Y G 1,Z G 1) and the end point coordinates Point-G2(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 B is the coordinate of the point on the steel bar where the shortest distance to the aggregate is located. The calculation formula is as follows:

[0218]

[0219]

[0220]

[0221]

[0222] Dist=|PointA-Point B|,

[0223] 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:

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

[0225] When Dist<(g-dia+0.01)·4 / 7, 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:

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

[0227]

[0228]

[0229] NewPointA=PointA+direction·0.4·(g-dia+R),

[0230] When (g-dia+0.01)·4 / 7≤Dist<(g-dia+0.01), a new radius is generated:

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

[0232] After the aggregate is processed by translation, etc. (including translation, diameter change and regeneration, etc.), it is necessary to re-judge whether it has invaded other steel bars until there is no invasion and intersection with all steel bars. The next aggregate information is extracted and the judgment is continued until the intrusion judgment and processing of all aggregates and all steel bars are completed.

[0233] Step 2: Construct a microscopic finite element model of reinforced concrete structure

[0234] After all the aggregates that meet the requirements are placed, each cell in the grid is traversed through a triple loop to calculate the difference between the distance from the sphere center coordinate to the cell center point and the radius R (referring to the aggregate radius, which can also be NewR). When its value is less than 0, its cell is divided into this aggregate, and it is marked and renamed with element update. The coordinates of the steel bar endpoints 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 the arrangement of the steel bar unit information is completed, it is defined as part-1, the aggregate unit information is defined as part-2, and the matrix unit information is defined as part-3. 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.

[0235] Step 2.1: Convert aggregate coordinate information into finite element simulation model information

[0236] 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.

[0237] Step 2.2: Convert the steel bar coordinate information into finite element simulation model information

[0238] Arrange the 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 number of rows D of node (that is, the number of nodes of node), create a column vector BH from 1 to D+O (the sum of the number of steel bar nodes and the number of aggregate nodes), which is used for subsequent node numbering. Initialize a zero matrix Nodea of ​​size (2·O, 5) to store the coordinates and other information of the steel bar nodes, arrange the processed node information in the required order to obtain the total node set, that is, fill the node information of g_node1 into the 1 to O rows of the zero matrix Nodea, and fill the node information of g_node2 into the 1+O to 2·O rows of the zero matrix Nodea.

[0239] Arrange the cell information, create a column vector BH1 from D+1 to D+O, create a column vector BH2 ​​from D+O+1 to D+2·O, create a column vector BH3 from D+2·O+1 to D+3·O, for element numbering, create a column vector PART of length O, define its row values ​​1 to O as 1. Combine PART, BH1, BH2, and BH3 into the overall cell matrix.

[0240] Step 2.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.

[0241] Step 2.4: Use the finite element analysis software ANSYS / LS-DYNA to directly open the obtained K file, where PART1 is the steel bar, PART2 is the aggregate, and PART3 is the matrix.

[0242] 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 the double-layer reinforced concrete structure is to obtain the finite element model of the reinforced concrete structure by simulating with MATLAB language and ANSYS / LS-DYNA finite element software. The characteristics are as follows: The steps include: S1. Determine the size and volume of cubic concrete: S2. Complete the setting and calculation of relevant parameters of the steel bars: arrange the first transverse and longitudinal steel bars in the delivery area and determine the coordinates, complete the arrangement of the first layer of steel bars according to the spacing, translate the first layer of transverse and longitudinal steel bars to obtain the second layer of steel bars, arrange the vertical steel bars alternately, and save the steel bar 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, and perform other processing such as translation on the intruded aggregates, ultimately ensuring that all aggregates and all steel bars do not intrude; S5. Construct a microscopic finite element model of reinforced concrete structure: After all the aggregates that meet the requirements are placed, each cell in the grid is traversed through a triple loop to calculate the difference between the distance from the sphere center coordinate to the cell center point and the aggregate radius. When the value is less than 0, the cell is divided into this aggregate, and it is marked and renamed with element update; the coordinates of the steel bar endpoints 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 the arrangement of the steel bar unit information is completed, it is defined as part-1, the aggregate unit information is defined as part-2, and the matrix unit information is defined as part-3; the fopen function in the 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, and the establishment of the microscopic finite element model is completed.

2. The method for establishing a mesoscopic finite element model of a double-layer reinforced concrete structure according to claim 1, characterized in that: After all the steel bars are arranged in S2, the matrices storing the starting coordinates of the horizontal, longitudinal and vertical steel bars are concatenated into g_node1, and the matrices storing the ending coordinates of the horizontal, longitudinal and vertical steel bars are concatenated into g_node2.

3. The method for establishing a microscopic finite element model of a double-layer reinforced ultra-high performance concrete structure according to claim 2, characterized in that: When arranging the node information in S5, 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 D of node and create a column vector BH from 1 to D+O for subsequent node numbering; Initialize a zero matrix Nodea of ​​size (2·O, 5) to store the coordinates and other information of the steel bar nodes. Arrange the processed node information in the required order to obtain the total node set, that is, fill the node information of g_node1 into rows 1 to O of the zero matrix Nodea, and fill the node information of g_node2 into rows 1+O to 2·O of the zero matrix Nodea.

4. The method for establishing a microscopic finite element model of a double-layer reinforced ultra-high performance concrete structure according to claim 3, characterized in that: When arranging the unit information in S5, create a column vector BH1 from D+1 to D+O, create a column vector BH2 ​​from D+O+1 to D+2·O, create a column vector BH3 from D+2·O+1 to D+3·O for element numbering, create a column vector PART with a length of O, define its row values ​​from 1 to O as 1, and combine PART, BH1, BH2, and BH3 into an overall unit matrix.

5. The method for establishing a mesoscopic finite element model of a double-layer reinforced concrete structure according to claim 4, characterized in that: The formula for calculating the delivery range in S3 is as follows: Lx=Nx·Msx, Ly=Ny·Msy, Lz=Nz·Msz, The total number of units Nmesh is calculated as follows: Nmesh=Nx·Ny·Nz, The calculation formula for the total number of nodes Nnoda is as follows: Nnoda=(Nx+1)·(Ny+1)·(Nz+1), Wherein, 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, and Msz is the unit size 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, 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: n=(i+1)+j·(Nx+1)+k·(Nx+1)·(Ny+1), Store the global number n of the current node in the first column of Node0, store the x coordinate of the current node in the second column, store the y coordinate of the current node in the third column, store the z coordinate of the current node in the fourth column, and perform a linear transformation on the matrix Node0 to adjust the node coordinates:

6. The method for establishing a mesoscopic finite element model of a double-layer reinforced concrete structure according to claim 5, characterized in that: After completing the adjustment of the node coordinates, the corresponding element information is generated. The model uses eight-node hexahedral elements, and each element corresponds to eight nodes. i, j, and k represent the indices 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 element, and the matrix Element is used to store the element node information. The calculation formula for n1 is as follows: n1 = (i + 1) + j·Nx + k·Nx·Ny, Calculation formula for element node 1: Element1 = (i + 1) + j·(Nx + 1) + k·(Nx + 1)·(Ny + 1), Calculation formula for element node 2: Element2 = (i + 2) + j·(Nx + 1) + k·(Nx + 1)·(Ny + 1), Calculation formula for element node 3: Element3 = (i + 2) + (j + 1)·(Nx + 1) + k·(Nx + 1)·(Ny + 1), Calculation formula for element node 4: Element4 = (i + 1) + (j + 1)·(Nx + 1) + k·(Nx + 1)·(Ny + 1), Calculation formula for element node 5: Element5 = (i + 1) + j·(Nx + 1) + (k + 1)·(Nx + 1)·(Ny + 1), Calculation formula for element node 6: Element6 = (i + 2) + j·(Nx + 1) + (k + 1)·(Nx + 1)·(Ny + 1), Calculation formula for element node 7: Element7 = (i + 2) + (j + 1)·(Nx + 1) + (k + 1)·(Nx + 1)·(Ny + 1), Calculation formula for element node 8: Element8 = (i + 1) + (j + 1)·(Nx + 1) + (k + 1)·(Nx + 1)·(Ny + 1), After completing the calculation of the element nodes, the element information of each unit is combined into the overall element matrix Element, and its global index number n1 of the element is added to complete the modeling of the element model.

7. The method for establishing a mesoscopic finite element model of a double-layer steel bar ultra-high performance concrete structure according to claim 6, characterized in that: In S3, when placing the aggregates, the aggregate radius is calculated first, and the specific steps are as follows: The initial value of Asize is 0, and the rand function rand(1) is used, which represents generating 1 random number between 0 and 1. When Asize < (Vpersen·3 / 4), the aggregate radius calculation formula is as follows: R = 2 / 3·Maxsize - a + rand(1)·1 / 3·Maxsize - a, When Vpersen·3 / 4 < Asize ≤ Vpersen·11 / 12, the aggregate radius calculation formula is as follows: R = 1 / 2·Maxsize - a + rand(1)·1 / 6·Maxsize - a, When Vpersen·11 / 12 < Asize, the aggregate radius calculation formula is as follows: R = Minsize - a + rand(1)·Minsize - a, Where Maxsize-a is the maximum radius of aggregate, Minsize-a is the minimum radius of aggregate, and Vpersen is the volume content of aggregate.

8. The method for establishing a microscopic finite element model of a double-layer reinforced ultra-high performance concrete structure according to claim 7, characterized in that: After the aggregate radius calculation is completed, the coordinates of the aggregate sphere center (Cx, Cy, Cz) are placed. The coordinates are randomly generated within the limit to ensure that the sphere does not exceed the boundary. 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 the center coordinates of the aggregates are placed, it is necessary to determine whether there is any intrusion or overlap between the aggregates. Each time an aggregate is randomly generated, it is necessary to determine whether there is any intrusion or overlap with all the previously generated aggregates. The correspondence between the center coordinate distance and the sum of the radius is used to determine whether there is any intrusion or overlap between the aggregates. If there is any overlap, it is necessary to regenerate and determine again. After the aggregate sphere center coordinates are placed, the aggregates need to be processed at the boundary to ensure that all aggregates meet the boundary conditions. A1 is taken as the minimum value of the aggregate sphere center coordinates (Cx, Cy, Cz). The judgment conditions are as follows: R <A1, (Cx+R)≤Lx, (Cy+R)≤Ly, (Cz+R)≤Lz, The volume formula of aggregate is: Vtotal=4 / 3·π·R 3 。 9. The method for establishing a microscopic finite element model of a double-layer reinforced ultra-high performance concrete structure according to claim 8, characterized in that: In S4, when determining the intrusion of aggregate and steel bars, the center coordinates of the spherical aggregate (PointA, Cx, Cy, Cz) are extracted, and the starting coordinates of the first steel bar (Point-G1, X G 1,Y G 1,Z G 1) and the end point coordinates Point-G2(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=|PointA-Point B|, 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.

10. The method for establishing a microscopic finite element model of a double-layer reinforced ultra-high performance concrete structure according to claim 9, characterized in that: Dist judgment conditions are as follows: When Dist<(g-dia+0.01), it is considered that the aggregate and the steel bar have invaded, and a secondary judgment is made. When Dist<(g-dia+0.01)·4 / 7, 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: NewR=(g-dia+R) / 2, NewPointA=PointA+direction·0.4·(g-dia+R), When (g-dia+0.01)·4 / 7≤dist<(g-dia+0.01), a new radius is generated: NewR=0.64·(g-dia+R), After the aggregate has been processed by translation, it is necessary to re-judge whether it has invaded other steel bars until there is no invasion and intersection with all steel bars. The next aggregate information is extracted and the judgment is continued until the intrusion judgment and processing of all aggregates and all steel bars are completed.

Images