A topological optimization method for concrete components based on the first principal tensile stress

Abstract

The present invention relates to a topological optimization method for concrete components based on the first principal tensile stress, and the steps are as follows: S1: define the initial design domain; S2: set the initial optimization parameters, volume rate constraint, evolution rate, maximum resurrection rate, sensitivity filtering radius and design variable filtering radius; S3: perform loop iteration operation: S3.1: calculate the target volume; S3.2: finite element analysis; S3.3: calculate the sensitivity number corresponding to the first principal tensile stress and the sensitivity number of the second principal stress; S3.4: sensitivity correction; S3.5: sensitivity filtering and averaging operation; S3.6: design variable update operation; S3.7: design variable filtering; S3.8: design variable post-processing; S4: determine whether the volume constraint and convergence criterion are met; if so, end; if not, return to S3. The present invention can well reduce the peak level of the first principal tensile stress so that the concrete component has the maximum bearing capacity under a given volume.

Description

Technical Field

[0001] The invention relates to the field of topological optimization of concrete components, and in particular to a topological optimization method of concrete components based on the first principal tensile stress. Background Art

[0002] For example, Chinese invention patent publication number CN 115525991 A discloses a three-main-beam UHPC prestressed cap beam and its topology optimization method, which includes the following steps: the prefabricated cap beam is designed using a topology optimization method, equipped with special-shaped steel bars, and cast using ultra-high performance concrete; according to given load conditions, performance indicators and constraints, the material distribution is topologically optimized in a given design area. Its optimization method does not involve reducing the peak level of the first principal tensile stress. Summary of the invention

[0003] The purpose of the present invention is to provide a topological optimization method for concrete components based on the first principal tensile stress, which can well reduce the peak level of the first principal tensile stress. When the target volume is reached, the stress peak is significantly reduced, the stress concentration is alleviated, the first principal stress tension and compression area in the structure is more obvious, and the overall structural stress is gradually uniform, so that the concrete component has the maximum bearing capacity under a given volume.

[0004] In order to achieve the above purpose, a topology optimization method of concrete components based on the first principal tensile stress is adopted. The specific steps are as follows:

[0005] S1: Define the initial design domain, divide the finite element mesh, and define the boundary conditions;

[0006] S2: Set the initial optimization parameters, volume rate constraint, evolution rate, maximum resurrection rate, sensitivity filter radius and design variable filter radius;

[0007] S3: Execute loop iteration operation. The process is as follows:

[0008] S3.1: Calculate the target volume of this iteration;

[0009] S3.2: Perform finite element analysis on the current topology;

[0010] S3.3: Calculate the sensitivity number corresponding to the first principal tensile stress and the sensitivity number of the second principal stress of the unit;

[0011] S3.4: Correct the sensitivity of the second principal stress according to the biaxial strength and obtain the final corrected sensitivity number;

[0012] S3.5: Filter and average the sensitivity data;

[0013] S3.6: Update the design variables according to the sensitivity;

[0014] S3.7: Filter the design variables;

[0015] S3.8: Post-process the design variables according to the design variables. After filtering the design variables, the values ​​of the design variables are any value from 0 to 1. And because the final topological solution is a discrete structure of either 0 or 1, it is necessary to update the design variables according to their values, i.e., post-processing operation.

[0016] S4: Determine whether the volume constraint and convergence criteria are met; if so, end the design; if not, return to step S3.

[0017] As a further improvement of the present invention, in S3.1, under the linear elasticity of the material, the optimization goal is to minimize the first principal tensile stress of the structural peak value, and the mathematical model for optimization based on volume constraints is described as follows:

[0018] Find: x = [x1, x2, ..., x i ,…,x n ](x i =0or1);

[0019] Minimize:

[0020] Obey:KU=F;

[0021] V(x)=∑x i v i =V req (1);

[0022] In the formula is the maximum first principal tensile stress in the structure, σ i,1 is the first principal stress at the center point of the i-th unit; when the first principal stress of the unit is tensile, the stress value is positive, and it is taken as itself when calculating the objective function; when the first principal stress of the unit is compressive, the stress value is negative, and it is taken as zero when calculating the objective function; V, v i 、V req are the total volume of the structure, the volume of the i-th unit and the target volume respectively; x is the design variable vector, that is, to find a set of solution vectors that minimizes the objective function f(x); i is the design variable value of the i-th unit, which can be 0 or 1. When it is 0, the unit is deleted, and when it is 1, the unit is retained; K and U are the overall stiffness matrix and displacement matrix of the structure respectively, and F is the external load vector.

[0023] As a further improvement of the present invention, in S3.2, it is assumed that the material constitutive and design variables x i The first principal stress at the center of the unit is calculated as follows:

[0024]

[0025] Where σ x , σ y are the principal stresses in the x and y directions, σ xy is the shear stress; A T Represents the matrix corresponding to the first principal stress; in the code, the principal stress of each unit is obtained by performing eigenvalue decomposition on the unit stress matrix after finite element analysis. When taking the first principal stress value, A T represents the eigenvector corresponding to the first principal stress as the eigenvalue;

[0026] The condensation function is used to approximate the maximum value function, and the p-norm in the condensation function is used. Then the objective function is approximately:

[0027] Minimize: In the formula, H(σ i,1 ) is the Heaviside function of the first principal stress of the i-th element, and its general expression is:

[0028]

[0029] In the formula, when the independent variable is less than zero, the function value is 0; otherwise, the function value is 1. When calculating formula (3), when the first principal stress is compressive stress, the unit takes zero when summing the objective function; when the first principal stress is tensile stress, the unit takes the pth power of the first principal stress of the unit when summing the objective function.

[0030] As a further improvement of the present invention, in S3.3, the sensitivity is derived using the adjoint method. First, a linear interpolation model is used to establish the relationship between the design variables and the unit material properties:

[0031] D i =x i D0 (5);

[0032] Where D i represents the actual elastic matrix of the i-th unit during finite element calculation; D0 represents the elastic matrix of the solid material;

[0033] Let the effective stress vector of the i-th element be:

[0034] σ i =D i B i u i / x i =D0B i u i (6);

[0035] Where σ i and ui are the stress vector and node displacement vector of the i-th element, B i represents the geometry matrix of the i-th unit;

[0036] According to equations (2) and (6), the partial derivative of the first principal stress with respect to the design variable is:

[0037]

[0038] Then the partial derivative of the objective function with respect to the design variables is:

[0039]

[0040] Substituting formula (7) into formula (8), we can further obtain:

[0041]

[0042] Where, L j is the index matrix from the global displacement vector to the local displacement vector of the j-th unit, satisfying u j =L j U;

[0043] The structure must satisfy the equilibrium condition, that is, KU = F. The adjoint method is used to obtain the partial derivatives of the design variables on both sides of the equation:

[0044]

[0045] In the formula, since the external load F is a known constant, the partial derivative of the right side of the equation with respect to the design variable is zero. Substituting the above formula into formula (9) yields:

[0046]

[0047] In order to simplify the above formula, a matrix factor λ is introduced to satisfy:

[0048]

[0049] Substituting equation (12) into equation (11), the derivative of the objective function with respect to the design variable can be obtained:

[0050]

[0051] The partial derivative of the overall stiffness matrix with respect to the design variables is:

[0052]

[0053] Where k j represents the stiffness matrix of the j-th element, Represents the stiffness matrix of the solid element. The specific expressions of the two are and

[0054] Then the sensitivity of the objective function to the design variables can be changed to:

[0055]

[0056] Multiply the formula (15) by a factor with opposite sign "-x i ", the original sensitivity expression of each unit is obtained:

[0057]

[0058] On the basis of the first principal tensile stress, the second principal stress is considered. Similar to the first principal tensile stress, the aggregation function of the second principal stress is obtained, and the sensitivity number is derived. However, it should be noted that the tension and compression state of the unit is not distinguished when considering the second principal stress.

[0059] As a further improvement of the present invention, in S3.4, according to the biaxial strength theory of concrete materials, the second sensitivity is reduced and summed with the sensitivity calculated by formula (16) to obtain the final sensitivity value, and the optimization is carried out based on this. The specific processing is as follows:

[0060] Similar to equation (16), the derivative of the second principal stress with respect to the design variable is:

[0061]

[0062] In the formula is the aggregation function of the second principal stress. Note that since the second principal stress is calculated without distinguishing between tension and compression, the second term on the far right of the equation is is a function of the actual second principal stress of each element. In addition, λ2 differs from λ in two ways: first, there is no Heaviside function; second, the eigenvector at this time is calculated as the eigenvalue corresponding to the second principal stress;

[0063] According to the biaxial strength curve of concrete, the second sensitivity is reduced. The biaxial strength curve function is as follows:

[0064] Lala:

[0065] Pressure:

[0066] Tension and compression:

[0067] Compression and tension:

[0068] In the formula, α s is the shear yield parameter, and its calculation formula is r is the biaxial compressive strength enhancement factor; ν is Poisson's ratio; f1 and f2 are the first and second principal stresses of the unit, respectively; f t and f c is the tensile strength and compressive strength of concrete;

[0069] First, the sensitivity of the second principal stress is calculated according to formula (17), and then divided by 10; then, it is reduced according to the actual bidirectional stress state of the unit. The specific reduction is as follows:

[0070] If the obtained bidirectional stress state of the unit is within the envelope diagram, the actual second principal stress of each unit is used to calculate the first principal stress on the envelope diagram by the above formula (18), and then the absolute value of the actual first principal stress of the unit is subtracted from the calculated first principal stress, and the minimum value is taken as the final reduction coefficient t of each unit; otherwise, no reduction is performed;

[0071] The final sensitivity of the processed unit is:

[0072]

[0073] Where t i It represents the minimum absolute value of the difference between the actual first principal stress value of the i-th element and the biaxial strength envelope of the concrete.

[0074] As a further improvement of the present invention, in S3.5, a filtering operation is performed on the unit sensitivity:

[0075]

[0076] In the formula is the linear weight factor, calculated as follows:

[0077]

[0078] Where r sen is the sensitivity filter radius, which is usually 2 to 3 times the length of the unit side to ensure that the deleted unit can be restored; Δ(i,j) is the distance between the center points of the i-th unit and the j-th unit;

[0079] At the same time, the filtered sensitivity is also averaged with its historical information, that is, the filtered sensitivity is averaged with the historical information of the previous two generations to further filter the sensitivity:

[0080]

[0081] Where (l) represents the number of optimization iteration cycles.

[0082] As a further improvement of the present invention, in S3.6, based on the filtered and averaged sensitivity values, the design variable update operation is performed as follows:

[0083]

[0084] In the formula and The upper and lower sensitivity thresholds are set for the desired number of cell deletions and restorations corresponding to the contemporary target volume.

[0085] As a further improvement of the present invention, in S3.7, a filtering operation is also performed on the design variables according to the following formula:

[0086]

[0087] In the formula is the linear weight factor, calculated as follows:

[0088]

[0089] Where r den Filter radius for the design variable.

[0090] As a further improvement of the present invention, in S3.8, the design variable update operation is performed again to ensure that the topological structure is a clear discrete structure:

[0091]

[0092] In the formula and Upper and lower thresholds for design variables set for the desired number of element deletions and restorations corresponding to the contemporary target volume.

[0093] As a further improvement of the present invention, when the volume rate remains unchanged, the iteration is continued according to the following formula:

[0094]

[0095] Where C i represents the i-th generation objective function value, that is, the p-norm of the first principal tensile stress of the i-th generation; k represents the current iteration number; N represents a constant; τ represents the allowable deviation;

[0096] After reaching the target volume, run 2N generations and compare the proportion of the difference between the sum of the objective functions of the last N generations and the sum of the objective functions of the first N generations. If it is less than the set limit, it means convergence. Otherwise, run one more generation and use the formula for another comparison. If the formula is still not satisfied after 2N generations, stop the iteration and take the topological solution corresponding to the minimum value of the objective function of these 4N generations as the optimal topological solution for output. BRIEF DESCRIPTION OF THE DRAWINGS

[0097] Figure 1 is the design domain of the L-shaped structure example.

[0098] Figure 2 The first principal tensile stress peak and volume diagram when p=4 varies with the number of iterations.

[0099] Figure 3 This is the topological diagram of some generations when p=4. 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] In the description of the present invention, it should be noted that the terms "center", "upper", "lower", "left", "right", "vertical", "horizontal", "inside", "outside" and the like indicate positions or positional relationships based on the positions or positional relationships shown in the drawings, and are only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the device or element referred to must have a specific position, be constructed and operated in a specific position, and therefore cannot be understood as limiting the present invention; the terms "first", "second", and "third" are only used for descriptive purposes, and cannot be understood as indicating or implying relative importance; in addition, unless otherwise clearly specified and limited, the terms "installed", "connected", and "connected" should be understood in a broad sense, for example, it can be a fixed connection, a detachable connection, or an integral connection; it can be a mechanical connection or an electrical connection; it can be a direct connection, or it can be an indirect connection through an intermediate medium, or it can be a connection between the two elements. For those of ordinary skill in the art, the specific meanings of the above terms in the present invention can be understood according to specific circumstances.

[0102] Example 1

[0103] A topology optimization method based on the first principal tensile stress

[0104] 1.1 Optimization goals

[0105] Under the linear elasticity of the material, the optimization goal is to minimize the peak first principal tensile stress of the structure. The mathematical model for optimization based on volume constraints can be described as follows:

[0106] Find: x = [x1, x2, ..., x i ,…,x n ](x i =0or1);

[0107] Minimize:

[0108] Obey:KU=F;

[0109] V(x)=∑x i v i =V req (1);

[0110] In the formula is the maximum first principal tensile stress in the structure, σ i,1 is the first principal stress at the center point of the i-th unit; in this embodiment, when the first principal stress of the unit is tensile, the stress value is positive, and it is taken as itself when calculating the objective function; when the first principal stress of the unit is compressive, the stress value is negative, and it is taken as zero when calculating the objective function; V, v i 、V req are the total volume of the structure, the volume of the i-th unit and the target volume respectively; x is the design variable vector, that is, to find a set of solution vectors that minimizes the objective function f(x); i is the design variable value of the i-th unit, which can be 0 or 1. When it is 0, the unit is deleted, and when it is 1, the unit is retained; K and U are the overall stiffness matrix and displacement matrix of the structure respectively, and F is the external load vector.

[0111] Taking 2D plane stress as an example, assuming that the material constitutive and design variables x i Regardless, the first principal stress at the center of the unit is calculated as follows:

[0112]

[0113] Where σ x , σ y are the principal stresses in the x and y directions, σ xy is the shear stress; A T represents the matrix corresponding to the first principal stress. For a two-dimensional plane, the size of the matrix is ​​1×3. Assume that A T =(abc), then the matrix satisfies a+b=1; in the code, the principal stress of each unit is obtained by performing eigenvalue decomposition on the unit stress matrix after finite element analysis. When taking the first principal stress value, A TRepresents the eigenvector corresponding to the first principal stress as the eigenvalue.

[0114] Since stress is a local quantity and the objective function is discontinuous and non-differentiable, the condensation function is used to approximate the maximum value function. This embodiment uses the p-norm, and the objective function is approximately:

[0115] Minimize: In the formula, H(σ i,1 ) is the Heaviside function of the first principal stress of the i-th element, and its general expression is:

[0116]

[0117] In the formula, when the independent variable is less than zero, the function value is 0; otherwise, the function value is 1. When calculating formula (3), when the first principal stress is compressive stress, the unit takes zero when summing the objective function; when the first principal stress is tensile stress, the unit takes the pth power of the first principal stress of the unit when summing the objective function.

[0118] 1.2 Optimizing sensitivity

[0119] The sensitivity is derived using the adjoint method. First, a linear interpolation model is used to establish the relationship between the design variables and the element material properties:

[0120] D i =x i D0 (5);

[0121] Where D i It represents the actual elastic matrix of the i-th unit during finite element calculation; D0 represents the elastic matrix of the solid material.

[0122] Let the effective stress vector of the i-th element be:

[0123] σ i =D i B i u i / x i =D0B i u i (6);

[0124] Where σ i and u i are the stress vector and node displacement vector of the i-th element, B i Represents the geometric matrix of the i-th element. In practical applications, in order to avoid mesh re-division and singularity of overall stiffness, extremely low elastic modulus is assigned to the empty element.

[0125] According to equations (2) and (6), the partial derivative of the first principal stress with respect to the design variable is:

[0126]

[0127] Then the partial derivative of the objective function with respect to the design variables is:

[0128]

[0129] Substituting formula (7) into formula (8), we can further obtain:

[0130]

[0131] Where, L j is the index matrix from the global displacement vector to the local displacement vector of the j-th unit, satisfying u j =L j U.

[0132] Since the structure must satisfy the equilibrium condition, that is, KU = F, in each finite element analysis and iteration process, the adjoint method is used to obtain the partial derivative of the design variables on both sides of the equation:

[0133]

[0134] In the formula, since the external load F is a known constant, the partial derivative of the right side of the equation with respect to the design variable is zero. Substituting the above formula into formula (9) yields:

[0135]

[0136] In order to simplify the above formula, a matrix factor λ is introduced to satisfy:

[0137]

[0138] Substituting equation (12) into equation (11), the derivative of the objective function with respect to the design variable can be obtained:

[0139]

[0140] The partial derivative of the overall stiffness matrix with respect to the design variables is:

[0141]

[0142] Where k j represents the stiffness matrix of the j-th element, Represents the stiffness matrix of the solid element. The specific expressions of the two are and

[0143] Then the sensitivity of the objective function to the design variables can be changed to:

[0144]

[0145] Since the traditional BESO algorithm deletes units with small sensitivity first, the greater the sensitivity calculated in the above formula, the smaller the contribution to the objective function. Therefore, in order to maintain the same increase and decrease monotonicity of sensitivity as the traditional BESO algorithm, a negative sign factor "-x" is multiplied on formula (15). i ", the original sensitivity expression of each unit is obtained:

[0146]

[0147] Since the stress value of the compressive unit is not considered when calculating the objective function, that is, the addition or deletion of units whose first principal stress is compressive stress has no effect on the change of the objective function. Therefore, the sensitivity of the objective function to the compressive unit is zero. However, the topological solution obtained by the sensitivity obtained in this way is deformed and has poor manufacturability, which is obviously unreasonable. Therefore, the initial sensitivity needs to be processed. The processing method of this paper is to consider the second principal stress on the basis of the first principal tensile stress. Similar to the first principal tensile stress, the aggregation function of the second principal stress is obtained, and the sensitivity number is derived. However, it should be noted that the tension and compression state of the unit is not distinguished when considering the second principal stress. Then, according to the biaxial strength theory of concrete materials, the second sensitivity is reduced and summed with the sensitivity calculated by formula (16) to obtain the final sensitivity value, and the optimization is carried out based on this. The specific processing is as follows:

[0148] Similar to equation (16), the derivative of the second principal stress with respect to the design variable is:

[0149]

[0150] In the formula is the aggregation function of the second principal stress. Note that since the second principal stress is calculated without distinguishing between tension and compression, the second term on the far right of the equation (i.e. ) is a function of the actual second principal stress of each element. In addition, λ2 is different from λ in two aspects: one is that there is no Heaviside function; the other is that the eigenvectors in this case are calculated as the eigenvalues ​​corresponding to the second principal stress.

[0151] According to the concrete biaxial strength curve, the second sensitivity is reduced. The biaxial strength curve function is as follows:

[0152] Lala:

[0153] Pressure:

[0154] Tension and compression:

[0155] Compression and tension:

[0156] In the formula, αs is the shear yield parameter, and its calculation formula is Where r is the biaxial compressive strength enhancement coefficient, the value range is 1.15-1.30, and the r value in this embodiment is 1.20; ν is Poisson's ratio; f1 and f2 are the first principal stress and the second principal stress of the unit respectively; f t and f c The tensile strength and compressive strength of concrete.

[0157] First, the sensitivity of the second principal stress is calculated according to formula (17), and then divided by 10; then, it is reduced according to the actual bidirectional stress state of the unit. The specific reduction is as follows:

[0158] If the obtained bidirectional stress state of the unit is within the envelope diagram, the actual second principal stress of each unit is calculated by the above formula to obtain the first principal stress on its envelope diagram, and then the absolute values ​​of the actual first principal stress of the unit and the calculated first principal stress are subtracted, and the minimum value is taken as the final reduction coefficient t of each unit; otherwise, no reduction is performed.

[0159] The final sensitivity of the processed unit is:

[0160]

[0161] Where t i It represents the minimum absolute value of the difference between the actual first principal stress value of the i-th element and the biaxial strength envelope of the concrete.

[0162] 3. Variable update mechanism and post-processing

[0163] In order to avoid common numerical problems such as mesh dependency and checkerboard phenomena in topology optimization, it is necessary to filter the element sensitivity:

[0164]

[0165] In the formula is the linear weight factor, calculated as follows:

[0166]

[0167] Where r sen is the sensitivity filter radius, which is usually 2 to 3 times the length of the unit side to ensure that the deleted units can be restored; Δ(i,j) is the distance between the center points of the i-th unit and the j-th unit.

[0168] At the same time, the filtered sensitivity is also averaged with its historical information to further stabilize the optimization process. That is, the filtered sensitivity is averaged with its previous two generations of historical information to further filter the sensitivity:

[0169]

[0170] Where (l) represents the number of optimization iteration cycles.

[0171] Then, based on the filtered and averaged sensitivity values, the design variable update operation is performed as follows:

[0172]

[0173] In the formula and The upper and lower sensitivity thresholds are set according to the number of required unit deletions and restorations corresponding to the contemporary target volume. To improve the stability and robustness of the algorithm, a filtering operation is also performed on the design variables:

[0174]

[0175] In the formula is the linear weight factor, calculated as follows:

[0176]

[0177] Where r den Filter radius for the design variable.

[0178] Since the design variables may contain any value between 0 and 1 at this time, it is necessary to perform another design variable update operation to ensure that the topological structure is a clear discrete structure:

[0179]

[0180] In the formula and Similar to the upper and lower thresholds determined previously based on sensitivity, here are the upper and lower thresholds of the design variable set based on the required number of unit deletions and restorations corresponding to the contemporary target volume.

[0181] When the target volume is reached, in order to avoid falling into a local solution and to find the global optimal solution as much as possible, the iteration is continued according to the following formula while the volume rate remains unchanged:

[0182]

[0183] Where C iRepresents the value of the objective function of the i-th generation. In this embodiment, it represents the p-norm of the first principal tensile stress of the i-th generation; k represents the current number of iterations; N represents a constant. In this embodiment, N=5; τ represents the allowable deviation. In this embodiment, τ=0.001. The overall meaning of the formula is that after reaching the target volume, run ten generations forward and compare the proportion of the sum of the objective functions of the last five generations and the sum of the objective functions of the first five generations in the sum of the objective functions of the last five generations. If it is less than the set limit, it means that convergence has been achieved; otherwise, run one more generation and use the formula for another comparison. If the formula is still not satisfied after ten generations of execution, stop the iteration and take the topological solution corresponding to the minimum value of the objective function of these twenty generations as the optimal topological solution for output.

[0184] 1.4 Iteration Flowchart

[0185] The whole process of the optimization problem based on the first principal tensile stress is as follows:

[0186] (1) Define the initial design domain, divide the finite element mesh, and define the boundary conditions;

[0187] (2) Setting initial optimization parameters, volume rate constraint, evolution rate, maximum resurrection rate, sensitivity filter radius, and design variable filter radius;

[0188] (3) Start loop iteration:

[0189] ① Calculate the target volume of this iteration;

[0190] ② Conduct finite element analysis on the current topological structure;

[0191] ③ Calculate the sensitivity number corresponding to the first principal tensile stress and the sensitivity number of the second principal stress of the unit;

[0192] ④ Correct the sensitivity of the second principal stress according to the biaxial strength and obtain the final corrected sensitivity number;

[0193] ⑤ Filter and average the sensitivity numbers;

[0194] ⑥ Update the design variables according to sensitivity;

[0195] ⑦Filter the design variables;

[0196] ⑧ Post-process the design variables according to the design variables;

[0197] (4) Determine whether the volume constraints and convergence criteria are met; if yes, end the design; if no, return to step (3).

[0198] 2 Numerical examples

[0199] 2.1 L-shaped beam example

[0200] For a concrete L-shaped beam, the design domain and other basic design conditions are as follows: Figure 1 As shown. Preset its Young's modulus E c =3.25×10 4 N / mm 2 , Poisson's ratio is 0.2. Using the plane stress assumption, the side length h e =1mm square bilinear elements are used to discretize the structure. The beam is subjected to a load of F = 4N, which is evenly distributed to eight nodes in two horizontal rows to avoid stress concentration. The optimized volume rate constraint is set to 40%, and the sensitivity filter radius r sen and design variable filter radius r den Take 3 and 5 times the unit side length respectively. In each iteration, the evolution rate c er 2%, the maximum percentage of units that can be revived 0.5%. Using C40 concrete, the tensile strength f t and compressive strength f c 1.71 and 19.1 N / mm respectively 2 The iteration base n after convergence is 10. The maximum first principal tensile stress in the structure at the beginning is 2.28N / mm 2 When the p value is 4, the corresponding objective function value is 4.25N / mm 2 .

[0201] like Figure 2 As shown, Figure 2 The first principal tensile stress peak and volume diagram when p=4 varies with the number of iterations.

[0202] Depend on Figure 2 It can be seen that before reaching the target volume, as the number of iterations increases, the volume of the remaining material in the structure continues to decrease, and the peak value of the first principal tensile stress fluctuates, especially around the 36th generation, when the first principal tensile stress oscillates more; after reaching the volume constraint, the volume rate in the structure remains unchanged, only the material distribution changes, and the peak value of the first principal tensile stress changes more slowly than before, and shows an upward trend. Therefore, the oscillation of the peak value of the first principal tensile stress indicates that a slight change in the material distribution may cause a drastic change in the local quantity including it.

[0203] like Figure 3 As shown, Figure 3 This is the topological diagram of some generations when p=4.

[0204] Figure 3 The figure shows some generation structure diagrams and the first principal stress cloud diagram. It should be noted that a positive principal stress indicates that the first principal stress of the unit is tensile stress. On the contrary, a negative value indicates compressive stress.

[0205] from Figure 3From the first iteration and the generation in which the tensile stress peak is minimized, it can be seen that the maximum first principal tensile stress of the initial generation structure is 2.28, and it reaches 2.13 in the final generation. When the volume is reduced by 60%, the stress peak is also decreasing, but the p-norm function value is increasing. The reason for this is that the stress at the corner of the initial structure is relatively concentrated, and the function value mainly depends on the unit of the stress concentration part; when the tensile stress peak is minimized, according to the principal stress cloud map, many units are close to the tensile stress peak, and the contribution of the maximum tensile stress unit to the objective function gradually decreases, while the number of units close to the maximum tensile stress increases. Therefore, the objective function value is large.

[0206] From the iteration process, it can be seen that the second generation structure deleted the stress concentration units at the corners and some units on the left side of the corners, causing the peak value of the first principal tensile stress in the structure to drop sharply; in the subsequent iterations, the structure continued to find its shape, and the peak value of the first principal tensile stress fluctuated up and down. The initial configuration has already appeared in the 22nd generation. The 36th generation is the maximum value of the peak value of the principal tensile stress in the entire iteration process; the 47th generation is the iteration diagram when the target volume is reached for the first time, and this is the minimum value of the peak value of the first principal tensile stress after reaching the target volume. In addition, from the optimal topology generation, it can be clearly seen that the tension and compression state of the first principal stress of the units in the structure, that is, the first principal stress of the units of the leftmost, lowermost and connecting loading point members of the L-shaped structure is less than zero, which is a compressive part. On the contrary, the first principal stress of the units corresponding to the inner members of the L-shaped structure is greater than zero, which is a tensile unit.

[0207] 3. Conclusion

[0208] (1) The topological optimization method based on the first principal tensile stress proposed in this paper can effectively reduce the peak level of the first principal tensile stress. When the target volume is reached, the stress peak is significantly reduced, the stress concentration is alleviated, the first principal stress tension and compression area in the structure is more obvious, and the overall structural stress is gradually uniform.

[0209] (2) As the volume of the structure continues to decrease, the overall stress level of the structure becomes more and more uniform, and the p-norm value gradually becomes less dependent on the stress peak, resulting in an increasing number of units in the structure that are close to the stress peak.

[0210] The above contents are further detailed descriptions of the present invention in combination with specific preferred embodiments, and it cannot be determined that the specific implementation of the present invention is limited to these descriptions. For those skilled in the art of the present invention, several equivalent substitutions or obvious variations can be made without departing from the concept of the present invention, and the performance or use is the same, which should be regarded as belonging to the protection scope of the present invention.

Claims

1. A topological optimization method for concrete components based on the first principal tensile stress, characterized in that The specific steps are as follows: S1: Define the initial design domain, divide the finite element mesh, and define the boundary conditions; S2: Set the initial optimization parameters, volume rate constraint, evolution rate, maximum resurrection rate, sensitivity filter radius and design variable filter radius; S3: Execute loop iteration operation. The process is as follows: S3.1: Calculate the target volume of this iteration; S3.2: Perform finite element analysis on the current topology; S3.3: Calculate the sensitivity number corresponding to the first principal tensile stress and the sensitivity number of the second principal stress of the unit; S3.4: Correct the sensitivity of the second principal stress according to the biaxial strength and obtain the final corrected sensitivity number; S3.5: Filter and average the sensitivity data; S3.6: Update the design variables according to the sensitivity; S3.7: Filter the design variables; S3.8: Post-process the design variables according to the design variables; S4: Determine whether the volume constraint and convergence criteria are met; if so, end the design; if not, return to step S3.

2. The method for topological optimization of concrete components based on the first principal tensile stress according to claim 1, characterized in that In S3.1, under the linear elasticity of the material, the optimization objective is to minimize the peak first principal tensile stress of the structure. The mathematical model for optimization based on volume constraints is described as follows: Search for: ; Minimize: ; obey: ; （1）； In the formula is the maximum first principal tensile stress in the structure, σ i,1 For the i The first principal stress at the center point of the unit; when the first principal stress of the unit is tensile, the stress value is positive and is taken as itself when calculating the objective function; When the first principal stress of the element is under compression, the stress value is negative and is taken as zero when calculating the objective function; V、v i 、V req They are the total volume of the structure, i The volume of the unit and the target volume; x is the design variable vector, that is, to find a set of solution vectors so that the objective function f ( x ) takes the minimum value; x i For the i The design variable value of the unit can be 0 or 1. When it is 0, the unit is deleted, and when it is 1, the unit is retained; K and U are the overall stiffness matrix and displacement matrix of the structure respectively, and F is the external load vector.

3. The method for topological optimization of concrete components based on the first principal tensile stress according to claim 2, characterized in that: In S3.2, it is assumed that the material constitutive and design variables x i The first principal stress at the center of the unit is calculated as follows: （2）； In the formula σ x , σ y They are x , y The principal stress in the direction σ xy is the shear stress; A T Represents the matrix corresponding to the first principal stress; in the code, the principal stress of each unit is obtained by performing eigenvalue decomposition on the unit stress matrix after finite element analysis. When taking the first principal stress value, A T represents the eigenvector corresponding to the first principal stress as the eigenvalue; The condensation function is used to approximate the maximum value function, and the condensation function is used to approximate the maximum value function. p norm, then the objective function is approximately: Minimize: (3) In the formula H ( σ i,1 ) is the i The general expression of the Heaviside function of the first principal stress of the unit is: （4）； In the formula, when the independent variable is less than zero, the function value is 0; otherwise, the function value is 1. When calculating formula (3), when the first principal stress is compressive stress, the unit takes zero when summing the objective function; when the first principal stress is tensile stress, the unit takes the first principal stress of the unit when summing the objective function. p Power.

4. The method for topological optimization of concrete components based on the first principal tensile stress according to claim 3, characterized in that In S3.3, the adjoint method is used to derive the sensitivity. First, a linear interpolation model is used to establish the relationship between the design variables and the unit material properties: （5）； Where D i Indicates i The actual elastic matrix of the unit D during finite element calculation; 0 represents the elastic matrix of the solid material; Order i The effective stress vector of element number is: （6）； Where σ i and u i Respectively i The stress vector and node displacement vector of the unit, B i Indicates i The geometric matrix of the unit number; According to equations (2) and (6), the partial derivative of the first principal stress with respect to the design variable is: （7）； Then the partial derivative of the objective function with respect to the design variables is: （8）； Substituting formula (7) into formula (8), we can further obtain: （9）； Where, L j For the j The index matrix of the number unit from the global displacement vector index to the local displacement vector satisfies ; The structure must satisfy the equilibrium condition, that is, KU=F. The adjoint method is used to obtain the partial derivatives of the design variables on both sides of the equation: （10）； In the formula, since the external load F is a known constant, the partial derivative of the right side of the equation with respect to the design variable is zero. Substituting the above formula into formula (9) yields: （11）； In order to simplify the above formula, a matrix factor is introduced λ ,satisfy: （12）； Substituting equation (12) into equation (11), the derivative of the objective function with respect to the design variable can be obtained: （13）； The partial derivative of the overall stiffness matrix with respect to the design variables is: （14）； Where k j Indicates The stiffness matrix of the element, Represents the stiffness matrix of the solid element. The specific expressions of the two are and ; Then the sensitivity of the objective function to the design variables can be changed to: （15）； Multiply the formula (15) by a factor with opposite sign" ", the original sensitivity expression of each unit is obtained: （16）； On the basis of the first principal tensile stress, the second principal stress is considered. Similar to the first principal tensile stress, the aggregation function of the second principal stress is obtained, and the sensitivity number is derived. However, it should be noted that the tension and compression state of the unit is not distinguished when considering the second principal stress.

5. The method for topological optimization of concrete components based on the first principal tensile stress according to claim 4, characterized in that In S3.4, according to the biaxial strength theory of concrete materials, the second sensitivity is reduced and summed with the sensitivity calculated by formula (16) to obtain the final sensitivity value, which is then used for optimization. The specific processing is as follows: Similar to equation (16), the derivative of the second principal stress with respect to the design variable is: （17）； In the formula is the aggregation function of the second principal stress. Note that since the second principal stress is calculated without distinguishing between tension and compression, the second term on the far right of the equation is is a function of the actual second principal stress for each element, and, λ 2 and λ There are two differences: first, there is no Heaviside function; second, the eigenvectors are calculated as the eigenvalues ​​corresponding to the second principal stress; According to the biaxial strength curve of concrete, the second sensitivity is reduced. The biaxial strength curve function is as follows: ； ； ； （18）； In the formula, α s is the shear yield parameter, and its calculation formula is , r is the biaxial compressive strength enhancement factor; ν is the Poisson’s ratio; f 1 and f 2 are the first principal stress and the second principal stress of the element respectively; f t and f c is the tensile strength and compressive strength of concrete; First, calculate the sensitivity of the second principal stress according to formula (17), and then divide it by 10; then reduce it according to the actual bidirectional stress state of the unit. The specific reduction is as follows: If the bidirectional stress state of the unit is within the envelope diagram, the actual second principal stress of each unit is calculated by the above formula (18) to calculate the first principal stress on the envelope diagram. Then, the absolute value of the actual first principal stress of the unit is subtracted from the calculated first principal stress, and the minimum value is taken as the final reduction factor of each unit. t ; Otherwise, no reduction will be made; The final sensitivity of the processed unit is: （19）； In the formula t i Indicates i The minimum absolute value of the difference between the actual first principal stress value of element No. and the biaxial strength envelope of concrete.

6. The method for topological optimization of concrete components based on the first principal tensile stress according to claim 5, characterized in that In S3.5, the unit sensitivity is filtered: （20）； In the formula is the linear weight factor, calculated as follows: （21）； In the formula r sen The sensitivity filter radius is set to 2 to 3 times the cell side length to ensure that deleted cells can be restored. For the i Unit No. and j The distance between the center points of the units; At the same time, the filtered sensitivity is also averaged with its historical information, that is, the filtered sensitivity is averaged with the historical information of the previous two generations to further filter the sensitivity: （22）； Where ( l ) represents the number of optimization iteration cycles.

7. The method for topological optimization of concrete components based on the first principal tensile stress according to claim 6, characterized in that In S3.6, based on the filtered and averaged sensitivity values, the design variable update operation is performed as follows: （23）； In the formula and The upper and lower sensitivity thresholds are set for the desired number of cell deletions and restorations corresponding to the contemporary target volume.

8. The method for topological optimization of concrete components based on the first principal tensile stress according to claim 7, characterized in that In S3.7, a filtering operation is also performed on the design variables as follows: （24）； In the formula is the linear weight factor, calculated as follows: （25）； In the formula r den Filter radius for the design variable.

9. The method for topological optimization of concrete components based on the first principal tensile stress according to claim 8, characterized in that In S3.8, the design variable update operation is performed again to ensure that the topological structure is a clear discrete structure: （26）； In the formula and Upper and lower thresholds for design variables set for the desired number of element deletions and restorations corresponding to the contemporary target volume.

10. The method for topological optimization of concrete components based on the first principal tensile stress according to claim 9, characterized in that In S4, when the volume rate remains unchanged, the iteration continues according to the following formula: （27）； In the formula C i Indicates i The objective function value is i The first principal tensile stress of p norm; k Indicates the current iteration number; N represents a constant; τ Indicates the allowable deviation; After reaching the target volume, run 2 N Generation, after comparison N The sum of the objective function and the previous N The difference between the sum of the objective functions is N If the proportion of the sum of the objective functions is less than the set limit, it means that convergence has been achieved; Otherwise, run one more generation and use the formula to make another comparison. If the execution continues for 2 N If the formula is still not satisfied, the iteration is stopped and these 4 N The topological solution corresponding to the minimum value of the objective function is output as the optimal topological solution.

Images