Dam creep parameter adaptive inversion method based on variable order fractional constitutive model

Abstract

The application discloses a dam creep parameter self-adaptive inversion method based on a variable-order fractional constitutive model, and comprises the following steps: acquiring measured displacement time series of multiple measuring points and load and environmental action time courses, and constructing a three-dimensional finite element model; constructing a variable-order fractional Burgers constitutive relationship reflecting the activity attenuation characteristics of concrete creep; discretizing a dam service period into multiple subintervals, and applying a superposition principle to calculate the partial strain increment of each subinterval in a segmented constant-order framework to obtain a variable-order fractional forward calculation model; constructing an inversion objective function, dividing the parameters to be inverted into a mechanical parameter group and a fractional-order evolution parameter group; alternately optimizing the mechanical parameter group and the fractional-order evolution parameter group to obtain and output a final inversion parameter set. The application solves the conflict between the aging material and the superposition principle, eliminates the inconsistency defects of the heterogeneous parameter joint inversion, and improves the physical rationality and numerical stability of the long-term service state prediction of the dam.

Description

Technical Field

[0001] This invention relates to dam safety monitoring and numerical simulation technology, and in particular to an adaptive inversion method for dam creep parameters based on a variable-order fractional-order constitutive model. Background Technology

[0002] Concrete dams undergo irreversible creep deformation under long-term water storage and environmental loads. Accurately identifying creep characteristic parameters is the technical foundation for assessing the structural health and long-term stability of dams. Extracting the underlying material mechanical properties from massive amounts of measured displacement data using inversion algorithms is a core technical step in constructing high-fidelity structural mechanics calculation models and realizing forward-looking extrapolation prediction of dam deformation behavior and long-term safety early warning.

[0003] Existing techniques for inverting dam creep parameters typically employ constant-order fractional-order rheological models to describe the viscoelastic mechanical behavior of concrete. In numerical forward modeling, current methods primarily utilize time-domain strain convolution integrals using creep compliance kernel functions of a fixed order. During parameter identification, a common approach is to construct a fitting objective function based on measured displacement residuals, apply uniformly weighted prior constraints using standard regularization methods, and employ optimization algorithms to iteratively search within a given fixed parameter space to seek the optimal combination of mechanical parameters such as the elastic modulus and constant-order viscosity coefficient.

[0004] Existing inversion techniques face several drawbacks when dealing with long-term dams. These include the difficulty in characterizing aging degradation with fixed-order models, solution divergence due to heterogeneous parameter dimensions and sensitivity, and the violation of the superposition principle by integrals in continuously aging materials. Concrete exhibits significant aging effects during long-term hydration and densification, with its creep activity decreasing over time. Existing fixed-order models struggle to map the time-varying degradation of this material's mechanical behavior. Directly introducing order evolution parameters introduces a mixture of pressure-based mechanical parameters and dimensionless evolution parameters with strictly limited values ​​into the inversion space. Traditional uniform weighting regularization methods lead to dimensional inconsistencies and sensitivity imbalances, resulting in over-constraint of high-sensitivity parameters and drastic drift of low-sensitivity parameters. Furthermore, directly using traditional stationary convolution kernels for full-time integration of continuously variable-order aging materials mathematically violates the Boltzmann superposition principle, which applies only to non-aging materials, causing logical conflicts in the underlying mechanical mechanisms. These related issues make it difficult for existing inversion methods to provide stable and reliable long-term extrapolation predictions. Summary of the Invention

[0005] The purpose of this invention is to provide an adaptive inversion method for dam creep parameters based on a variable-order fractional-order constitutive model, in order to solve at least one of the aforementioned problems in the existing technology.

[0006] Technical solution: An adaptive inversion method for dam creep parameters based on a variable-order fractional-order constitutive model, comprising:

[0007] Obtain the time series of measured displacements at multiple measuring points and the time history of load and environmental interaction, and construct a three-dimensional finite element model;

[0008] Based on the material partitioning information of the three-dimensional finite element model, a variable-order fractional-order Burgers constitutive relation reflecting the creep activity decay characteristics of concrete is constructed, which includes an order evolution model.

[0009] The service life of the dam corresponding to the three-dimensional finite element model is obtained, and it is discretized into multiple sub-intervals. Under the piecewise constant order framework, the superposition principle is applied to each sub-interval. Based on the variable order fractional order Burgers constitutive relation, the time history of load and environmental action, and the three-dimensional finite element model, the deviatoric strain increment of each sub-interval is calculated to obtain the variable order fractional order forward model.

[0010] Construct an inversion objective function and divide the parameters to be inverted into a mechanical parameter set and a fractional-order evolution parameter set;

[0011] Based on the measured displacement time series at multiple measurement points, the variable-order fractional-order forward model and the inversion objective function, the mechanical parameter set and the fractional-order evolution parameter set are alternately optimized to obtain and output the final inversion parameter set.

[0012] Beneficial effects: This invention resolves the conflict between aging materials and the superposition principle, eliminates the inconsistency defects in the joint inversion of heterogeneous parameters, and improves the physical rationality and numerical stability of the prediction of the long-term service performance of dams. Attached Figure Description

[0013] Figure 1 This is a schematic diagram of the overall process of the adaptive inversion method for dam creep parameters based on a variable-order fractional-order constitutive model provided in this application embodiment.

[0014] Figure 2 This is a schematic diagram of the process for calculating the deviatoric strain increment of each sub-interval based on the variable-order fractional-order Burgers constitutive relation, the time history of load and environmental action, and a three-dimensional finite element model, as provided in the embodiments of this application.

[0015] Figure 3 This is a schematic diagram of the process for determining the endpoints of non-uniform sub-intervals using an equal error allocation criterion based on the order of change rate, provided in an embodiment of this application.

[0016] Figure 4 This is a schematic diagram of the process for performing dimensionless normalization on each parameter to be inverted in the mechanical parameter set and the fractional-order evolution parameter set, provided in an embodiment of this application. Detailed Implementation

[0017] This invention addresses the challenge of accurately predicting the creep characteristics of dam concrete by proposing an adaptive inversion method based on a variable-order fractional-order constitutive model. A variable-order fractional-order Burgers constitutive model reflecting the creep activity decay characteristics of concrete is constructed. By introducing an evolutionary equation of order, the aging behavior of the material during long-term service is accurately described. In numerical computation, the complex variable-order problem is transformed into a piecewise constant-order problem. Combining the superposition principle, the elastic and creep strain increments of each sub-interval are solved, overcoming the difficulty of solving variable-order calculus in finite element analysis. Regarding the parameter inversion strategy, the proposed method divides the parameters to be inverted into a mechanical parameter group and a fractional-order evolutionary parameter group, and performs alternating optimization, reducing the complexity of high-dimensional nonlinear inversion and achieving accurate identification of dam creep parameters and high-precision prediction of long-term performance.

[0018] Example 1: This section details the overall architecture and macroscopic data flow process of the adaptive inversion method for dam creep parameters based on a variable-order fractional-order constitutive model, such as... Figure 1 As shown.

[0019] Step 101: Obtain the measured displacement time series of multiple measurement points and the time history of load and environmental interaction, and construct a three-dimensional finite element model.

[0020] Specifically, the multi-point measured displacement time series refers to deformation observation records extracted from the dam safety monitoring system, covering multiple complete water storage and discharge cycles. Its physical representation is a multidimensional matrix or data stream containing spatial coordinates (x, y, z) and timestamp sequences. The load and environmental action time history includes external excitation data strictly aligned with the displacement observation period. Its representation is a discrete time series array, specifically including daily recorded upstream and downstream reservoir water level elevations, temperature vectors of grid nodes inside and on the dam body, and self-weight density parameters of the dam body concrete and foundation rock. During data acquisition, the system performs data alignment and denoising operations to ensure data quality. The three-dimensional finite element model is the underlying spatial physical carrier generated by meshing the dam body and foundation system using the finite element numerical method based on the dam's geometric configuration and geological zoning information. Material zoning markings and displacement constraint boundary conditions are applied internally, providing a standardized geometric boundary framework for subsequent parameter calculations.

[0021] In some alternative implementations, the load and environmental action time history can further include time-varying meteorological data such as reservoir water temperature field distribution and ambient wind speed. By enriching the dimensions of external environmental observation, the physical completeness of the input boundary of the finite element space can be effectively improved.

[0022] Step 102: Based on the material partitioning information of the three-dimensional finite element model, construct a variable-order fractional Burgers constitutive relation that reflects the creep activity decay characteristics of concrete. The variable-order fractional Burgers constitutive relation includes an order evolution model.

[0023] In this embodiment, the variable-order fractional-order Burgers constitutive relation is used to describe the viscoelastic-plastic mechanical behavior of dam materials under long-term stress at a macroscopic level. Specifically, to overcome the difficulty of traditional constant-order models in reflecting the aging effects of concrete materials caused by hydration reactions and microstructural densification, this scheme extends the fractional-order number of the Maxwell branch in the classic Burgers model from a fixed constant to a function that dynamically changes with service time. This dynamic change is then mathematically quantified by the order evolution model. By introducing the order evolution model, this constitutive relation can accurately characterize the real physical process of the concrete dam body gradually weakening its viscous flow characteristics and continuously decaying its creep activity as its service life increases.

[0024] Furthermore, the specific mathematical expression of the order evolution model is not limited to a single function type. It can be implemented by adopting an exponential decay analytical function driven solely by absolute service time, or by combining an equivalent hydration age driving function with the spatial temperature variation history of the dam, depending on the specific degradation characteristics of the engineering materials.

[0025] Step 103: Obtain the dam service life corresponding to the three-dimensional finite element model. Discretize the dam service life into multiple sub-intervals. Apply the superposition principle to each sub-interval within a piecewise constant-order framework. Calculate the deviatoric strain increment of each sub-interval based on the variable-order fractional-order Burgers constitutive relation, the load and environmental action time histories, and the three-dimensional finite element model to obtain the variable-order fractional-order forward model. In other words, the variable-order fractional-order forward model can also be considered a variable-order fractional-order forward analysis model.

[0026] Specifically, this step constructs the forward modeling physics engine of this invention. Discretizing the dam's service life into multiple sub-intervals transforms the continuous variable-order calculus problem in the time domain into a finite number of piecewise constant-order sub-problems. Under the piecewise constant-order framework, the fractional order within each independent sub-interval is frozen as a constant, indicating that the material strictly satisfies the non-aging time-invariant condition within a single sub-interval, allowing the Boltzmann superposition principle in classical mechanics to be legally applied within each sub-interval. Based on this theory, the algorithm applies the load and environmental action time histories to the three-dimensional finite element model and, combined with the constitutive parameters after constant-order processing, calculates the partial strain increments at each node through piecewise convolution integration. The above-mentioned complete numerical computation pipeline, integrating mesh boundaries, environmental input, piecewise constitutive model, and strain integration rules, together constitutes the variable-order fractional-order forward model.

[0027] Optionally, the division of multiple sub-intervals can adopt a uniform division strategy of the time axis, or, in view of the nonlinear characteristics of severe creep activity decay in the early stage of dam service, a non-uniform adaptive grid division strategy based on error distribution principles such as order change rate can be adopted to maximize the time domain integration accuracy while controlling the global computation time cost.

[0028] Step 104: Construct the inversion objective function and divide the parameters to be inverted in the variable-order fractional-order Burgers constitutive relation and the variable-order fractional-order forward calculation model into a mechanical parameter group and a fractional-order evolution parameter group.

[0029] In this embodiment, to transform parameter identification into a quantifiable optimization mathematical problem, the system constructs an inversion objective function. This function primarily measures the fitting error between the forward model's calculated output and the actual monitoring data, as well as the penalty for parameter values ​​deviating from prior knowledge. The mathematical expression of the inversion objective function is a weighted combination of the data fitting term and the prior constraint terms for each parameter group, specifically in the form:

[0030] ;

[0031] Where J(p) is the inversion objective function, p is the parameter vector to be identified, and F d λ is a dimensionless data fitting term used to measure the degree of fit between the calculated output of the forward model and the actual monitoring data. g F is the grouping regularization parameter corresponding to the g-th parameter group. prior_g This is a dimensionless quadratic form of the parameter deviation for the g-th parameter group. The specific calculation formulas for each component are detailed in subsequent embodiments.

[0032] Meanwhile, the system strictly divides all unknown material properties into two parameter groups with distinctly different physical properties. The mechanical parameter group includes traditional mechanical scalars such as elastic modulus and viscosity coefficient, whose values ​​typically range from thousands to tens of thousands of megapascals and have a direct and highly sensitive impact on the structural stiffness matrix. The fractional evolution parameter group, on the other hand, contains various dimensionless coefficients that define the order evolution model. Their values ​​are strictly constrained within a physical open interval of 0 to 1 and have lower sensitivity to the final displacement response.

[0033] As an improvement to the above scheme, when constructing the inversion objective function, a Bayesian maximum a posteriori estimation fitting term based on the data error covariance matrix can be introduced to incorporate the sensor observation noise level and statistical uncertainty carried by the measured displacement time series of multiple measurement points into the optimization penalty system.

[0034] Step 105: Based on the measured displacement time series of multiple measurement points, the variable-order fractional forward model and the inversion objective function, the mechanical parameter set and the fractional evolution parameter set are alternately optimized to obtain the final inversion parameter set.

[0035] Specifically, this step executes a closed-loop optimization solution for parameter inversion. The alternating optimization strategy refers to the system first fixing the current estimated values ​​of the fractional-order evolving parameter set in a multidimensional, highly non-convex parameter search space, and then using a single driving force to perform an extremum search within the subspace of the mechanical parameter set. After convergence, the updated mechanical parameter set is fixed, and the system then optimizes the fractional-order evolving parameter set. This layered, decoupled, alternating iterative process effectively avoids the local minima trap that easily occurs when optimizing mixed parameters of high dimensions and magnitudes. Within each iteration, the system repeatedly calls the variable-order fractional-order forward model to output simulated displacements and compares the residuals with the measured displacement time series from multiple measurement points until the relative changes in parameters and the decrease in the objective function between two adjacent iterations are both below the pre-configured convergence threshold. At this point, the corresponding complete set of parameters is extracted as the final inversion parameter set.

[0036] In another specific embodiment, during the cycle of performing alternating optimization, the dimensionless dynamic regularization weights can be independently applied to adjust the extreme differences in sensitivity between the two sets of parameters, ensuring that the high-sensitivity parameter can freely follow the measured data, while the low-sensitivity parameter is reasonably anchored within the physical prior interval.

[0037] In other alternative implementations, the global search engine used in the alternating optimization process is not limited to genetic algorithms; it can also employ other global optimization methods such as differential evolution, particle swarm optimization, or simulated annealing. Similarly, the local fine-tuning engine is not limited to quasi-Newton algorithms; it can also employ the Levenberg-Marquardt algorithm, conjugate gradient method, or trust region method. The combination strategy of the above-mentioned global search and local fine-tuning can be flexibly selected according to the parameter dimensions and non-convexity of the specific problem.

[0038] Step 106: Output the final inversion parameter set.

[0039] In this embodiment, the final inversion parameter set represents the globally optimal material mechanics profile extracted from massive dam safety monitoring data. The output of the final inversion parameter set provides external engineering application terminals with specific mechanical modulus values ​​and precise coefficients of the evolution equations characterizing long-term aging.

[0040] Correspondingly, after the output operation is completed, the final inversion parameter set can be directly loaded into large-scale commercial finite element software or structural health assessment system as the underlying material card benchmark. This is used to make forward-looking extrapolation predictions of the deformation evolution trend of the dam under future extreme water level conditions or ultra-long service cycles, providing data decision-making basis with clear physical mechanism support for dam reinforcement and full life cycle safety early warning.

[0041] Example 2 details the specific implementation process of applying the superposition principle interval by interval within a piecewise constant-order framework. It resolves the theoretical conflict between variable-order aging materials and the classical non-aging superposition principle, supplements the underlying computational details of the finite element time-stepping algorithm and the numerical approximation of predetermined mathematical functions, and ensures the feasibility and computational stability of the physics engine. Figure 2 As shown.

[0042] Step 201: Based on the variable-order fractional-order Burgers constitutive relation, determine the constant-order creep compliance kernel function for each subinterval.

[0043] Specifically, this step is used to eliminate the time-varying parameter barrier in the global integration. Traditional variable-order models represent material parameters that change continuously over time, falling under the category of aging materials, while the integral form in classical mechanics strictly relies on the time-invariant premise of non-aging. Therefore, within each pre-divided discrete sub-interval, the algorithm extracts the order value at the starting point of the interval as a global representative value, substitutes it into the compliance equation of the constitutive relation, and reduces the dimensionality of the complex variable-order time-varying kernel function, instantiating it into a constant-order creep compliance kernel function specific to that interval. In this embodiment, the specific form of the creep compliance kernel function of the fractional-order Burgers model is the compliance expression of a fractional-order four-element rheological model, which is well-known in the art. Those skilled in the art can directly write the corresponding analytical expression of creep compliance based on the selected specific variant of the Burgers model.

[0044] By freezing the order, the material approximately degenerates into a non-aged material within this short time span, providing a rigorous physical and mathematical basis for the subsequent introduction of integral superposition. Furthermore, the constant-order creep compliance kernel function contains a highly complex single-parameter Mittag-Leffler function. To achieve efficient numerical solutions in the underlying finite element program and avoid the risk of long-term integral divergence, the Padé approximation algorithm is used to truncate the function through rational polynomial expansion during the specific calculation.

[0045] Specifically, for the single-parameter Mittag-Leffler function E α The numerical calculation of (-z) employs the [M / N]-order Padé rational approximation method, where the typical configurations of M and N are M=N=4 or M=N=5. This is achieved by subtracting E at z=0. α The Taylor series expansion of (-z) takes up to the first (M+N) terms, and the numerator is a polynomial of degree M, P. M (z) and the Nth degree polynomial in the denominator Q N Solving the system of equations with undetermined coefficients for (z) yields the rational approximation E. α (-z)≈P M (z) / Q N(z). Where α is the order parameter, [M / N] is the order configuration of the Pad approximation, M and N are the highest degree of the polynomial, and (M+N) terms are the number of truncation terms in the Taylor series expansion.

[0046] In some alternative implementations, the numerical inverse Laplace transform method based on the optimal parabolic contour integral, proposed by Garrappa et al., can also be used to perform high-precision calculations of the Mittag-Leffler function.

[0047] By approximating the infinite series form of the Mittag-Leffler function as a ratio of the numerator to the denominator polynomial, the computational efficiency and floating-point stability of the kernel function are improved when it is repeatedly called frequently within millions of discrete time steps.

[0048] Step 202: Based on the time history of load and environmental effects and the three-dimensional finite element model, determine the deviatoric stress increment in each sub-interval.

[0049] In this embodiment, this step transforms the external physical boundary conditions into the actual stress response characteristics within the structure. Specifically, the three-dimensional finite element model constitutes the spatial skeleton for stiffness transfer, while the load and environmental interaction time histories provide the external nodal force vectors that dynamically evolve over time. Under the standard creep incremental method solution framework, the system adopts an implicit time-stepping scheme, gradually applying the external load increments of the current sub-interval to the overall equilibrium equations. For the complex nonlinear loading process experienced by the dam structure during the impoundment period, the nonlinear iteration of the Jacobian matrix accurately solves for the three-dimensional stress tensor changes at each Gaussian integral point within the current time step. Based on the tensor decomposition principle, the algorithm strips the volumetric stress component from the calculated total stress tensor increment, extracting the pure shear component driving the viscous flow of the material, thereby obtaining the deviatoric stress increments within each sub-interval.

[0050] For example, in a large-scale mesh model containing 100,000 elements, the deviatoric stress increment extracted in a certain iteration can be structured into a dense matrix of dimension (L,6), where L represents the total number of Gaussian integration points, and each row of the matrix strictly corresponds to the six independent stress component increments of a single spatial integration point in the current time step. In some optional implementations, for extreme temperature drops or sudden drops in water level, the deviatoric stress increment can also be locally plastically reduced and corrected according to preset tensile damage envelope parameters to more realistically reflect the stress redistribution state in the microcrack development region.

[0051] Step 203: For any k-th subinterval, perform a constant-order convolution integral using the constant-order creep compliance kernel function of each preceding subinterval and the deviatoric stress increment in each preceding subinterval, and sum the results to obtain the preceding cumulative creep deviatoric strain.

[0052] Specifically, this step accurately quantifies the residual memory impact of historical loads on the deformation of the dam material at the current moment. Within a rigorous mathematical piecewise framework, the Boltzmann superposition principle holds locally because the material properties within each sub-interval have been normalized. For the currently solving k-th sub-interval, the computational engine traverses all historical load states from the initial moment to the (k-1)-th sub-interval. For any j-th historical sub-interval, the system retrieves the constant-order creep compliance kernel function specific to that interval and performs a rigorous convolution integration operation with the actual deviatoric stress increment within that interval to calculate the residual creep effect induced by that historical load increment at the current moment. By linearly summing the local integration results of all preceding historical sub-intervals, the preceding cumulative creep deviatoric strain representing the total historical memory can be obtained. To intuitively demonstrate this historical memory accumulation and transformation process, the computational logic is expressed as the following linear formula:

[0053] ;

[0054] Among them, e history (t) represents the cumulative creep strain from the preceding sequence. This is the summation operator for the preceding 1 to (k-1)th subintervals. For the definite integral operator in the time domain of the j-th subinterval, J αj (t-τ) represents the fractional order α of the j-th subinterval. j The kernel function is a constant-order creep compliance function with constant-order parameters. Its time variable (t-τ) represents the actual creep duration from the historical loading time τ to the current observation time t. j Let t be the endpoint time of the j-th sub-interval. j-1 Let τ be the starting time of the j-th subinterval, and s be the internal variable of the integral. j (τ) represents the deviatoric stress component within the j-th historical sub-interval, ds j (τ) / dτ is the derivative of the deviatoric stress component with respect to time at the corresponding historical moment within the j-th subinterval.

[0055] The above piecewise accumulation formula implicitly adopts the classic loading age hypothesis when dealing with long-term evolution spanning different intervals. In other words, the stress increment excited in the j-th historical sub-interval remains rigidly anchored to the material order state at the moment of application throughout its subsequent long creep evolution process, without changing with the latest order at the current moment. This physical assumption effectively alleviates the integral singularity curse that can easily lead to abrupt parameter changes across intervals, enabling numerical calculations at the engineering scale to maintain the rationality of the physical mechanism while possessing high computational efficiency.

[0056] Step 204: Perform a constant-order convolution integral using the constant-order creep compliance kernel function of the k-th sub-interval and the deviatoric stress increment in the k-th sub-interval to obtain the creep deviatoric strain increment of the current interval.

[0057] In this embodiment, this step is used to capture the instantaneous delayed deformation response excited by the newly added physical load in the latest time step. Specifically, the finite element solver directly calls the constant-order creep compliance kernel function frozen for the current k-th sub-interval in step 201, using it as the sole integration operator, and performs a convolution operation with the deviatoric stress increment extracted in step 202 within the k-th sub-interval. This is based on the fact that, since the time span of the current discrete time step is usually controlled within an extremely small range, the degree of creep activity decay of the dam material is negligible within this small span. Using the constant-order kernel function of the current interval for calculation can reconstruct the true continuous-order transient response with high approximation accuracy. The value obtained through this integration operation is the creep deviatoric strain increment of the current interval.

[0058] ;

[0059] Among them, e current (t) represents the creep strain increment in the current interval. J is the definite integral operator in the time domain of the current k-th subinterval. αk (t-τ) is the constant-order creep compliance kernel function corresponding to the k-th subinterval, and t is the current specific observation time. k-1 Let τ be the starting time of the current k-th subinterval, and s be the internal variable of the integral. k (τ) represents the deviatoric stress component within the current k-th subinterval, ds k (τ) / dτ is the derivative of the deviatoric stress component with respect to time for the current k-th subinterval.

[0060] The method of strictly separating the historical long-term cumulative terms from the current transient incremental terms effectively avoids the unreasonable theoretical approximation of the full-time-domain stationary convolution kernel for variable-order aging materials, ensuring that the derivation of the entire mechanical forward modeling engine converges at the mathematical level and conforms to the superposition constitutive law of aging viscoelastic materials in a physical sense.

[0061] Step 205: Based on the deviatoric stress increment in the k-th sub-interval and the pre-configured elastic shear modulus, calculate the instantaneous elastic deviatoric strain. Then, superimpose the previous cumulative creep deviatoric strain, the current interval creep deviatoric strain increment, and the current interval instantaneous elastic deviatoric strain to obtain the deviatoric strain increment of each sub-interval.

[0062] Specifically, this step can also involve calculating the instantaneous elastic deviatoric strain based on the deviatoric stress increment in the k-th sub-interval and the pre-configured elastic shear modulus, superimposing the previous cumulative creep deviatoric strain, the current interval creep deviatoric strain increment, and the current interval instantaneous elastic deviatoric strain to obtain the total deviatoric strain at the end of the k-th sub-interval, and determining the deviatoric strain increment of the sub-interval based on the difference between the total deviatoric strain at the end of the k-th sub-interval and the total deviatoric strain at the end of the (k-1)-th sub-interval.

[0063] This is the final closed-loop stage of strain state update in the forward modeling engine. The system first calculates the purely elastic instantaneous deviatoric strain, independent of time accumulation effects, based on the deviatoric stress increment extracted in step 202 within the current k-th sub-interval of the material and the pre-configured elastic shear modulus, using the generalized Hooke's law. The specific calculation formula is as follows:

[0064] e elastic =Δs k / (2G');

[0065] Among them, e elastic For elastic instantaneous deviatoric strain, Δs k G' is the deviatoric stress increment in the k-th sub-interval determined in step 202, and G' is the pre-configured elastic shear modulus, which can be determined by the elastic modulus E' in the mechanical parameter set and the pre-configured Poisson's ratio ν through the relationship G'=E' / [2(1+ν)].

[0066] The cumulative creep deviatoric strain representing the long historical memory, the creep deviatoric strain increment representing the current interval excited by recent loads, and the instantaneous elastic deviatoric strain representing the transient geometric response are linearly added along the tensor dimension. The sum of these three components constitutes the complete deviatoric strain state evolution quantity of the material element at the end of the current sub-interval, i.e., the deviatoric strain increment of each sub-interval. During the assembly of the underlying finite element code, this deviatoric strain increment is recombined with the independent volumetric strain increment into a full strain tensor, which is mapped to the three-dimensional spatial displacement of each mesh node through the geometric compatibility equation. Finally, it is fed back as the output result to the upper-level optimization interface for comparison and verification of the inverted objective function. Furthermore, when dealing with the complex strongly constrained boundary state of the dam foundation rock mass, a pre-calibrated Poisson's ratio evolution matrix can be introduced to dynamically correct the lateral deformation coupling of the superimposed final strain deviatoric quantity, more accurately simulating the real three-dimensional nonlinear shear slip characteristics at the dam interface.

[0067] Example 3 describes the generalized mathematical boundary constraints of the core component order evolution model in the variable-order fractional-order Burgers constitutive relation, as well as three parallel specific driving evolution methods and non-monotonic alternatives for predetermined adverse conditions.

[0068] Step 301: The order evolution model is used to characterize the evolution of the fractional order of the Maxwell branch in the variable-order fractional Burgers constitutive relation as a function of the driving variables.

[0069] Specifically, this step defines the functional positioning of the order evolution model from a macroscopic physical mechanism perspective. In variable-order fractional-order rheology, the Maxwell branch describes the irreversible permanent viscous creep behavior of materials. As a mathematical operator characterizing viscous flow properties, the magnitude of the fractional order directly determines the strength of the material's creep activity. When the order approaches a unit value, the material exhibits strong viscous flow characteristics; when the order approaches zero, the material behavior degenerates into pure elasticity. The physical aging process of dam concrete gradually hardening and its microporous structure becoming denser under long-term loads and environmental influences is mathematically equivalent to a decay process with a continuously decreasing fractional order. The driving variables are the independent time scales or cumulative state measures that propel this decay process forward.

[0070] Step 302: The order evolution model satisfies the boundary conditions: when the driving variable is zero, the fractional order is equal to the initial order; when the driving variable approaches infinity, the fractional order asymptotically approaches the limit order; the initial order is greater than the limit order.

[0071] In this embodiment, this step defines a unified generalized mathematical boundary constraint for the various possible functional forms of order evolution. Specifically, when the dam has just completed impoundment and is at the start of its service life, the driving variable is equal to zero, the cement hydration reaction inside the dam body is relatively intense, the creep activity of the material is at its highest level, and the fractional order is anchored to a relatively high initial order. As the dam continues to operate, physical aging and hydration reactions gradually deplete. When the driving variable extends to infinity, the creep activity of the material will eventually stabilize at a residual level. At this point, the fractional order will smoothly converge to a lower limiting order. The strict inequality constraint that the initial order is greater than the limiting order conforms to the normal engineering aging law of the gradual degradation and stabilization of concrete material properties, ensuring that the subsequent inversion optimization algorithm does not diverge directionally within the optimization space.

[0072] Step 303, the driving variable is absolute service time; the order evolution model is defined as an exponential decay function:

[0073] α(t)=α ∞ +(α0-α ∞ )e (-βt) ;

[0074] Where α(t) is the fractional order corresponding to the absolute service time t, α0 is the initial order, and α ∞ β is the limiting order, and β is a positive parameter characterizing the aging rate.

[0075] Specifically, this step presents a basic and easily deployable parameterized implementation. The abstract driving variable is directly instantiated as the absolute service time measured in calendars. The system assumes that the dam material follows a macroscopically uniform aging degradation rate throughout its entire lifespan. In this basic implementation, the model introduces only three parameters that need to be identified through an inversion algorithm: the initial and limiting orders that determine the upper and lower bounds, and the positive aging rate parameter that controls the steepness of the curve's decline from high to low. This approach, while maintaining highly stable parameter identifiability, can capture the rapid early creep and stable convergence characteristics of most conventional hydraulic structures in their later stages.

[0076] In some alternative implementations, the aforementioned monotonically decreasing aging physics assumption may have limitations for deteriorating conditions such as alkali-aggregate reaction in dam concrete. In such cases, the order evolution model can be extended to a non-monotonic double-exponential substitution function. Specifically, the expansion of alkali-aggregate reaction products and the local propagation of microcracks may cause the creep activity of the material to increase in reverse at a certain service stage, exhibiting a non-monotonic rebound trend of first decreasing and then increasing on the evolution curve. Correspondingly, the substitution formula for such deteriorating conditions is:

[0077] α(t)=α ∞ +(α0-α ∞ )*exp(-β1*t)+α r *(1-exp(-β2*t));

[0078] Where α(t) is the fractional order corresponding to the absolute service time t, α ∞ Let α be the first limiting order, α0 be the initial order, and β1 be the aging decay rate in the first stage of service. r β1 represents the order rebound amplitude triggered by the second-stage degradation, and β2 represents the evolution rate of the second-stage degradation. By introducing this alternative formula, the system can be effectively endowed with the underlying capability to diagnose abnormal degradation mechanisms.

[0079] Step 304: The driving variable is the equivalent hydration age; obtain the measured temperature history in the time history of load and environmental effects; calculate the equivalent hydration age based on the measured temperature history using the Arrhenius temperature rate function, the calculation formula is: t eq (t)=∫0 t exp[(E a / R)*(1 / T ref -1 / T(τ))]dτ;

[0080] Among them, t eq (t) represents the equivalent hydration age corresponding to the absolute service time t, T(τ) represents the measured temperature history at time τ, and E aLet T be the apparent activation energy, R be the universal gas constant, and T be the apparent activation energy. ref The preset reference temperature is used; the order evolution model is defined as: α(t) = α ∞ +(α0-α ∞ )exp(-β*t eq (t));

[0081] Where α(t) is the fractional order, α0 is the initial order, and α ∞ β is the limiting order, and β is the decay rate constant of creep activity at the reference temperature.

[0082] Based on the above basic embodiments, this step provides a preferred implementation method that involves interdisciplinary coupling. The actual physical aging rate of dam concrete is strongly dependent on ambient temperature. The cement hydration rate is higher in summer high temperatures or in the deep, non-heat-dissipating areas of the dam than in winter or in the surface areas of the dam. Using a uniform calendar time would mask this critical spatial heterogeneity. This solution embeds the mature equivalent age theory from concrete materials science into a variable-order calculus operator. The system extracts long-term measured temperature history from existing temperature monitoring instruments on the dam and uses the Arrhenius formula to nonlinearly convert it into the virtual time required to achieve the same degree of curing under standard isothermal conditions, i.e., the equivalent hydration age.

[0083] Furthermore, in the specific engineering parameter configuration, the general gas constant R is fixed at 8.314 joules per mole Kelvin, and the preset reference temperature T is... ref The apparent activation energy E is typically set to 293 Kelvin. a Typical values ​​range from 30,000 to 45,000 joules per mole. These values ​​can be extracted directly from historical reports of cement adiabatic temperature rise tests or retrieved from a predefined database of ordinary silicate cement material properties.

[0084] Through this operation, different physical coordinate points in the finite element mesh will automatically evolve differentiated order decay trajectories based on their independent temperature histories, thus achieving a better time-scale nonlinear mapping between temperature stress effects and fractional-order constitutive models.

[0085] Step 305: The driving variable is normalized time; normalized time is calculated using absolute service time and the pre-configured total length of the observation period; the order evolution model is a semi-parametric expression based on Bernstein polynomial basis functions, the formula is:

[0086] ;

[0087] Where α(t) is the fractional order, α0 is the initial order, and α ∞ Let K be the limit order, K be the preset expansion order, ξ be the normalization time, and B be the limit order. k,K(ξ) are K-order Bernstein polynomial basis functions, c k This is the k-th Bernstein control coefficient. To ensure the physical plausibility of the monotonically decreasing fractional order, the Bernstein control coefficients satisfy the following constraints: c0=1; c0≥c1≥c2≥...≥c k ≥0.

[0088] In another specific embodiment, to address the potential rigidity in the shape preset of the aforementioned exponential decay function, this step provides a semi-parametric optimization scheme for shape identification driven by the depth of monitoring data. Since a dam may experience a brief creep plateau period or an S-shaped transition in its early service life, a purely exponential curve-driven forced fitting can easily lead to systematic biases. Therefore, the system first divides the absolute service time by the total length of the entire observation period to obtain a dimensionless normalized time, and then introduces a Bernstein polynomial basis, which possesses excellent endpoint preservation and convex hull properties, as a spatial generator.

[0089] Specifically, the semi-parametric expression discretizes the unknown continuous function of order evolution into a linear combination of a set of discrete control coefficients and corresponding basis functions. This construction replaces the highly complex functional constraints of determining the negation of derivatives for continuous functions by forcibly constraining a series of monotonically decreasing constant linear inequalities. In practical applications, the preset expansion order K is generally configured as 4, enabling the model to possess the expressive flexibility to fit various complex real aging paths, such as fast-then-slow, slow-then-fast, and even step-like pauses, with only a small increase in optimization dimensions.

[0090] Step 306: The parameters to be inverted also include a set of Bernstein control coefficients independent of the mechanical parameter set and the fractional evolution parameter set; in constructing the inversion objective function, a second-order difference smoothness penalty term is also constructed for the Bernstein control coefficient set; the second-order difference smoothness penalty term is added to the inversion objective function to suppress the oscillatory overfitting of the Bernstein control coefficients.

[0091] In this embodiment, this step is activated as a necessary defense mechanism for step 305. Because the semi-parametric model grants the inversion algorithm excessively high morphological freedom, when there are large areas of missing displacement monitoring data for a certain period of time or a low signal-to-noise ratio, the inversion program is prone to causing high-frequency fluctuations in the Bernstein control coefficients that do not conform to physical laws in order to accommodate local noise points. To eliminate this potential problem, an independent set of Bernstein control coefficients is extracted, and a second-order difference smoothness penalty term is custom-designed for it. The specific calculation formula is as follows:

[0092] ;

[0093] Among them, R smoothThis is a second-order difference smoothness penalty term. For a summation operator from subscript 1 to K-1, c k-1 Forward Bernstein control coefficients, c k c is the current central Bernstein control coefficient. k+1 These are the backward Bernstein control coefficients.

[0094] Understandably, this formula, based on the principle of approximation using the second derivative of discrete sequences, accurately measures the degree of curvature of the broken line formed by three adjacent control coefficients. When this penalty term is forcibly injected into the system's inversion objective function, the inversion engine, facing blind spots due to insufficient data constraints, will automatically force adjacent control coefficients to be arranged in an arithmetic sequence, causing the local order evolution curve to degenerate into the safest and smoothest linear interpolation form, reducing the risk of long-term prediction divergence induced by overfitting.

[0095] Example 4 details how to discretize the dam's service life in the most efficient way during the forward modeling physics engine construction phase using intelligent reconstruction technology of the temporal grid. The time step is dynamically adjusted using an equal error allocation criterion, and a fallback constraint is introduced in extreme cases to minimize local truncation errors while ensuring the stability of the underlying finite element numerical solution.

[0096] The dam's service life is discretized into multiple sub-intervals, including using an equal error allocation criterion based on the order of change rate to determine the endpoints of non-uniform sub-intervals, such as... Figure 3 As shown, the specific steps include the following:

[0097] Step 401: Calculate the absolute value of the rate of change of the order evolution model over the dam's service life.

[0098] Specifically, this step is the mathematical prerequisite for achieving adaptive control of constant-order approximation error. In the early stages of normal dam service, due to the rapid cement hydration reaction and microstructural densification process, the creep activity of concrete materials decays rapidly, causing a sharp drop in the fractional-order metric of time, resulting in a large absolute value of its rate of change. If a conventional fixed-step partitioning strategy is used, a high level of piecewise constant-order approximation local truncation error will inevitably accumulate during periods of drastic order changes. To accurately quantify the degree of change in material properties at different service stages, a first-order differential operation is performed on the predetermined order evolution model with respect to absolute service time, and the absolute value is taken to obtain the absolute value of the rate of change characterizing the dynamic intensity of aging. This continuous rate signal constitutes the driving force for the subsequent construction of a non-uniform time-discrete grid.

[0099] Step 402: Integrate the absolute value of the rate of change over the absolute service time to calculate the order arc length function and the total order arc length corresponding to the dam's service cycle.

[0100] In this embodiment, the system utilizes integral transformation to map the rate of change index obtained in the previous step into a geometric metric space with cumulative characteristics. By performing a definite integral on the absolute value of the rate of change starting from the zero point of service, a monotonically increasing order arc length function with time is obtained. This function establishes a novel virtual time scale: during periods of dramatic order changes in real physical time, the virtual order arc length increases rapidly; while in the later stages of service where the order tends to level off in real physical time, the order arc length almost stops increasing. Extending the upper limit of the integral to the end of the entire dam service cycle allows for the calculation of the total order arc length. The specific integral formula is expressed as follows:

[0101] ;

[0102] Where s(t) is the calculated order arc length function, ∫0 t Let v be a definite integral operator from the origin of time to the current absolute service time t. rate (τ) represents the first derivative of the order evolution model corresponding to the integral variable τ with respect to time, and dτ is the integrand differential. As an example, when the system adopts the basic exponential decay model, the above integral has an exact analytical solution, and its order arc length function can be directly expressed as the difference between the initial order and the limiting order multiplied by the exponential decay complement term controlled by the decay rate constant.

[0103] Step 403: Divide the total order arc length equally according to the preset total number of sub-intervals to obtain multiple arc length division points.

[0104] Specifically, this step performs uniform mesh partitioning within the virtual geometric space constructed by the order arc length function. To achieve a balanced distribution of global computational errors, the system requires that the magnitude of local truncation errors caused by the constant-order approximation be equal in each discrete time step. This error distribution criterion is mathematically rigorously and equivalently proven to be a uniform division of the total order arc length. The system obtains the total number of sub-intervals pre-configured by the engineer or inversion framework, and uses this number as the partitioning base to divide the total order arc length into several line segments of equal length. The partitioning nodes of each line segment constitute multiple arc length partitioning points. The calculation formula is:

[0105] s j =j*(S total / N');

[0106] Among them, s j S represents the value of the j-th arc length division point, where j is the index number of the preset sub-intervals, increasing from zero. total Let N' be the total order arc length, and N' be the preset total number of subintervals. Through this ingenious transformation domain division operation, the system effectively avoids the difficulty of directly performing complex nonlinear error calculations on the physical time axis.

[0107] Step 404: Using the inverse function of the order arc length function, multiple arc length division points are mapped back to the physical time coordinate to obtain the endpoints of the non-uniform sub-intervals of each sub-interval.

[0108] In this embodiment, the system reverse-engineers the virtual arc length nodes configured with equal errors into a real time scale usable in engineering. Utilizing the strictly monotonically increasing property of the order arc length function, the algorithm constructs its unique mathematical inverse function and substitutes each arc length point obtained in the previous step as an input parameter into this inverse function for solution. The resulting physical time scalar sequence after the reverse mapping represents the endpoints of the non-uniform sub-intervals. Under this adaptive allocation mechanism, in the early stages, due to the need for more physical time to accumulate sufficient arc length increments, the generated endpoints are relatively dense on the time axis; while in the later stages, even small order fluctuations correspond to long real times, and the spacing between the endpoints is automatically and significantly increased.

[0109] Taking the exponential decay pattern as an example, the inverse function mapping formula is derived from steps 402-403 as follows:

[0110] t j =(-1 / β)*ln(1-j / N'(1-exp(-βT)));

[0111] Among them, t j Let be the endpoint of the j-th non-uniform sub-interval mapped back to physical time coordinates, and β be a positive parameter controlling the aging rate. Let be the natural logarithm operator, j be the step index corresponding to the arc length division point, N' be the preset total number of subintervals, and T be the end time of the dam's service life. When j=N', t N' =T, ensuring a complete service life of the dam with discrete time-domain coverage.

[0112] In some alternative implementations, when βT is much greater than 1, i.e., the dam's service life is much longer than the characteristic time of material aging 1 / β, the above formula can be simplified to:

[0113] t j =(-1 / β)*ln(1-j / N').

[0114] Step 405: Calculate the time span between the endpoints of two adjacent non-uniform sub-intervals to obtain the width of each sub-interval.

[0115] Specifically, as a preparatory step for introducing a numerical stability fallback, the system needs to examine the geometric span of the newly generated adaptive mesh in physical space. The algorithm traverses all endpoints of the non-uniform subintervals, calculating the absolute difference between the time coordinates of the next endpoint and the time coordinates of the previous endpoint. A series of differences constitute the actual duration of each subinterval in the physical time domain, i.e., the width of each subinterval.

[0116] Step 406: When the width of any sub-interval is greater than the pre-configured upper limit of the interval width, the width of the corresponding super-wide sub-interval is divided into two equal parts, and an endpoint is added in the corresponding interval.

[0117] In this embodiment, this step is used to prevent excessive divergence of the adaptive algorithm under extreme conditions. As mentioned earlier, in the very later stages of a dam's service life, the rate of change of the order approaches zero, and the equal error allocation criterion tends to assign an excessively large time step spanning decades. An excessively large time step may cause the finite element solver to have difficulty converging in nonlinear, unsteady load increment iterations due to a deterioration in the condition number of the Jacobian matrix. To eliminate this potential engineering hazard, a mandatory rigid truncation rule is introduced. The algorithm compares the width of each sub-interval with the pre-configured upper limit of the interval width. If it finds a situation exceeding this safety threshold, it forcibly strips the adaptive exemption of that interval, takes the arithmetic mean of the start and end endpoints of the interval as the insertion coordinates, and performs bisection processing to forcibly add new endpoints. Furthermore, the pre-configured upper limit of the interval width can be flexibly set according to the specific service life of the project; for example, it can be strictly limited to one-third of the total service time.

[0118] Step 407: Incorporate the newly added endpoints into the original sequence to obtain the updated non-uniform sub-interval endpoints.

[0119] Specifically, after traversing and locally segmenting all abnormally wide intervals, the system merges all newly generated endpoints with the original adaptive endpoint set. The merged time series is then rearranged in ascending order to generate the updated non-uniform sub-interval endpoint set ultimately used for the underlying integration engine. This updated grid system maintains the accuracy of maximizing the capture of early, severe creep while also possessing a safety margin against long-term integration divergence.

[0120] Based on the above embodiments, to accommodate scenarios in engineering practice where computing power is limited or the inversion accuracy requirements are not high, an alternative time-domain discretization implementation is provided. In this alternative, the system can directly skip all complex nonlinear calculations related to the rate of change of order and arc length integral, degenerating into the most basic uniform interval division strategy. Specifically, the total physical duration of the dam's service life is directly divided by the preset total number of sub-intervals to obtain a fixed constant step size; starting from the absolute service time zero point, this constant step size is used as the accumulation base to successively superimpose, quickly generating a sequence of equidistant sub-interval endpoints. It is understandable that this uniform discretization alternative scheme greatly simplifies the pre-calculation logic, but at the cost of artificially setting a total number of sub-intervals several times larger than the adaptive strategy to achieve the same early-stage approximation accuracy as the adaptive mesh, which will lead to a large amount of redundancy in the finite element integration calculations in the later stages of the dam's service life.

[0121] Example 5 details an adaptive grouping regularization optimization engine based on dual normalization, constructed to address the challenges of dimensional inconsistency and extremely poor sensitivity caused by joint inversion of multiple heterogeneous parameters. It demonstrates the complete mathematical derivation process from parameter normalization to closed-loop parameter update via grouping regularization, and supplements the physical meaning of the underlying switching mechanism of the optimization algorithm.

[0122] Step 501, before constructing the inversion objective function, also includes dimensionless normalization of each parameter to be inverted in the mechanical parameter set and the fractional-order evolution parameter set, such as... Figure 4 As shown.

[0123] Specifically, this step is a fundamental preparation for eliminating the dimensional curse of mixed physical quantities in inversion. In variable-order fractional-order creep inversion problems, the mechanical parameter set includes traditional mechanical parameters such as the elastic modulus, whose numerical magnitudes can reach tens of thousands of megapascals and have a pressure dimension. In contrast, the fractional-order evolution parameter set includes parameters such as the initial order and the limit order, which are purely dimensionless values ​​within the open interval 0 to 1. If these two sets of parameters with drastically different physical properties are directly placed in the same error functional based on Euclidean distance for regularization penalty calculation, it will lead to inconsistent dimensions in the regularization term, causing mathematical non-rigor and severe oscillations and divergences in the optimization process. Therefore, the system must perform a thorough dimensionless normalization and dimensionality reduction mapping on all variables to be inverted before constructing the objective function.

[0124] Step 502: Obtain the pre-configured reference scale for each parameter to be inverted.

[0125] In this embodiment, the reference scale is the baseline anchor point for performing dimensionless mapping. For most parameters, the absolute value of the prior mean set based on engineering experience or geological survey data is directly extracted as the reference scale. However, for some evolutionary parameters with a prior mean of zero, the system intelligently switches the baseline and extracts half the width of the physically permissible range of the parameter as the reference scale. This pre-configuration mechanism ensures that each physical variable is assigned a normalized denominator that strictly matches its true magnitude.

[0126] Step 503: Divide each parameter to be inverted by its corresponding reference scale to obtain the corresponding dimensionless parameters.

[0127] Specifically, this operation transforms the parameter space into a pure mathematical optimization space. The transformation formula is:

[0128] p bar_i =p i / p ref_i ;

[0129] Where, p bar_i p is the dimensionless parameter corresponding to the i-th parameter to be inverted.i For the i-th parameter to be inverted, carrying the original physical dimensions, p ref_i is the reference scale corresponding to the i-th parameter to be inverted.

[0130] After the above division operation, all parameters are stripped of their original physical dimensions, and in the normalized parameter space, the magnitude of the numerical fluctuation of each parameter is compressed to a similar benchmark level, laying an unbiased mathematical foundation for subsequent fairness sensitivity assessment and penalty weight allocation.

[0131] Step 504: Calculate the partial derivatives of the displacement of each measuring point with respect to each parameter to be inverted using a variable-order fractional forward model, and multiply the partial derivatives by their respective reference scales to obtain the normalized sensitivity vectors corresponding to each parameter to be inverted.

[0132] In this embodiment, this step quantitatively evaluates the relative control weight of each dimensionless parameter on the macroscopic displacement output. The finite difference method is used to apply a small perturbation to the variable-order fractional-order forward model to approximate the partial derivatives. To balance the truncation error of the difference method with the rounding error of floating-point calculations, the perturbation step size of the finite difference method is pre-configured to a sufficiently small threshold, such as a parameter reference scale of one ten-thousandth. After obtaining the original partial derivatives, they must be multiplied by their respective reference scales to complete the chain rule mapping of the sensitivity vector, as shown in the formula:

[0133] ;

[0134] in, Let p be the normalized sensitivity vector corresponding to the i-th parameter to be inverted. ref_i Let i be the reference scale corresponding to the i-th parameter to be inverted. Calculate the minute changes in displacement at each measuring point. Let be the small perturbation of the i-th parameter to be inverted. Through this operation, the dimensions of the normalized sensitivity vector are unified into pure displacement dimensions, clearly characterizing the influence of each parameter on the displacement field of the dam measuring point when a relative fluctuation of the same proportion occurs, thereby achieving direct comparability of sensitivity across physical quantity categories.

[0135] Step 505: In constructing the inversion objective function, dimensionless grouped weighted regularization matrices corresponding to the mechanical parameter set and the fractional evolution parameter set are constructed based on the normalized sensitivity vector.

[0136] Specifically, using the relative sensitivity information obtained in the previous step, a differentiated parameter adaptive penalty constraint system is constructed. Since the order evolution parameters indirectly affect displacement through a highly nonlinear convolution path, their sensitivity to data is typically much lower than that of the mechanical parameters that directly affect the stiffness matrix. To prevent unreasonable parameter drift due to insufficient data constraints caused by low-sensitivity parameters, weighting matrices must be designed independently for different parameter groups.

[0137] Step 506: Calculate the L2 norm of the normalized sensitivity vector corresponding to each parameter to be inverted in the mechanical parameter set and the fractional evolution parameter set.

[0138] In this embodiment, the L2 operator is used to scalarize and reduce the dimensionality of the vector. The algorithm traverses the set of mechanical parameters and the set of fractional evolution parameters, and calculates the square root of the sum of squares of all elements within each normalized sensitivity vector to obtain a comprehensive scalar eigenvalue representing the global influence intensity of the parameter, i.e., the L2 norm. The larger the L2 norm, the easier it is for the parameter to be identified by measured deformation data.

[0139] Step 507: Extract the maximum value of the L2 norm from the mechanical parameter group and the fractional evolution parameter group respectively, and use it as the sensitivity norm reference scale for each parameter group.

[0140] Specifically, the system searches for the most sensitive leading parameter within each group. The algorithm executes a maximum value comparison instruction in both the mechanical parameter group and the fractional evolution parameter group, defining the maximum value of the selected L2 norm as the sensitivity norm reference scale. This maximum value represents the parameter feature with the strongest voice in the data fitting within the current group, and serves as a standard anchor point for measuring the relative sensitivity deficiency of other parameters within the same group.

[0141] Step 508: Divide the square of the sensitivity norm reference scale corresponding to each parameter group by the square of the L2 norm corresponding to each parameter to be inverted in that group to obtain the dimensionless weighting coefficients corresponding to each parameter to be inverted.

[0142] In this embodiment, this is the core feedback control logic in the entire regularization mechanism. The calculation formula is:

[0143] ;

[0144] Among them, w i S is the dimensionless weighting coefficient corresponding to the i-th parameter to be inverted. g This serves as the sensitivity norm reference scale for the group to which this parameter belongs. Let L be the L2 norm of the normalized sensitivity vector corresponding to the i-th parameter to be inverted.

[0145] Since the dimensions of the numerator and denominator are both the square of the displacement, the dimensionless weighting coefficient obtained after division is a purely dimensionless value. This formula constructs a reciprocal decay mechanism: the parameter with the highest sensitivity within the group will receive the minimum regularization weight equal to the unit value, and is allowed to freely deviate from prior knowledge to follow the observed data signal; conversely, the lower the sensitivity of the parameter, the smaller its denominator, and the larger the dimensionless weighting coefficient it will receive, and will be subject to a more stringent prior mean anchoring constraint in the inversion. This strategy of automatically allocating constraint strength based on the parameter's own identifiability effectively solves the problems of over-constraint of strong parameters and under-constraint of weak parameters that are easily caused by traditional uniform constant regularization methods.

[0146] Step 509: Using the dimensionless weight coefficients corresponding to each parameter to be inverted as diagonal elements, assemble the dimensionless grouping weighted regularization matrices corresponding to the mechanical parameter group and the fractional evolution parameter group, respectively.

[0147] Specifically, the discrete weight scalars generated in the previous step are integrated into a linear algebra matrix structure. The algorithm creates two initial zero matrices, and then sequentially fills the main diagonal with the dimensionless weight coefficients of each parameter within the mechanical parameter group to form the dimensionless grouped weighted regularization matrix of the first parameter group; similarly, the dimensionless weight coefficients of each parameter within the fractional evolution parameter group are filled into the main diagonal of another zero matrix to form the matrix of the second parameter group. The off-diagonal elements of the matrix are strictly zero, representing the assumption that the parameters are independent and uncorrelated in their prior distributions.

[0148] Step 510: During the process of alternately optimizing the mechanical parameter set and the fractional evolution parameter set, the regularization parameters of each parameter set are adaptively updated using the relative deviation between the data fitting state and the prior constraint state.

[0149] In this embodiment, this step provides a dynamic control valve for the entire optimization loop. As explained in the previous framework planning, the alternating optimization strategy essentially comprises two decoupled loops, inner and outer. Specifically, in the early stages of iterative optimization, the system preferentially calls a genetic algorithm as the global search engine. The genetic algorithm introduces an elite retention mechanism and adaptive crossover mutation probability to widely disperse the population in the vast multidimensional non-convex parameter space, giving the system the physical detection capability to escape the trap of local pseudo-convergent minima.

[0150] Specifically, the key hyperparameters of the genetic algorithm are configured as follows: the population size is set to 10 to 20 times the dimension of the parameters to be inverted; for example, when the total dimension of the parameters to be inverted is 10, the population size is configured to 100 to 200 individuals. The selection strategy uses a tournament selection method, with a tournament size of 3. The initial value of the crossover probability is set to 0.8 to 0.9, and the initial value of the mutation probability is set to 0.01 to 0.05. The elite retention ratio is set to 5% to 10% of the population size. The fitness function is directly the negative value of the inverted objective function. The maximum number of iterations is set to 500 to 2000 generations. When the standard deviation of the fitness of individuals in the population is lower than the preset convergence threshold within 50 consecutive generations, for example, one-thousandth of the absolute value of the objective function, it is determined that the basin region where the global optimum is located has been locked.

[0151] When the fitness variance of individuals in the genetic population falls below a preset convergence threshold, the system determines that it has locked the basin region where the global optimum is located, immediately shuts down the genetic algorithm, and smoothly switches to the quasi-Newton algorithm for relay execution. The quasi-Newton algorithm utilizes the gradient of the objective function and the Hessian approximation information of the Hessian matrix to perform fine and fast convergence optimization within a local region. After each round of alternating inner and outer optimization loops, a dynamic shift operation is performed on the grouping regularization parameter that adjusts the weight of the data fitting term and the prior constraint term based on the current fitting state.

[0152] Step 511: Calculate the displacement of each measuring point corresponding to the dimensionless parameters under the current iteration using the variable-order fractional forward model, and subtract the calculated displacement of each measuring point from the measured displacement time series of multiple measuring points to obtain the displacement residual vector.

[0153] Specifically, the forward modeling engine is invoked to predict the theoretical deformation response of the dam based on the newly identified dimensionless parameter array. By comparing the data with actual monitoring data through vector subtraction, an absolute error column vector reflecting the current model's prediction deviation is extracted, namely the displacement residual vector. This residual vector intuitively reflects the degree to which the current parameter combination is insufficient in reproducing the real physical process.

[0154] Step 512: Calculate the dimensionless data fitting term using the displacement residual vector and the pre-configured data error covariance matrix.

[0155] In this embodiment, a statistical perspective is further introduced to weight the error measurement. Considering the varying service environments of different sensor nodes at the dam site, leading to drastically different measurement accuracies, the system reads a pre-configured data error covariance matrix. In most engineering scenarios, assuming that the errors at each observation point are independent, this matrix can be degenerated into a diagonal matrix containing the variance of the displacement observation errors at each measuring point. The calculation formula is:

[0156] F d =r n T*C d -1 *r n ;

[0157] Among them, F d For dimensionless data fitting term, r n Let r be the displacement residual vector. n T C is the transpose of the displacement residual vector. d Let C be the data error covariance matrix. d -1 It is its inverse matrix. Since the data error covariance matrix has the dimension of the square of the displacement, through the quadratic form operation after weighting the inverse matrix, the dimensionless data fitting term loses its physical dimension and becomes a statistical scalar to measure the overall data fidelity and purity of the model.

[0158] Step 513: Using the dimensionless parameters under the current iteration, the pre-configured prior mean of the dimensionless parameters, and the corresponding dimensionless grouping weighted regularization matrix, calculate the dimensionless parameter deviation quadratic form of the mechanical parameter group and the fractional evolution parameter group, respectively.

[0159] Specifically, this step involves penalizing the parameter search span. The algorithm extracts the deviation distance vectors between the current iteration parameters and the pre-configured initial prior assumptions for both the mechanical parameter set and the evolution parameter set. Then, using the dimensionless grouped weighted regularization matrix generated in step 509 as the intermediate kernel metric matrix, it calculates the deviation metric in the Mahalanobis distance formula. The calculation formula is as follows:

[0160] F prior_g =(p bar_g n -p bar_g 0 ) T *W g *(p bar_g n -p bar_g 0 );

[0161] Among them, F prior_g For the dimensionless parametric deviation quadratic form of the g-th parameter group, p bar_g n p is the dimensionless parameter vector for the current iteration. bar_g 0 For the pre-configured dimensionless parameter prior mean vector, (p bar_g n -p bar_g 0 ) T W is the transpose of the parameter deviation vector. gLet g be the dimensionless grouping weighted regularization matrix corresponding to the g-th parameter group. The result of this quadratic operation is also a completely dimensionless positive real scalar, representing the degree of risk of each parameter group deviating from the a priori engineering knowledge zone in the dimensionless weighted space.

[0162] Step 514: Using the ratio of the dimensionless data fitting term to the dimensionless parameter deviation quadratic form corresponding to each parameter group, calculate the updated grouping regularization parameters of the mechanical parameter group and the fractional evolution parameter group in the next iteration.

[0163] In other words, this step can also be to use the dimensionless data fitting term, the ratio of the number of dimensions of each parameter group to the total amount of observed data, and the quadratic form of the dimensionless parameter deviation to calculate the grouping regularization parameters of the mechanical parameter group and the fractional evolution parameter group after the update in the next iteration.

[0164] In this embodiment, this step completes the mathematical closed loop of the adaptive update logic. A seesaw-like balancing mechanism is established: the ratio of two dimensionless quadratic forms is used as the adjustment multiplier. The specific calculation formula is as follows:

[0165] λ g (n+1) =((m g / M)*F d ) / F prior_g ;

[0166] Where, λ g (n+1) m represents the updated grouping regularization parameter for the next iteration. g Let F be the number of dimensions of the g-th parameter group, M be the total amount of observed data, and F be the number of dimensions of the g-th parameter group. d For dimensionless data fitting term, F prior_g Let g be the dimensionless parametric deviation quadratic form corresponding to the g-th parameter group.

[0167] Driven by this update rule, if the current dimensionless data fitting term is too large, it indicates that the current model lacks the accuracy of tracking the measured data. The algorithm will automatically increase the numerator of the ratio, reduce the regularization constraint strength, and allow the parameters to have a greater degree of adjustment redundancy to fit the data. Conversely, if the calculated dimensionless parameter deviation quadratic form is large, it indicates that the parameters have deviated too much from the prior reasonable range and there is a risk of distortion. The surge in the denominator will lead to a smaller updated regularization parameter, which will enhance the centripetal pull.

[0168] In addition, the dimensional adjustment coefficient (m) introduced in the formula g / M) eliminates the artificial statistical bias caused by the fact that the dimension of the mechanical parameter group is usually much higher than that of the fractional evolution parameter group, and ensures that the high-dimensional group and the low-dimensional group are always in an equal dialogue system in the dynamic allocation of regularization tightness.

[0169] Example 6 details how, after the parametric inversion optimization converges, the temporal characteristics of the residual sequence are used to adaptively correct the order evolution model using nonparametric methods, and how strict physical boundary constraints ensure the rationality of the correction results. This example constitutes a two-level identification framework consisting of parametric baseline identification and nonparametric local correction, effectively compensating for potential structural biases in the preset parametric function form.

[0170] Step 601, before outputting the final inversion parameter set, also includes a step of performing order evolution nonparametric correction using residual time series diagnostics.

[0171] Specifically, this step establishes a post-processing closed-loop mechanism for model self-checking and repair. In conventional inversion processes, once the global optimization algorithm determines that the parameter set has reached its optimum, the system usually outputs the results directly. However, if the preset order evolution model function form itself differs from the actual material aging law of the dam, even if parameter optimization reaches its limit, there will still be a systemic structural mismatch between the model calculation results and the measured data. To overcome the limitations of this single parameterization form, this embodiment introduces a local correction module based on residual time series before outputting the final results. By treating the values ​​of the order evolution function in each sub-interval as independent adjustable quantities, the residual information of the data itself is used to refine the theoretical model and capture the fine features that deviate from the main trend.

[0172] Step 602: Calculate the displacement of each measuring point corresponding to the final inversion parameter set using the variable-order fractional forward model. Subtract the calculated displacement of each measuring point from the measured displacement time series of multiple measuring points to obtain the final displacement residual vector.

[0173] In this embodiment, this step extracts the current optimal fitting deviation of the model. The final inversion parameter set obtained after alternating global and local optimization is substituted back into the forward modeling physics engine to generate a theoretically predicted deformation time history covering the entire observation period. The system extracts the actual measured displacement time series of multiple measurement points point by point, and subtracts the aforementioned theoretically predicted deformation time history, thereby extracting the pure deviation sequence. The deviation data of all spatial measurement points and time nodes are flattened and arranged in chronological order in one dimension, thus assembling the final displacement residual vector. This vector contains all physical dynamic signals that were not successfully explained by the current parameterized model.

[0174] Step 603: Calculate the residual moving average sequence and Durbin-Watson statistic of the final displacement residual vector, and determine whether there is structural mismatch based on the Durbin-Watson statistic and the residual moving average sequence.

[0175] Specifically, the system needs to perform professional time-series feature mining on the extracted error signals. The algorithm selects a preset time window length, performs smoothing filtering on the final displacement residual vector, and calculates the mean error within each window to generate a smooth residual moving average sequence. This sequence can effectively filter out high-frequency random observation noise and expose the low-frequency drift trend of the error during long-term service. At the same time, the first-order autocorrelation test index, namely the Durbin-Watson statistic, is calculated. This statistic rigorously detects whether there is a temporal dependency in the error sequence by quantifying the ratio of the sum of squares of the differences between residuals at two adjacent time steps to the sum of squares of residuals.

[0176] Step 604, the specific criteria for determining the existence of structural mismatch are as follows: when the Durbin-Watson statistic is less than a preset lower threshold or greater than a preset upper threshold, or when the absolute mean of the residual moving average sequence exceeds a preset trend threshold, a structural mismatch is determined to exist.

[0177] In this embodiment, this step establishes a quantitative judgment rule for triggering the nonparametric correction mechanism. In classical statistical theory, the ideal pre-set threshold for no autocorrelation is usually set to a value of 2. If the calculated Durbin-Watson statistic is significantly lower than or significantly higher than a value of 2, it proves that there is obvious positive or negative autocorrelation within the residual sequence, respectively. On the other hand, if the residual moving average sequence does not oscillate randomly around the zero axis with white noise, but exhibits a systematic time trend such as a sustained positive bias in the early stage and a sustained negative bias in the later stage, it is also a clear signal of unreasonable fitting. When either of the above two time series anomalies is triggered, the system makes a judgment, determining that the current parametric function form fails to accurately fit the true nonlinear decay trajectory, judging that there is a structural mismatch, and immediately activating the subsequent correction calculation procedure.

[0178] Step 605: When structural mismatch is determined to exist, calculate the sensitivity matrix of the order of the fractional order of each sub-interval in the variable-order fractional forward model to the order of the displacement calculated at each measuring point.

[0179] Specifically, this step provides gradient guidance information for calculating the correction. The system decouples the parameterized functions between the orders of each sub-interval, treating the constant-order value of each sub-interval within the piecewise constant-order framework as an independent variable. For each independent sub-interval, the sensitivity of the displacement response of all observation points to a unit small perturbation of the fractional-order order of that sub-interval is calculated using the finite difference method or analytical derivation, forming a column vector. The sensitivity column vectors of all sub-intervals are then horizontally concatenated in chronological order to construct a complete order sensitivity matrix. This matrix serves as a linear transformation operator that inversely maps the displacement space residual to the order space correction.

[0180] Step 606: Using the final displacement residual vector, the order sensitivity matrix, the pre-configured second-order difference matrix, and the pre-configured smoothness regularization parameter, construct and solve the regularized least squares problem to obtain the order correction vector.

[0181] In this embodiment, this step provides a pure algebraic closed-form solution for nonparametric correction. To prevent severe high-frequency oscillations in the correction amount caused by relying on residuals for reverse derivation, a pre-configured second-order difference matrix is ​​introduced to apply a smoothing penalty to the order correction amount of adjacent sub-intervals, and a pre-configured smoothness regularization parameter is used to balance the weight between data approximation and the smoothness of the correction curve. Furthermore, the smoothness regularization parameter can be automatically determined by the generalized cross-validation criterion. The problem of finding the optimal correction amount is transformed into a typical regularized least squares mathematical optimization problem. The closed-form solution calculation formula for this problem is:

[0182] δ α_opt =(G T *G+μ*D2 T *D2) -1 *G T *r;

[0183] Where, δ α_opt To obtain the optimal order correction vector, G is the order sensitivity matrix. T Let D2 be its transpose matrix, μ be a pre-configured smoothness regularization parameter, and D2 be a pre-configured second-order difference matrix. T Its transpose matrix, (G T *G+μ*D2 T *D2) -1 The expression represents the inversion of the combined matrix within the parentheses, where r is the final displacement residual vector. Using the above linear algebraic formula, the system can calculate the optimal order compensation value required for each time sub-interval in a single, non-iterative operation, improving the algorithm's execution efficiency and robustness.

[0184] Step 607: Superimpose the order correction vector onto the piecewise fractional order sequence corresponding to the final inversion parameter set, and update the final inversion parameter set.

[0185] Specifically, the theoretical parameterized order sequence calculated from the original parameterized evolution equation at the center points of each sub-interval is extracted. The order correction vectors calculated in the previous step are then arithmetically added term by term according to the corresponding interval indices. The linear superposition operation seamlessly injects the data-driven fine-tuning signal into the baseline main trend, generating a new corrected piecewise order sequence that integrates macroscopic decay patterns and microscopic real fluctuations. The system then uses this corrected sequence to replace the original pure parameterized configuration, completing the substantial update of the final inversion parameter set.

[0186] Step 608, after superimposing the order correction vector onto the piecewise fractional order sequence corresponding to the final inversion parameter set and updating the final inversion parameter set, also includes a physical boundary constraint step.

[0187] In this embodiment, this step constitutes a safety valve for the entire correction process. Since the regularized least squares solution process is a purely unconstrained convex optimization mathematical process, the theoretical correction values ​​derived from it, when superimposed on the original order, may cause the final order value in some intervals to exceed the theoretical limits of mechanics of materials. To ensure that the updated parameter set can be directly used for subsequent safety assessments and finite element predictions without causing the underlying calculation program to crash, a physical boundary constraint step must be forcibly inserted before output.

[0188] Step 609: Determine whether the values ​​of each interval of the segmented fractional-order sequence in the updated final inversion parameter set exceed the pre-configured physical allowable open interval.

[0189] Specifically, the order physical meaning of dam concrete under fractional-order rheological theory requires that its viscoelastic properties lie between pure elasticity and pure viscosity. The pre-configured physical allowable open interval is strictly set as a set of real numbers greater than 0 and less than 1. The system iterates through each constant value in the updated piecewise fractional-order sequence one by one, and uses a logical discrimination rule to verify whether the value is lower than the lower limit of the interval or higher than the upper limit of the interval.

[0190] Step 610: If any interval value exceeds the physically allowed open interval, the excess part is truncated to the corresponding interval boundary value to obtain the final inversion parameter set updated with boundary constraints.

[0191] In this embodiment, when the system detects that the order value of a certain sub-interval has exceeded the limit, it will immediately trigger a forced pruning procedure. For example, if the corrected order of a certain interval decreases to a negative value, the system will forcibly reset it to a small positive number; if the corrected order of a certain interval soars beyond the unit 1, it will be truncated to a slightly less than 1 limit positive number. After a comprehensive one-dimensional scan and forced pruning, all illegal abnormal parameters are cleared, and the final output to the outside is the updated inversion parameter set after boundary constraints, which meets the physical self-consistency requirements with high confidence. Furthermore, as an optional supporting scheme, after completing the correction and pruning of the order sequence, the latest order sequence can be fixed and a simplified regression inversion can be re-executed only for the mechanical parameter set, thereby prompting the elastic modulus and other mechanical characteristics to adaptively adapt to the adjusted local order fluctuations, achieving deep microscopic coordination between parameter sets.

[0192] If no structural mismatch is determined, the final inversion parameter set remains unchanged.

[0193] Example 7 describes a dual-benchmark comparison verification method designed to verify the superiority and physical rationality of the proposed variable-order fractional-order model in long-term extrapolation prediction. It also provides a set of typical numerical examples of dam creep parameter inversion, showcasing the specific implementation details and parameter magnitude characteristics of the underlying algorithm.

[0194] Step 701: Output the final inversion parameter set and predict creep characteristics.

[0195] In this embodiment, this step serves as the final output and engineering application interface of the parameter identification process. Before the system outputs the final inversion parameter set, rigorous numerical simulation and prediction comparison are required to ensure that the parameter set not only fits well within the known historical observation data range but also possesses non-divergent stability within the future prediction range. To make the technical details more specific, a set of virtual numerical implementation examples of a high dam model is introduced. The system sets the prior mean of the elastic modulus of the dam concrete to 25,000 MPa, the prior mean of the Maxwell branch viscosity coefficient to a fractional power of 10,000 MPa multiplied by days, the prior mean of the initial fractional order to 0.6, and the limiting order to 0.2. The first 5 years after the dam impoundment are selected as the historical observation period, and the aforementioned dual dimensionless and adaptive grouping regularization optimization engine is used for inversion. Assuming the final inversion parameter set from the final converged output contains an optimized elastic modulus of 26,500 MPa, an initial order of 0.55, a limiting order of 0.15, and a decay rate parameter of 0.01 ± 1 per day, outputting this parameter set lays the numerical foundation for subsequent long-term extrapolation verification.

[0196] To clearly demonstrate the specific process of the inversion calculation, a simplified numerical example is given below. Assume the dam mesh model contains 5000 elements, three characteristic measuring points are selected at the dam crest, and the observation period is from day 0 to day 1825 (5 years) after impoundment. The initial parameters are set as follows: the prior value of the elastic modulus is 25000 MPa, and the reference scale is 25000 MPa; the prior value of the initial order α0 is 0.6, and the reference scale is 0.6; the limiting order α... ∞ The prior value of is 0.2, and the reference scale is 0.2; the prior value of the decay rate β is 0.005 to the power of negative one per day, and the reference scale is 0.005.

[0197] After dimensionless normalization, the initial dimensionless value of all parameters is 1.0. After the first round of forward modeling, the L2 norm of the displacement residual vector at 1825 observation times across three measuring points is 12.5 mm. The calculated normalized sensitivity L2 norm of the elastic modulus is 8.3 mm, while that of α0 is 0.15 mm, reflecting a 55-fold difference in sensitivity. Therefore, the dimensionless weighting coefficient of the elastic modulus in the mechanical parameter set is 1.0, while the dimensionless weighting coefficient of α0 in the fractional-order evolution parameter set is approximately (0.15 / 0.15). 2 =1.0, the largest within the group, and the weighting coefficient of β is larger due to lower sensitivity.

[0198] After 350 generations of global search using a genetic algorithm, the fitness variance converged. The algorithm was then switched to a local fine-tuning algorithm (L-BFGS quasi-Newton algorithm) with limited memory for the Broyden-Fletcher-Goldfarb-Shannon algorithm. After three rounds of alternating outer-layer iterations, the relative changes of each parameter were all below 0.1%, and the objective function decreased by less than 0.01%, indicating convergence. The final output elastic modulus was 26500 MPa, the initial order was 0.55, the limiting order was 0.15, and the decay rate was 0.01 m / s² -1. At this point, the L2 norm of the displacement residual vector decreased to 1.2 mm, a 90.4% reduction from the initial value.

[0199] Step 702: Obtain the optimal constant order numerical value and its corresponding mechanical parameters determined in advance through pre-inversion of the constant order fractional order model, and construct the global optimal constant order model as the first constant order benchmark.

[0200] Specifically, this step establishes the first reference frame for comparative verification, namely, the benchmark model representing the highest fitting level of the traditional technical path. During the initialization phase of the entire inversion process, the system forcibly disabled the time evolution mechanism of the order, locking the order decay rate to zero, degenerating into a classic constant-order fractional-order model. Using this degenerated model, global optimization was performed on 5 years of historical observation data to find a single constant-order value and its corresponding elastic modulus and other mechanical characteristics that best fit the measured displacement within that observation period. In the numerical example above, the optimal constant-order value obtained from the pre-inversion optimization was 0.35, with a corresponding elastic modulus of 24000 MPa. The system extracted this parameter combination, representing the best performance limit of the traditional constant-order model, solidified it, and defined it as the first constant-order benchmark. This benchmark is used to rigorously evaluate whether the variable-order model of this invention has achieved a substantial improvement in fitting accuracy compared to the globally optimal constant-order fitting under the same historical training data, avoiding unfairness in comparative verification due to inappropriate benchmark selection.

[0201] Step 703: Calculate the time-weighted average of the order evolution model corresponding to the final inversion parameter set over the entire observation period, and use the time-weighted average as a fixed order to construct an equivalent average constant-order model as the second constant-order benchmark.

[0202] In this embodiment, the system continues to construct a second reference frame for separating the verification variables. The final inversion parameter set is extracted, keeping the mechanical parameter set unchanged, and only the order evolution model it contains is subjected to time-domain averaging. Over the entire observation period from 0 to 5 years, the algorithm performs definite integral calculations on the determined parameterized or semi-parameterized order evolution function, and divides the total area of ​​the integral by the total duration of the observation period to obtain the mathematical expectation of the dynamically changing order sequence in the time dimension, i.e., the time-weighted average. Using the time-weighted average as a fixed order to replace the original variable-order evolution mechanism, combined with the unchanged mechanical parameters, the system assembles an equivalent average constant-order model, which is defined as the second constant-order benchmark. This benchmark is established to isolate the time-varying variable of the order under the stringent condition of consistent mechanical background parameters, thereby isolating the contribution of the action of simply changing the order from a constant to an evolution function to the improvement of the long-term prediction trend.

[0203] Step 704: Using the variable-order fractional-order forward model corresponding to the final inversion parameter set, the first constant-order benchmark, and the second constant-order benchmark, the displacement prediction comparison and verification of the long-term service conditions exceeding the observation period are carried out, and the long-term displacement prediction comparison results are obtained and output.

[0204] Specifically, this is the verification step for performing extrapolation simulations. The prediction timeline is significantly extended from year 5 to the 20th year of the dam's service life or even longer, and the reservoir's routine water level history for the next 15 years and long-term environmental temperature cycle fluctuations are used as external load conditions input into the computational environment. Three independent computation engines—the variable-order fractional-order forward model, the first constant-order benchmark, and the second constant-order benchmark—are activated to perform deformation calculations for future periods in parallel, generating three theoretical prediction curves of dam crest displacement that extend over time.

[0205] By comparing the development trends of the three curves, the technical effect of the order evolution mechanism introduced in this invention can be clearly observed. Since both the first and second constant-order benchmarks mathematically employ fixed fractional derivative kernel functions, their creep compliance maintains an unbounded growth trend resembling an approximately power function over an infinitely long time axis. This causes the predicted curves of traditional constant-order models to exhibit a continuous, monotonically increasing displacement value over time beyond the observation period, violating the fundamental physical and engineering understanding that the creep deformation of concrete materials inevitably tends towards convergence and stability after sufficient aging and hardening.

[0206] Conversely, because the order of the variable-order model of this invention gradually decreases and smoothly approaches the limiting order in the later stages of service, the activity of Maxwell viscous flow is dynamically suppressed. The predicted displacement curve of the variable-order model closely matches the measured data at the end of the historical observation period, while in the long-term prediction stage of extrapolation, the slope of the curve gradually flattens, the displacement growth rate slows down significantly over time, and eventually asymptotically converges to a physically reasonable upper limit of deformation. This comparative result strongly demonstrates that the variable-order model has irreplaceable technical advantages in overcoming the long-term extrapolation divergence defects and improving the reliability of structural life-cycle safety assessments.

[0207] This application introduces a variable-order fractional constitutive relation and its corresponding order evolution model. By using absolute time, equivalent hydration age, or normalized time as driving variables, the fractional order dynamically decays with the service life, mapping the real physical process of the gradual weakening of the material's viscous flow properties and improving the physical fidelity of the model's description of the long-term deformation state of the dam. It solves the problem that fixed-order models are difficult to characterize the aging and decay of concrete.

[0208] This application employs a piecewise constant-order time-domain discretization framework and a piecewise convolution superposition algorithm. It freezes the order within sub-intervals adaptively divided based on an equal-error criterion, ensuring that the material strictly satisfies the non-aging time-invariant condition locally. It legally introduces the Boltzmann superposition principle for partial strain integration, effectively resolving theoretical paradoxes at the underlying mathematical mechanism and enhancing the theoretical rigor of the forward modeling engine. It also resolves the theoretical conflict arising from the direct integration of continuously variable-order aged materials violating the superposition principle.

[0209] This application constructs an adaptive grouping regularization optimization engine based on dual dimensionless normalization. It performs dimensionality reduction by normalizing the parameter reference scale and sensitivity norm, and dynamically allocates penalty weights in the pure mathematical space based on the inverse of relative sensitivity. This eliminates the dimensional inconsistency defect caused by the mixing of stress dimensions and dimensionless variables, balances the optimization trajectories of high and low sensitivity parameters, and achieves efficient convergence of the inversion solution and high stability of extrapolation prediction. It solves the problems of dimensional curse and sensitivity imbalance caused by joint inversion of multiple heterogeneous parameters.

[0210] The preferred embodiments of the present invention have been described in detail above. The various specific technical features described in the above embodiments can be combined in any suitable manner without contradiction. To avoid unnecessary repetition, the present invention will not further describe the various possible combinations.

Claims

1. An adaptive inversion method for dam creep parameters based on a variable-order fractional-order constitutive model, characterized in that, include: Obtain the time series of measured displacements at multiple measuring points and the time history of load and environmental interaction, and construct a three-dimensional finite element model; Based on the material partitioning information of the three-dimensional finite element model, a variable-order fractional-order Burgers constitutive relation reflecting the creep activity decay characteristics of concrete is constructed, which includes an order evolution model. The service life of the dam corresponding to the three-dimensional finite element model is obtained, and it is discretized into multiple sub-intervals. Under the piecewise constant order framework, the superposition principle is applied to each sub-interval. Based on the variable order fractional order Burgers constitutive relation, the time history of load and environmental action, and the three-dimensional finite element model, the deviatoric strain increment of each sub-interval is calculated to obtain the variable order fractional order forward model. Construct an inversion objective function and divide the parameters to be inverted into a mechanical parameter set and a fractional-order evolution parameter set; Based on the measured displacement time series at multiple measurement points, the variable-order fractional-order forward model and the inversion objective function, the mechanical parameter set and the fractional-order evolution parameter set are alternately optimized to obtain and output the final inversion parameter set; Within a piecewise constant-order framework, the superposition principle is applied inter-interval. Based on the variable-order fractional-order Burgers constitutive relation, load and environmental action time histories, and a three-dimensional finite element model, the deviatoric strain increments of each sub-interval are calculated, including: Based on the variable-order fractional-order Burgers constitutive relation, the constant-order creep compliance kernel function of each subinterval is determined; Based on the time history of load and environmental action and the three-dimensional finite element model, the deviatoric stress increment in each sub-interval is determined; For any kth subinterval, the cumulative creep deviatoric strain of the preceding subinterval is obtained by performing a constant-order convolution integral and summing the constant-order creep compliance kernel function of each preceding subinterval and the deviatoric stress increment in each preceding subinterval. The creep partial strain increment of the current interval is obtained by performing a constant-order convolution integral using the constant-order creep compliance kernel function of the k-th sub-interval and the partial stress increment in the k-th sub-interval. Based on the deviatoric stress increment in the k-th sub-interval and the pre-configured elastic shear modulus, the instantaneous elastic deviatoric strain is calculated. The cumulative creep deviatoric strain of the previous sequence, the creep deviatoric strain increment of the current interval, and the instantaneous elastic deviatoric strain of the current interval are superimposed to obtain the deviatoric strain increment of each sub-interval and the total deviatoric strain at the end of the k-th sub-interval.

2. The method according to claim 1, characterized in that, The order evolution model is used to characterize the evolution of the fractional order of Maxwell branches in the variable-order fractional Burgers constitutive relation as a function of the driving variables. The order evolution model satisfies the boundary condition: when the driving variable is zero, the fractional order is equal to the initial order; As the driving variable approaches infinity, the fractional order asymptotically approaches the limit order. The initial order is greater than the limiting order.

3. The method according to claim 2, characterized in that, The driving variable is the absolute service time; The order evolution model is defined as an exponential decay function: α(t)=α ∞ +(α0-α ∞ )e (-βt) ; Where α(t) is the fractional order corresponding to the absolute service time t, α0 is the initial order, and α ∞ β is the limiting order, and β is a positive parameter characterizing the aging rate.

4. The method according to claim 2, characterized in that, The driving variable is the equivalent hydration age; Obtain the measured temperature history during the time history of load and environmental effects; The equivalent hydration age is calculated using the Arrhenius temperature rate function based on the measured temperature history. The calculation formula is as follows: t eq (t)=∫0 t exp[(E a / R)*(1 / T ref -1 / T(τ))]dτ; Among them, t eq (t) represents the equivalent hydration age corresponding to the absolute service time t, T(τ) represents the measured temperature history at time τ, and E a Let T be the apparent activation energy, R be the universal gas constant, and T be the apparent activation energy. ref This is the preset reference temperature; The order evolution model is defined as: α(t)=α ∞ +(α0-α ∞ )exp(-β*t eq (t)); Where α(t) is the fractional order, α0 is the initial order, and α ∞ β is the limiting order, and β is the decay rate constant of creep activity at the reference temperature.

5. The method according to claim 2, characterized in that, The driving variable is normalized time; Normalized time is calculated using absolute service time and the total length of pre-configured observation periods; The order evolution model is a semi-parametric expression based on Bernstein polynomial basis functions, and the formula is: ； Where α(t) is the fractional order, α0 is the initial order, and α ∞ Let K be the limit order, K be the preset expansion order, ξ be the normalization time, and B be the limit order. k,K (ξ) are K-order Bernstein polynomial basis functions, c k This is the k-th Bernstein control coefficient; To ensure the physical plausibility of the monotonically decreasing order of fractional orders, the Bernstein control coefficients satisfy the following constraints: c0=1; c0≥c1≥c2≥...≥c k ≥0.

6. The method according to claim 1, characterized in that, The dam's service life is discretized into multiple sub-intervals, including determining the endpoints of non-uniform sub-intervals using an equal error distribution criterion based on the order of change rate: Calculate the absolute value of the rate of change of the order evolution model over the dam's service life; By integrating the absolute value of the rate of change over the absolute service time, the order arc length function and the total order arc length corresponding to the dam's service cycle are calculated. The total order arc length is divided equally according to the preset total number of sub-intervals, resulting in multiple arc length division points; By using the inverse function of the order arc length function, multiple arc length division points are mapped back to physical time coordinates to obtain the endpoints of non-uniform sub-intervals for each sub-interval.

7. The method according to claim 1, characterized in that, Before constructing the inversion objective function, the process also includes dimensionless normalization of each parameter to be inverted in the mechanical parameter set and the fractional evolution parameter set: Obtain the pre-configured reference scale for each parameter to be inverted; Divide each parameter to be inverted by its corresponding reference scale to obtain the corresponding dimensionless parameter. The partial derivatives of the displacement of each measuring point with respect to each parameter to be inverted are calculated using a variable-order fractional forward model. The partial derivatives are then multiplied by their respective reference scales to obtain the normalized sensitivity vectors for each parameter to be inverted.

8. The method according to claim 1, characterized in that, Output the final inversion parameter set and perform creep characteristic prediction, including: Obtain the optimal constant-order numerical values ​​and their corresponding mechanical parameters determined in advance through pre-inversion of the constant-order fractional-order model, and construct the globally optimal constant-order model as the first constant-order benchmark; Calculate the time-weighted average of the order evolution model corresponding to the final inversion parameter set over the entire observation period, and use this as a fixed order to construct an equivalent average constant-order model as the second constant-order benchmark; The variable-order fractional-order forward model corresponding to the final inversion parameter set, the first constant-order benchmark, and the second constant-order benchmark were used to conduct displacement prediction comparison and verification for long-term service conditions that exceed the observation period, and the long-term displacement prediction comparison results were obtained and output.

9. The method according to claim 5, characterized in that, The parameters to be inverted also include the Bernstein control coefficient set, which is independent of the mechanical parameter set and the fractional evolution parameter set; The construction of the inversion objective function also includes the construction of a second-order difference smoothness penalty term for the Bernstein control coefficient set; A second-order difference smoothness penalty term is added to the inversion objective function to suppress oscillatory overfitting of the Bernstein control coefficients.

Images