Dam temperature prediction method based on finite element analysis and thermal parameter inversion

Abstract

The present invention relates to a dam temperature prediction method based on finite element analysis and thermal parameter inversion, comprising: establishing a finite element model of each concrete pouring bin; determining thermal parameters and establishing an objective function; performing temperature simulation calculations using the finite element model and calculating the objective function value; performing iterative inversion calculations on the thermal parameters using an artificial bee colony algorithm to determine the quality of the obtained thermal parameter values; repeating the inversion calculations to obtain thermal parameter values ​​that meet the requirements, and comparing the calculated temperature values ​​with the actual concrete temperature values, stopping the iterative inversion after the expected target is achieved; and predicting the temperature of the dam concrete based on the obtained thermal parameter values. The present invention utilizes the ABC algorithm to perform inversion calculations on the thermal parameters, thereby reducing the error between the calculated thermal parameter values ​​and the actual values, improving the accuracy of dam temperature prediction, and facilitating the implementation of refined real-time temperature control of the dam.

Description

Technical Field

[0001] The invention belongs to the field of temperature control of water conservancy and hydropower concrete dams, and particularly relates to a dam temperature prediction method based on finite element analysis and thermal parameter inversion. Background Art

[0002] During the construction of concrete dams, significant temperature stresses are easily generated within the dam body due to factors such as the inherent heat of hydration and external temperature. When the concrete's tensile strength is insufficient to withstand these temperature stresses, thermal cracks are highly likely to form. To accurately and in real time understand the dam's temperature behavior, relying solely on current temperature data from existing instruments is insufficient. Software simulation is required to predict future dam temperature trends and obtain an accurate concrete temperature field for analyzing the dam's stress state. However, the concrete temperature field simulation process is susceptible to numerous uncertainties, such as external temperature, internal cooling water flow, and concrete thermal parameters. While external temperature and cooling water flow can be monitored on-site, thermal parameters must be determined through laboratory testing. However, due to the limitations of laboratory testing, there is a significant discrepancy between the experimental and actual values ​​of thermal parameters. Therefore, it is necessary to use field-measured data from the dam construction site and intelligent optimization analysis methods, comprehensively considering the influence of external temperature and cooling water flow, and using inverse calculations to determine the true thermal parameter values ​​during dam construction.

[0003] The inverse calculation process for dam thermal parameters essentially involves multiple iterations of a temperature field simulation program using an optimization algorithm. Existing optimization algorithms for temperature field simulations have poor local search capabilities and a strong dependence on initial values. For complex problems, they are prone to falling into local optima and cannot guarantee the global optimal solution. Furthermore, they suffer from high computational complexity and long iteration times.

[0004] The Artificial Bee Colony (ABC) algorithm is a biomimetic intelligent computing method that simulates the behavior of bee colonies in their quest for high-quality nectar sources. Proposed by Turkish scholar Karaboga in 2005, it was developed by applying the collective intelligence of bee colonies to solving optimization problems, combining it with artificial intelligence techniques. This algorithm falls under the umbrella of swarming algorithms, similar to the fish swarm algorithm and the ant colony algorithm. The ABC is essentially a generalized neighborhood search algorithm that divides a bee colony into three types: leader bees, follower bees, and scout bees, and two basic behaviors: sharing nectar sources and abandoning them. The location of each nectar source represents a feasible solution to the optimization problem, the nectar volume at the source represents the quality of the corresponding solution, and the speed at which the bees search for the source represents the speed at which the corresponding optimization problem is solved. By switching between leader, follower, and scout bees in different situations and leveraging heuristic strategies, the ABC algorithm not only effectively performs local search but also possesses the ability to find global optimization.

[0005] Compared with other intelligent algorithms, the ABC algorithm has the advantages of fewer control parameters, simple and easy implementation of calculations, fast search speed, and fitness function as the main evolutionary basis. However, ABC still has problems such as slow convergence speed in the later stage and easy to fall into local optimal solutions.

[0006] Therefore, the improved artificial bee colony algorithm is studied and applied to the inverse calculation of the dam thermal parameters to achieve accurate prediction of the dam temperature. Summary of the Invention

[0007] The purpose of the present invention is to address the above-mentioned problems and provide a dam temperature prediction method based on finite element analysis and thermal parameter inversion. In the search process of the ABC algorithm solution space, the guiding role of the current global optimal value is increased to accelerate the convergence speed of the ABC algorithm, and a cross operation is added in the leading bee field search process to improve the global optimization ability; the improved ABC algorithm is used to perform inverse calculation of the thermal parameters to reduce the error between the calculated thermal parameter values ​​and the true values; and the dam temperature is predicted based on the calculated thermal parameter values ​​to facilitate the implementation of refined temperature control of the dam.

[0008] The technical solution of the present invention is a dam temperature prediction method based on finite element analysis and thermal parameter inversion, comprising the following steps:

[0009] Step 1: Establish the finite element model of each concrete pouring bin and set the initial conditions and boundary conditions;

[0010] Step 1-1: In the finite element software, establish a finite element model of the casting silo and mesh the finite element model; at the same time, assign values ​​to the concrete material and unit properties;

[0011] Step 1-2: The initial conditions of the finite element model include the initial temperature of the concrete, which is obtained based on the measured data on site;

[0012] Steps 1-3: Boundary conditions mainly include the external environment temperature and water cooling of concrete;

[0013] Step 2: Determine the thermal parameters and their value ranges used for dam temperature prediction, and establish the objective function of thermal parameters;

[0014] Step 2-1: Determine the inverse thermal parameters and their value ranges. The inverse thermal parameters include the concrete surface heat release coefficient, concrete temperature rise law parameters, and the final adiabatic temperature rise of concrete.

[0015] Step 2-2: Establish the objective function of thermal parameter inversion and transform the concrete thermal parameter inversion problem into an optimization problem;

[0016] Step 3: Randomly generate a set of initial thermal parameter values ​​within the range of thermal parameter values ​​and substitute them into the finite element model for temperature simulation calculation to obtain the temperature values ​​at the nodes in the concrete pouring bin. Substitute them into the objective function together with the actual temperature of the concrete at the nodes to calculate the objective function value;

[0017] Step 3-1: Simulate the concrete pouring sequence, atmospheric temperature, pouring temperature, and hydration heat;

[0018] Step 3-2: Use the "element generation" and "element killing" functions of the finite element software to simulate the concrete pouring in compartments. Before the simulation of the compartment pouring, all dam body elements are "killed". During the simulation, the dam body elements are activated in sequence according to the pouring progress.

[0019] Step 3-3: Use the parametric programming language in the finite element software to apply the concrete heat generation rate at different time periods to the temperature field of the casting bin, then calculate the temperature field of the concrete in the casting bin to obtain the calculated temperature of the concrete;

[0020] Step 4: Use the global artificial bee colony algorithm to iteratively invert the thermal parameters, and use the objective function to judge the quality of the obtained thermal parameter values;

[0021] Step 5: Repeat step 4 to obtain the value of the thermal parameters that meet the objective function value requirements, and compare the calculated temperature value with the actual temperature value of the concrete. Stop the iterative inversion when the expected target is achieved;

[0022] Step 6: Based on the obtained thermal parameter values, predict the temperature of the dam concrete.

[0023] Preferably, in the steps 1-2, the bottom surface of the casting bin and the side surface of the transverse seam are used as adiabatic boundaries, and the top surface of the casting bin and the upstream and downstream surfaces are used as third-type temperature boundary conditions;

[0024] When concrete is in contact with air, it is assumed that the heat flux q passing through the concrete surface is proportional to the concrete surface temperature T and the air temperature T a The expression of the third type of heat transfer boundary condition is as follows:

[0025]

[0026] Where q is the heat flow; λ is the thermal conductivity; T a* is the comprehensive equivalent air temperature; β is the heat release coefficient of the concrete surface; T s is the daily average temperature of the concrete surface; n is the unit vector in the normal direction.

[0027] Preferably, in steps 1-3, the calculation formula for the concrete temperature during cooling and water flow is as follows:

[0028] T(t)=T ωi +(T i -T ωi )f i (t)+θ0ψ i (t) (2)

[0029]

[0030] Where T(t) represents the average temperature of concrete; e is a natural constant; T ωi T is the cooling water temperature of the i-th gear; i is the concrete temperature when the i-1 level water flow ends and the i-th level water flow begins; φ i (t) is the water cooling function when the water is passed through the i-th gear; θ0 is the final adiabatic temperature rise; ψ i (t) is the water-cooling temperature rise function when the water is flowing in the i-th gear; p i is the water cooling parameter when the water is flowing in the i-th gear; t i is the moment when the flow rate or water temperature changes; t is the cooling time; s, m1, and m2 are all coefficients to be determined;

[0031] The calculation formula of water cooling parameters is as follows:

[0032]

[0033] k=2.09-1.35ξ+0.320ξ 2 (7)

[0034] ξ=λL / c w ρ w q w (8)

[0035] The calculation formula for the equivalent thermal conductivity of non-metallic cooling water pipes is:

[0036]

[0037] Where D is the equivalent cooling column diameter; b is the equivalent cooling column radius; S1 and S2 are the horizontal and vertical spacing of the water pipe arrangement respectively; ρ w is the density of water, q w is the water flow rate, L is the length of the water pipe, c w is the specific heat of water, λ is the thermal conductivity of concrete; a' is the equivalent thermal conductivity, c is the outer radius of the water pipe, r0 is the inner radius of the water pipe, λ1 is the thermal conductivity of the water pipe; ξ represents the intermediate variable related to the specific heat of water, water density, water flow rate, and water pipe length.

[0038] Preferably, in step 2-2, the expression of the objective function f(x) is:

[0039]

[0040] min f(x)x=(x1,x2,x3,…,x D )∈[L d ,U d ] (14)

[0041] Where x represents the feasible solution of the optimization problem, x * represents the optimal solution of the optimization problem, x is a D-dimensional vector, and D represents the number of parameters of the optimization problem; T ij '、T ij They represent the calculated value and measured value of concrete temperature at the jth position i at time; p represents the total number of measuring points, q is the total duration; L d 、U d They represent the lower and upper bounds of the solution space respectively.

[0042] Preferably, in step 4, a global artificial bee colony algorithm based on crossover operation is used to iteratively invert the thermal parameters.

[0043] The global artificial bee colony algorithm based on crossover operation addresses the problem of slow convergence of the artificial bee colony algorithm by improving the search equation of the solution space as follows:

[0044]

[0045] Where v id represents the velocity of the i-th particle; x id represents the position of the i-th particle in the search space; represents the dth component in the global optimal solution vector; α is a random number in [-1,1]; β is a random number in [0,C], where C is a non-negative constant;

[0046] pass The item can balance the exploration and development capabilities of the artificial bee colony algorithm, but reduces the global optimization capability. To address this problem, the crossover operation in the genetic algorithm is introduced to further improve the artificial bee colony algorithm.

[0047]

[0048] Where v id ' represents the velocity of the i-th particle after the introduction of the crossover operation; cr represents the crossover operator; rand is a random number in [0,1]; when cr takes a smaller value, the development ability of the improved artificial bee colony algorithm is enhanced, conversely, when cr takes a larger value, the exploration ability of the improved artificial bee colony algorithm is enhanced; corresponding to different optimization problems, cr can take different values ​​to enable the algorithm to achieve the best optimization ability, thereby improving the adaptability of the algorithm to different optimization problems.

[0049] Compared with the prior art, the present invention has the following beneficial effects:

[0050] 1) Taking into account the influence of cooling water flow and external ambient temperature, the present invention applies the ABC algorithm to the inverse calculation of concrete thermal parameters based on the temperature measurement data of the optical fiber buried in the concrete silo, thereby reducing the error between the calculated thermal parameter values ​​and the actual thermal parameter values ​​of the concrete during the construction period.

[0051] 2) The inversion results are of great significance for clarifying the relationship between temperature changes and thermal parameters. The dam temperature is then predicted based on the calculated thermal parameter values, which improves the accuracy of dam temperature prediction and is conducive to the implementation of refined real-time temperature control of the dam.

[0052] 3) In view of the weak development capability of the ABC algorithm, the present invention increases the convergence speed of the ABC algorithm and improves the global optimization capability of the ABC algorithm by adding the guiding role of the current global optimal value in the search process of the ABC algorithm solution space and adding a crossover operation in the search process of the leading bee field. BRIEF DESCRIPTION OF THE DRAWINGS

[0053] The present invention will be further described below with reference to the accompanying drawings and examples.

[0054] Figure 1 Schematic diagram of a distributed optical fiber temperature measurement system according to an embodiment of the present invention.

[0055] Figure 2 Schematic diagram of a thermal parameter inversion framework according to an embodiment of the present invention.

[0056] Figure 3 This is a typical layout diagram of distributed optical fibers in a concrete pouring bin according to an embodiment of the present invention.

[0057] Figure 4 This is a plan view of the distributed optical fiber layout of the 021 concrete pouring bin according to an embodiment of the present invention.

[0058] Figure 5 This is a finite element model of the concrete pouring bin according to an embodiment of the present invention.

[0059] Figure 6 This is a flow chart of the concrete pouring bin temperature simulation pre-processing according to an embodiment of the present invention.

[0060] Figure 7 This is a curve diagram of water temperature change in the concrete pouring bin according to an embodiment of the present invention.

[0061] Figure 8 This is a curve diagram of water flow change in a concrete pouring bin according to an embodiment of the present invention.

[0062] Figure 9This is a curve diagram of the ambient temperature change of the concrete pouring bin according to an embodiment of the present invention.

[0063] Figure 10 This is a graph showing how the objective function of the ABC algorithm according to an embodiment of the present invention changes with the number of iterations.

[0064] Figure 11 This is a graph showing how the objective function of different nectar source quantities changes with the number of iterations in an embodiment of the present invention.

[0065] Figure 12 This is a comparison chart of the temperature prediction value and the measured value of the 021 concrete pouring bin in an embodiment of the present invention.

[0066] Figure 13 This is a comparison chart of the temperature prediction value and the measured value of the 022 concrete pouring bin in an embodiment of the present invention.

[0067] Figure 14 Schematic diagram of component crossover operation of the CGABC algorithm according to an embodiment of the present invention. DETAILED DESCRIPTION

[0068] In this example, the Baihetan concrete double-curvature arch dam, located on the border between Sichuan and Yunnan, is the second hydropower station in the cascade development of the lower Jinsha River. It consists of 31 sections. Located in a subtropical monsoon region, the dam site has an average annual temperature of 21.9°C, with large temperature swings between extremes and significant daytime and nighttime variations. The highest temperature reached 42.7°C in September, while the lowest reached 2.1°C in December. The average water temperature was 17.4°C, and the average wind speed was 1.9 m / s.

[0069] The concrete grades used in this example are primarily C30, C35, and C40. To analyze the temperature variations of C40 concrete in the Baihetan Project under the influence of the external pouring environment, this example selected the 020-022 pouring bins of the 27# dam section. The 020-022 pouring bins are three consecutive bins from bottom to top, and an improved artificial bee colony algorithm was used to invert the concrete thermal parameters.

[0070] In the embodiment, optical fibers are laid out in each section of the dam to measure concrete temperature data. The optical fibers in each section are laid out in an L-shaped pattern. The starting and ending points of the optical fiber laying lines in the casting chamber are 3 to 5 meters away from the horizontal joint template and 3 meters away from the upstream and downstream templates. Figures 3-4 As shown in the figure, the average temperature of the downstream optical fiber in each casting bin is taken as the measured temperature at the center point. Thermometers are also placed at typical measurement points for cross-verification.

[0071] like Figure 1 As shown in FIG, the dam temperature prediction method based on finite element analysis and thermal parameter inversion includes the following steps:

[0072] Step 1: Establish a finite element model for the concrete pouring silo and set initial and boundary conditions;

[0073] Step 1-1: Create a solid model of the casting silo in CAD software. Extract coordinates from this CAD model and import it into ANSYS to create a finite element model of the casting silo. Mesh this finite element model. Assign concrete material and element properties. In ANSYS, select solid70, which supports birth / death elements, as the concrete element.

[0074] Step 1-2: The initial condition mainly refers to the initial temperature of the concrete, which is obtained from the measured data on site. In the embodiment, the bottom surface of the casting bin and the side surface of the transverse joint are used as adiabatic boundaries, and the top surface of the casting bin and the upstream and downstream surfaces are used as third-type temperature boundary conditions.

[0075] When concrete is in contact with air, it is assumed that the heat flux q passing through the concrete surface is proportional to the concrete surface temperature T and the air temperature T a The difference is proportional to the third type of heat transfer boundary condition expression:

[0076]

[0077] In formula (1), q is the heat flow; λ is the thermal conductivity; T a* is the comprehensive equivalent air temperature; β is the heat release coefficient of the concrete surface, unit is kJ / (m 2 ·h·℃); T s is the daily average temperature of the concrete surface; n is the unit vector in the normal direction.

[0078] Steps 1-3: External conditions mainly include ambient temperature and water cooling of concrete. The first stage of water cooling typically uses multiple levels of water flow, with varying water temperatures and flow rates.

[0079] The calculation formula for the average temperature of concrete when cooling water is considered is as follows:

[0080] T(t)=T ωi +(T i -T ωi )φ i (t)+θ0ψ i (t) (2)

[0081]

[0082]

[0083] Where e is a natural constant in mathematics, T ωi is the water temperature of the i-th gear; T iis the concrete temperature when the i-1 level water flow ends and the i-th level water flow begins; φ i (t) is the water cooling function when the water is passed through the i-th gear; θ0 is the final adiabatic temperature rise; ψ i (t) is the water-cooling temperature rise function when the water is flowing in the i-th gear; p i is the water cooling parameter when the water is flowing in the i-th gear; t is the cooling time; t i is the moment when the flow rate or water temperature changes. When the flow rate or water temperature changes, t must start from 0; s, m1, and m2 are all coefficients to be determined.

[0084] The calculation formula of water cooling parameters is as follows:

[0085]

[0086] k=2.09-1.35ξ+0.320ξ 2 (7)

[0087] ξ=λL / c w ρ w q w (8)

[0088] The calculation formula for the equivalent thermal conductivity of non-metallic cooling water pipes is:

[0089]

[0090] In formulas (5) to (8), k, a', and D are constants; D is the equivalent cooling column diameter; b is the equivalent cooling column radius; S1 and S2 are the horizontal and vertical spacings of the water pipes, respectively; ρ w is the density of water, q w is the water flow rate, L is the length of the water pipe, c w is the specific heat of water, λ is the thermal conductivity of concrete; a′ is the equivalent thermal conductivity, which can be obtained by substituting the known equivalent cooling column radius b, water pipe outer radius c, water pipe inner radius r0, concrete thermal conductivity λ, and water pipe thermal conductivity l1 into formula (9).

[0091] In the embodiment, the finite element model of the 27# dam section 020-022 casting bins was established using ANSYS finite element software. The concrete grade of these three continuous casting bins was C40 and the gradation was the same. Each casting bin was divided into 6 0.5m thick layers, and the bin thickness was 3m. The finite element model of the casting bin was meshed using hexahedral eight-node isoparametric elements, with a total of 74976 elements and 83655 nodes. Figure 5 As shown. The density of C40 concrete is 2601kg / m 3 , specific heat c is 0.851kJ / (kg·℃), thermal conductivity α is 0.0078m 2 / d.

[0092] The sides of the silo's transverse joints were treated as adiabatic boundaries, and the silo's top and upstream and downstream surfaces were treated as third-class temperature boundary conditions. The cooling water flow time, water temperature, and flow rate were calculated using a time-based water flow scheme. The ambient temperature was the measured average temperature over the 0.5-day period. To eliminate the impact of silo 019 on the concrete temperature of silo 020, the measured concrete temperature of silo 020 was excluded from the inversion calculation. Following the bottom-up pouring sequence used in actual construction, utilizing the birth and death unit method in ANSYS's thermal analysis process—"killing" all units before calculation—only silo 020 units were activated during pouring in silo 020, only silo 020 and silo 021 units were activated during pouring in silo 021, and all three silo units were activated during pouring in silo 022. The temperature load was applied with a time step of 0.5 days, for a total duration of 30 days.

[0093] The temperature simulation pre-processing flow chart of ANSYS software is as follows Figure 6 The cooling effect of the water pipe is simulated and calculated using the equivalent heat conduction method. The actual cooling water temperature, water flow rate and external ambient temperature information of the water pipe are shown as follows: Figure 7 、 8 , as shown in 9.

[0094] Step 2: Determine the thermal parameters to be inverted and their value ranges;

[0095] Step 2-1: Use the concrete surface heat release coefficient β, the concrete temperature rise law parameter n, and the final adiabatic temperature rise θ0 of concrete as the parameters to be inverted. Considering that the hydration heat of cement is a key factor affecting concrete thermal stress, the actual temperature field calculation uses the adiabatic temperature rise θ of concrete. The calculation of adiabatic temperature rise requires the use of concrete's thermal parameters.

[0096] The adiabatic temperature rise expression of concrete is usually expressed in hyperbolic form, which is as follows:

[0097]

[0098] Where θ(τ) represents the adiabatic temperature rise function of concrete, and τ represents the age of concrete.

[0099] Although the hyperbola formula is highly consistent with the test data, it is inconvenient to calculate. Therefore, Zhu Bofang proposed a combined exponential formula for calculating the adiabatic temperature rise of concrete:

[0100]

[0101] The relationship between the parameters of formula (11) is as follows:

[0102]

[0103] Where θ0 is the final adiabatic temperature rise of concrete, τ is the age of concrete, s, m1, and m2 are three coefficients to be determined. In the embodiment, the value of s is 0.60.

[0104] In the embodiment, a combined exponential formula is selected as the expression for the adiabatic temperature rise of concrete. Furthermore, considering that the main thermal parameters of concrete include the equivalent surface heat release coefficient β, the temperature rise law parameter n, the final adiabatic temperature rise θ0, the thermal conductivity λ, the thermal conductivity α, the density ρ of concrete, and the specific heat c. Among them, the density ρ, specific heat c, and thermal conductivity α of concrete can be obtained with relatively accurate results through indoor tests; the thermal conductivity λ can be determined according to the formula α = λ / cρ. After the concrete is poured, the surface needs to be covered with an insulation blanket and maintained in running water, and the transverse joints are directly exposed to the air, so the equivalent surface heat release coefficient β needs to be inverted.

[0105] Therefore, the present invention selects three thermal parameters, namely, the final adiabatic temperature rise θ0, the temperature rise law parameter n, and the equivalent surface heat release coefficient β, as parameters to be inverted.

[0106] Referring to the concrete design value of the area and the indoor experimental values ​​of the Yangtze River Scientific Research Institute, the range of thermal parameters to be inverted is: θ0∈[18,30], n∈[1.2,6], β∈[100,1200].

[0107] Step 2-2: Establish the objective function;

[0108] In order to minimize the difference between the temperature value of the typical node calculated by the artificial bee colony algorithm and the actual temperature value of the node, that is, when x * =[θ0 * ,n * ,β * ], f(x * )=minf(x), and establish the following objective function expression (13), thereby transforming the inversion problem of concrete thermal parameters into the optimization problem of finding the f(x) function.

[0109]

[0110] min f(x)x=(x1,x2,x3,…,x D )∈[L d ,U d ] (14)

[0111] Where f(x) is the objective function of the intelligent optimization problem; x is the feasible solution of the optimization problem, x * is the optimal solution to the optimization problem, each solution is a D-dimensional vector; D is the number of parameters of the intelligent optimization problem; T′ ij 、T ijThey represent the calculated value and measured value of concrete temperature at the jth position at time i, respectively; p is the total number of measuring points, and q is the total duration. d 、U d are the lower and upper bounds of the solution space respectively.

[0112] Step 3: Substitute a set of randomly selected thermal parameter values ​​within the value range into ANSYS for temperature simulation calculation to obtain the temperature calculation value T′ of each node at different times in the concrete pouring bin. ij Then, the actual temperature and calculated temperature of the concrete are brought into the objective function (13) for calculation to obtain the objective function value.

[0113] Step 3-1: The main task of solving the concrete temperature problem is to simulate the concrete pouring sequence, atmospheric temperature, pouring temperature, and hydration heat. However, since ANSYS software cannot directly apply the hydration heat of concrete, it is necessary to convert the hydration heat into a heat generation rate (HGEN), which is the amount of heat generated per unit volume of concrete per unit time. This can be obtained by differentiating the cement hydration heat formula. Once the heat generation rate is calculated, it is applied to the unit and the temperature field is calculated.

[0114] The calculation formula for the heat generation rate is:

[0115]

[0116] Where HGEN is the heat generation rate of concrete hydration; Q is the heat generated in concrete; p is the density of concrete, unit Kg / m 3 ; c is the specific heat of concrete, unit is kJ / (kg·℃).

[0117] According to formula (15), the heat generation rate of concrete can be determined by the given adiabatic temperature rise expression.

[0118] Step 3-2: Because the concrete pouring bins are not poured all at once but in separate bins, the "element life and death" function in ANSYS is used to implement this process. Before the calculation, all dam elements are "killed." During the calculation, the dam elements are then activated sequentially according to the pouring progress.

[0119] Step 3-3: Apply the concrete heat generation rate at different time periods to the temperature field using the parametric programming language (APDL) in ANSYS software, and then calculate the temperature field of the concrete in the pouring bin to obtain the calculated temperature T' of the concrete. ij '.

[0120] Step 4: The thermal parameters that need to be inverted are iterated through the improved global artificial bee colony algorithm, and the process of step 4 is repeated. The values ​​of the thermal parameters that meet the requirements of the objective function value are obtained by relying on the optimization ability of the improved global artificial bee colony algorithm.

[0121] The standard process of the artificial bee colony algorithm ABC is as follows:

[0122] 1) Initialization phase;

[0123] The value range of the parameter to be inverted is taken as the solution space, and SN nectar sources are randomly generated in the solution space. Each nectar source x i (x i1 ,x i2 ,…,x iD ) T The position (i = 1, 2, ..., SN) is a potential solution. Each nectar source only has one leader bee collecting nectar at a time, where the number of leader bees = the number of follower bees = SN / 2. The initial position of the nectar source i = 1, 2, ..., SN is randomly generated according to the following formula:

[0124] x id =L d +rand(0,1)(U d -L d ) (16)

[0125] Where rand(0,1) represents a random number between 0 and 1.

[0126] At the same time, the maximum number of mining times is initialized to facilitate escaping the local optimal solution later.

[0127] 2) Leading bee stage;

[0128] The leading bee searches the neighborhood of nectar source i according to formula (17) and generates a new nectar source location v i , calculate the fitness value of the new nectar source according to formula (18) i ), using greedy selection method and fitness(x i ) for comparison, if v i The function value is better than x i , the location information of the new nectar source is retained.

[0129]

[0130] In formulas (16)-(18), d is a random integer in [1, D], indicating that i (x i1 ,x i2 ,…,x iD ) TThe random one-dimensional disturbance in is generated; a is the acceleration coefficient, which is usually 1; is a random number in [-1,1], representing the perturbation amplitude; j≠i and j∈{1,2,…,SN}, representing the selection of a nectar source different from nectar source i in the neighborhood. As the search gets closer to the optimal solution, the range of the neighborhood will gradually decrease; x i represents the location of nectar source i; fitness(x i ) is the fitness value of nectar source i; f(x i ) is the function value of nectar source i.

[0131] 3) Follow the bee stage

[0132] The following bee calculates the following probability p through formula (18) i , and the roulette wheel method is used to select the honey source, that is, a random number r is randomly generated in [0,1]. If r<p i , then the follower bee selects nectar source i and performs the same steps as the previous leading bee, that is, turns to formula (16) to generate a new nectar source in the neighborhood and performs greedy selection; if r>p i , then proceed to the next step.

[0133]

[0134] 4) Scouting Bee Stage

[0135] Determine whether the objective function is reached within the maximum number of iterations. If so, the loop ends and the optimal solution is output. Otherwise, the nectar source is abandoned, and the leading bee that collected nectar at the nectar source is transformed into a scout bee. Return to equation (16) and randomly generate a new nectar source to replace it. The loop continues until the optimal solution is obtained. The optimal solution at this time is the optimal value of the inversion parameter.

[0136] The global artificial bee colony algorithm (GABC) of the embodiment is as follows:

[0137] In the algorithm cycle from step 1) to step 4), there is no guidance for the memory of the current global optimal value, only a random component x i , which leads to insufficient development capability of ABC algorithm and slow convergence speed. To solve this problem, the present invention uses the global optimal solution (gbest) in particle swarm algorithm to improve the search equation (17) in standard artificial bee colony algorithm.

[0138]

[0139] In formula (20) is the dth component of the algorithm's current global optimal solution vector; α is a random number in [-1, 1]; β is a random number in [0, C], where C is a non-negative constant. When C is 0, Equation (20) is equal to Equation (17). The algorithm performance analysis test with different values ​​of C shows that the algorithm performance is better when C is 1.5, and the global optimization ability of the algorithm is enhanced.

[0140] Compared with formula (17), formula (20) is The item can balance the exploration and development capabilities of the ABC algorithm, but it reduces the global optimization capability to a certain extent. To address this problem, the present invention introduces the crossover operation in the genetic algorithm, further improves the GABC algorithm, and proposes a global artificial bee colony algorithm (Cross-Global Artificial Bee Colony algorithm, CGABC) based on the crossover operation.

[0141] The crossover operation is to exchange the partial code values ​​of two parents bit by bit to obtain a new individual. In the algorithm, each component first generates a uniform random number rand between 0 and 1. If rand < the set value cr, the target component is accepted, otherwise the current component is retained, such as Figure 14 shown.

[0142] The present invention combines the GABC algorithm with the two-term crossover in the genetic algorithm. By allowing the leading bee to search the neighborhood and then perform a crossover operation with the current global optimal value, that is, adding formula (21) after formula (20) to improve the global optimization ability of the algorithm:

[0143]

[0144] In formula (21), cr is the crossover operator set in the algorithm, which generally takes a value of 0.3 to 0.6; rand is a uniform random number between [0, 1]. When cr takes a smaller value, the development ability of the algorithm is enhanced. Conversely, when cr takes a larger value, the exploration ability of the algorithm is enhanced.

[0145] Corresponding to different optimization problems, cr can take different values ​​to enable the algorithm to achieve the best optimization ability, thereby improving the adaptability of the algorithm to different optimization problems.

[0146] Figure 10 The objective function comparison values ​​of the standard ABC algorithm and the CGABC algorithm of the embodiment at the same number of iterations in the inversion calculation of the thermal parameters of C40 concrete are given. It can be seen that the objective function values ​​of both algorithms gradually decrease with the increase of the number of iterations, but the CGABC algorithm converges faster than the ABC algorithm. At the same number of iterations, the objective function values ​​are mostly smaller than those of the ABC algorithm, that is, the inversion calculation results are better, indicating that the CGABC algorithm has better performance in the inversion of thermal parameters.

[0147] Figure 11 The performance of the improved CGABC algorithm in the inversion calculation of thermal parameters of C40 concrete at different colony sizes is given. The maximum number of iterations maxCycle=200, the maximum number of mining times limit=100, and the dimension D=3 are set. By selecting different numbers of artificial bee colony nectar sources SN for parameter inversion calculation, the change of the objective function value with the number of iterations is obtained when the number of nectar sources is 20, 40, and 60, respectively.

[0148] Step 5: Bring back the thermal parameters obtained by inversion for inverse analysis;

[0149] In the embodiment, the thermal parameters of C40 concrete obtained by inversion when the number of honey sources is 40 are: θ0 = 24 ° C, n = 2.86, β = 489.2 kJ / (m 2 ·d·℃), as shown in Table 1, and substituted it into the temperature field of 021 and 022 casting bins of 27# dam section for simulation calculation, and the calculated temperature values ​​of concrete at different ages were obtained. The calculated temperature and measured temperature change trends of 021 casting bin and 022 casting bin are compared as shown in Figure 12 、 13 As shown. Figure 12 and Figure 13 It can be seen that the measured values ​​of the concrete temperature in the 021 and 022 casting bins are in good agreement with the calculated values, and the change trends are consistent. The differences at individual time points in the 022 casting bin are due to the unstable cooling water flow during the actual construction process.

[0150] Table 1 Thermal parameter inversion value table

[0151]

[0152] Step 6: Based on the obtained thermal parameter values, predict the temperature of the dam concrete.

[0153] To verify that the CGABC algorithm of the present invention has stronger optimization capabilities than the pre-improved ABC and GABC algorithms, five commonly used performance test functions, Rastrigin, Sphere, Rosenbrock, Ackely, and Schaffer, were selected to perform optimization performance tests on the ABC, GABC, and CGABC algorithms, respectively. The Sphere function is a unimodal function, and the other four test functions are multimodal functions. The theoretical optimal value of all test functions is 0. The expressions and search ranges of these five performance test functions are shown in Table 2.

[0154] Table 2 Performance test function table

[0155]

[0156] In the performance comparison test, the number of nectar sources = 50, the number of leading bees = the number of following bees = 25, the maximum number of iterations = 5000, the maximum number of mining times = 300, the dimension D of the Rosenbrock function = 2, the dimension D of the Sphere, Rastrigin, Ackely, and Schaffer functions are all 30, and the optimization performance is better when the crossover operator in the CGABC algorithm is 0.5. All test functions are repeated 10 times using each algorithm under the above-mentioned set parameters. The algorithm programs are all written in Matlab language, and the mean value and standard deviation SD of each function are obtained.

[0157] Table 3 compares the performance of the standard artificial bee colony algorithm ABC, the improved global artificial bee colony algorithm GABC, and the crossover-based global artificial bee colony algorithm CGABC on five different test functions. The test results in Table 3 show that the CGABC algorithm has a stronger optimization capability than the ABC and GABC algorithms under the same set conditions. By introducing the global optimal solution gbest into the steps of the ABC algorithm, the CGABC algorithm accelerates the algorithm's convergence and, combined with the crossover operation, significantly improves the algorithm's search breadth and depth.

[0158] Table 3 Comparison of three artificial bee colony algorithms

[0159]

[0160] In summary, it is shown that the artificial bee colony algorithm can be applied to the inversion of dam thermal parameters, and the inversion results are highly reliable, thus laying the foundation for the subsequent real-time control of dam concrete temperature.

Claims

1. A dam temperature prediction method based on finite element analysis and thermal parameter inversion is characterized by: The following steps are involved: Step 1: Establish the finite element model of each concrete pouring bin and set the initial conditions and boundary conditions; Step 1-1: In the finite element software, establish a finite element model of the casting silo and mesh the finite element model; at the same time, assign values ​​to the concrete material and unit properties; Step 1-2: The initial conditions of the finite element model include the initial temperature of the concrete, which is obtained based on the measured data on site; Steps 1-3: Boundary conditions include ambient temperature and water cooling of concrete; Step 2: Determine the thermal parameters and their value ranges used for dam temperature prediction, and establish the objective function of thermal parameters; Step 2-1: Determine the inverse thermal parameters and their value ranges. The inverse thermal parameters include the concrete surface heat release coefficient, concrete temperature rise law parameters, and the final adiabatic temperature rise of concrete. Step 2-2: Establish the objective function of thermal parameter inversion and transform the concrete thermal parameter inversion problem into an optimization problem; Step 3: Randomly generate a set of initial thermal parameter values ​​within the range of thermal parameter values ​​and substitute them into the finite element model for temperature simulation calculation to obtain the temperature values ​​at the nodes in the concrete pouring bin. Substitute them into the objective function together with the actual temperature of the concrete at the nodes to calculate the objective function value; Step 4: Use the global artificial bee colony algorithm to iteratively invert the thermal parameters, and use the objective function to judge the quality of the obtained thermal parameter values; Step 5: Repeat step 4 to obtain the value of the thermal parameters that meet the objective function value requirements, and compare the calculated temperature value with the actual temperature value of the concrete. Stop the iterative inversion when the expected target is achieved; Step 6: Based on the obtained thermal parameter values, predict the temperature of the dam concrete.

2. The dam temperature prediction method according to claim 1, characterized in that: In the steps 1-2, the bottom surface of the casting bin and the side surface of the transverse seam are used as adiabatic boundaries, and the top surface of the casting bin and the upstream and downstream surfaces are used as third-type temperature boundary conditions; When concrete is in contact with air, the heat flux through the concrete surface is assumed to be Concrete surface temperature and temperature The expression of the third type of heat transfer boundary condition is as follows: ； (1) In the formula is the heat flow; is the thermal conductivity; is the comprehensive equivalent temperature; is the heat release coefficient of the concrete surface; is the daily average temperature of the concrete surface; is the normal direction unit vector.

3. The dam temperature prediction method according to claim 2, characterized in that: In steps 1-3, the calculation formula for the concrete temperature during cooling and water flow is as follows: ；(2) ；(3) ； (4) In the formula Indicates the average temperature of concrete; is a natural constant; For the Cooling water temperature; for The water flow of the gear is completed and the Concrete temperature when water flow begins; For the Water cooling function when the gear is open; is the final adiabatic temperature rise; For the Water cooling temperature rise function when the gear is open to water; For the Water cooling parameters when the gear is open; The moment when the flow rate or water temperature changes; For cooling time; All are coefficients to be determined; The calculation formula of water cooling parameters is as follows: ；(5) ；(6) ；(7) ；(8) The calculation formula for the equivalent thermal conductivity of non-metallic cooling water pipes is: ； (9) In the formula is the equivalent cooling column diameter; is the equivalent cooling column radius; They are the horizontal and vertical spacing of the water pipe arrangement respectively; is the density of water, is the water flow rate, is the length of the water pipe, is the specific heat of water, is the thermal conductivity of concrete; is the equivalent thermal conductivity, is the outer radius of the water pipe, is the inner radius of the water pipe, is the thermal conductivity of the water pipe; Represents intermediate variables related to water specific heat, water density, water flow rate, and water pipe length.

4. The dam temperature prediction method according to claim 3, characterized in that: In step 2-2, the objective function The expression is: (13) ； (14) In the formula x represents a feasible solution to the optimization problem, x for dimensional vector, The number of parameters representing the optimization problem; 、 Respectively represent j Locations i Calculated and measured values ​​of concrete temperature at the moment; p Indicates the total number of measuring point locations, q is the total duration; 、 They represent the lower and upper bounds of the solution space respectively.

5. The dam temperature prediction method according to claim 4, characterized in that: Step 3 includes the following sub-steps: Step 3-1: Simulate the concrete pouring sequence, atmospheric temperature, pouring temperature, and hydration heat; Step 3-2: Use the "element generation" and "element killing" functions of the finite element software to simulate the concrete pouring in compartments. Before the simulation of the compartment pouring, all dam body elements are "killed". During the simulation, the dam body elements are activated in sequence according to the pouring progress. Step 3-3: Use the parametric programming language in the finite element software to apply the concrete heat generation rate in different time periods to the temperature field of the casting bin, then calculate the temperature field of the concrete in the casting bin to obtain the calculated temperature of the concrete.

6. The dam temperature prediction method according to claim 5, characterized in that: In step 3-1, the hydration heat is converted into a heat generation rate. The heat generation rate is the amount of heat generated per unit volume of concrete per unit time and is obtained by deriving the cement hydration heat formula. After the heat generation rate is calculated, the heat generation rate is applied to the unit and the temperature field is calculated. The heat generation rate is calculated as follows: ； (15) In the formula is the heat generation rate of concrete hydration; is the heat generated in the concrete; is the density of concrete; is the specific heat of concrete; It represents the adiabatic temperature rise of concrete; represents the derivative of the adiabatic temperature rise function of concrete; Indicates the age of concrete; According to formula (15), the heat generation rate of concrete is determined by the adiabatic temperature rise expression.

7. The dam temperature prediction method according to claim 6, characterized in that: In step 4, a global artificial bee colony algorithm based on crossover operation is used to iteratively invert the thermal parameters. The global artificial bee colony algorithm based on crossover operation addresses the problem of slow convergence of the artificial bee colony algorithm by improving the search equation of the solution space as follows: ； (20) In the formula Indicates the i The speed of a particle; Indicates the i The position of each particle in the search space; represents the first d Quantity is a random number in [-1,1]; for A random number in C is a non-negative constant; Formula (20) is obtained by The item can balance the exploration and development capabilities of the artificial bee colony algorithm, but reduces the global optimization capability. To address this problem, the crossover operation in the genetic algorithm is introduced to further improve the artificial bee colony algorithm. ； (21) In the formula Indicates that after the crossover operation is introduced i The speed of a particle; Represents the crossover operator; rand is a random number in [0,1]; when When the value is smaller, the development ability of the improved artificial bee colony algorithm is enhanced. On the contrary, when When the value is larger, the exploration ability of the improved artificial bee colony algorithm is enhanced; corresponding to different optimization problems, Different values ​​can be taken to enable the algorithm to achieve the best optimization ability, thereby improving the adaptability of the algorithm to different optimization problems.

Images