Method for drawing jetty structure vulnerability curve based on random load combination

Abstract

This invention relates to the field of disaster prevention and mitigation technology for marine engineering structures, and discloses a method for plotting vulnerability curves of guide seawall structures based on random load combinations. The method includes: constructing a nonlinear finite element model of a deep-water guide seawall; introducing a two-dimensional phase coupling vector to generate random dynamic load samples containing temporal coupling characteristics of earthquakes and waves; inputting the samples into the model, tracking the structural frequency drift in real time and considering pore flow abrupt change damping, and calculating the multidimensional dynamic response of the structure; constructing a multidimensional dynamic limit state surface that shrinks with cumulative plastic damage, and determining the failure state based on the relative relationship between the dynamic response and the state surface; estimating the failure probability using extreme value statistics theory, and establishing a probabilistic demand model containing cross-coupling terms of earthquake and wave intensity; and plotting a vulnerability surface describing the change of the failure probability with bivariate strength indices based on this model. This invention can accurately assess the dynamic reliability of deep-water guide seawalls under the coupling effects of multiple disasters.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of disaster prevention and mitigation of marine engineering structures, in particular to a method for drawing a breakwater structure vulnerability curve based on random load combination. BACKGROUND

[0002] As an important marine infrastructure, deepwater breakwaters have long been serving in complex dynamic environments and are threatened by both earthquakes and waves. Vulnerability analysis is an effective means to assess the failure probability of such structures under random dynamic loads. Existing methods for vulnerability assessment of breakwater structures are usually based on finite element numerical simulation, using Monte Carlo method or incremental dynamic analysis to establish the probability relationship between environmental load intensity and structural damage index.

[0003] When dealing with the problem of combined action of earthquake and wave, the existing technology usually assumes that the two are statistically independent or only performs simple linear superposition, ignoring the phase coupling mechanism in the time domain. This simplified processing method is difficult to cover the extreme unfavorable working condition when the peak of ground motion and the trough of wave overlap, that is, it cannot accurately simulate the critical state of minimum buoyant force and maximum horizontal thrust of the structure, thus leading to the evaluation result of breakwater anti-sliding stability being unsafe. In terms of numerical model, the existing method uses equivalent linearization or constant parameter model to deal with the fluid-structure coupling problem, which fails to reflect the nonlinear damping changes caused by the transition of pore flow state of rockfill from laminar to turbulent under strong earthquake action, and also fails to fully consider the dynamic amplification effect caused by the resonance locking of self-vibration frequency drift due to stiffness degradation and wave eigenfrequency caused by cumulative damage of the structure.

[0004] In addition, the traditional failure criterion is based on fixed displacement or stress threshold, ignoring the material property that the ultimate resistance of the structure will gradually decay with the accumulation of plastic damage, resulting in the inability to truly evaluate the residual carrying capacity of the structure in the later stage of dynamic time history. Existing probability regression models also often lack quantitative description of multi-disaster source cross-coupling terms, limiting the applicability and accuracy of multi-dimensional vulnerability surface in complex marine environments. SUMMARY

[0005] The application provides a jetty structure vulnerability curve drawing method based on random load combination, which comprises constructing a deep water jetty nonlinear finite element analysis model and generating a random dynamic load sample sequence containing time domain coupling characteristics. In the generation of the load sample sequence, a reference earthquake and a wave power spectral density function are selected, a two-dimensional phase coupling vector covering the most unfavorable phase difference working condition of waves and earthquakes is defined, and a phase coupling sample pair is generated by using a Latin hypercube sampling method; the phase coupling sample pair is linearly combined with an independent micro random phase sequence, and the earthquake acceleration and wave force time history after time window alignment are synthesized by using a trigonometric series superposition method, so that the random dynamic load reflecting the time domain overlap effect of the earthquake peak value and the wave trough is constructed.

[0006] The random dynamic load sample sequence is batch input into the nonlinear finite element analysis model to calculate the multi-dimensional dynamic response of the structure. In the dynamic analysis process, the relative seepage velocity of the pore fluid is calculated in real time to identify the pore flow state transition, and the damping matrix is adjusted accordingly to reflect the nonlinear energy dissipation characteristics of the fluid inside the rockfill under strong earthquake action; at the same time, the natural frequency drift of the structure caused by cumulative damage is monitored, the overlap of the frequency drift trajectory and the high energy band of the wave spectrum is identified, and the effective load vector is dynamically modified according to the wave spectrum energy capture function to reflect the nonlinear influence of the fluid-structure coupling resonance effect on the structure response.

[0007] In the failure judgment stage, the method sets a structure limit state threshold and constructs a multi-dimensional dynamic limit state surface. The normalized settlement of the embankment crest, the maximum slip amount of the armor layer, and the average pore pressure ratio of the key area of the foundation are selected as independent state components to construct an outer convex closed surface composed of points with a limit state function value of zero. In order to reflect the degradation of the structure resistance with damage, a decay function that monotonically decreases with the increase of the cumulative plastic damage variable is defined, and the initial static threshold is reduced in real time, so that the multi-dimensional dynamic limit state surface shrinks to the origin with the accumulation of damage. Based on the relative relationship between the multi-dimensional dynamic response of the structure and the dynamic limit state surface, the algebraic safety margin is calculated and the structure failure state is judged.

[0008] For samples that do not directly fail, the method uses the extreme value statistical theory to perform parameter fitting on the sequence composed of maximum risk indicators by using a generalized extreme value distribution function, and calculates the exceedance probability of the maximum risk indicator greater than or equal to zero as the estimated failure probability. Subsequently, regression analysis is performed on the multi-dimensional dynamic response of the structure to establish a probability demand model containing a cross-coupling term, which is the product of the logarithmic values of the earthquake intensity index and the wave intensity index. Finally, according to the probability demand model and the structure limit state threshold, a vulnerability surface describing the change of the failure probability with the earthquake intensity index and the wave intensity index is drawn, and the annual average failure frequency of the jetty structure is calculated by full probability integration combined with the seismic hazard curve of the target sea area and the wave long-term probability density function.

[0009] This invention provides a method for plotting the vulnerability curve of a guide embankment structure based on random load combinations. It has the following beneficial effects: 1. This invention generates random dynamic load samples by introducing a two-dimensional phase coupling vector, achieving specific phase combination control of seismic acceleration time history and wave force time history in the time domain. This method can simulate the most unfavorable condition where the peak ground motion and wave trough coincide, i.e., the moment when the maximum hydrostatic pressure difference on the seaward side of the guide dike and the minimum effective self-weight occur simultaneously. This overcomes the deficiency of traditional methods where independent random phase superposition cannot cover extreme load combinations, improving the conservatism and safety of the assessment of the guide dike's anti-slip stability.

[0010] 2. This invention constructs a flow regime change damping model and frequency locking mechanism in finite element analysis, improving the physical realism of the calculation of the nonlinear dynamic response of the structure. By updating the pore fluid damping in real time to reflect the flow regime transition from laminar to turbulent, and dynamically correcting the wave load based on the natural frequency drift caused by the cumulative damage to the structure, this method can capture the resonance amplification effect between structural stiffness degradation and the spectral characteristics of external loads, thereby more accurately predicting the failure evolution process of deep-water guide dikes under the coupled action of strong earthquakes and large waves.

[0011] 3. This invention establishes a multidimensional dynamic limit state surface that shrinks with cumulative damage and employs a probabilistic demand model including cross-coupling terms to improve the accuracy of vulnerability assessment. The limit state threshold is defined as a function that monotonically decreases with the accumulation of plastic damage, consistent with the resistance decay characteristics of geotechnical materials under cyclic loading. Simultaneously, a coupling term between earthquake and wave intensity is introduced into the regression model, quantifying the nonlinear synergistic influence of the two disaster sources on the structural response, enabling the final vulnerability surface to more objectively reflect the probability of structural failure under complex environments. Attached Figure Description

[0012] Figure 1 A flowchart illustrating a method for plotting the vulnerability curve of a guide embankment structure based on random load combinations, provided in an embodiment of the present invention. Figure 2 This is a flowchart illustrating the generation of a three-dimensional vulnerable hypersurface of a guide embankment structure incorporating coupling effects, as described in an embodiment of the present invention. Detailed Implementation

[0013] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0014] Please see the appendix Figure 1This invention provides a method for plotting the vulnerability curve of a guide embankment structure based on random load combinations, comprising the following steps: S10, construct a multidimensional non-stationary random load model with phase coupling mechanism, and generate a time-domain correlated wave and earthquake joint time sample set by introducing phase coupling vector; S20. Establish a nonlinear dynamic equilibrium equation based on fluid regime change, discretize the guide embankment structure into a liquid-containing porous medium system, and define the fluid regime-dependent nonlinear damping matrix and the damage-driven time-varying stiffness matrix. S30 calculates the spectral migration trajectory and fluid-structure interaction resonance locking effect, tracks the instantaneous natural frequency of the structure in real time, and dynamically corrects the effective load vector based on the wave spectrum energy capture function. S40, construct the dynamic limit state surface and perform the first-pass failure determination, define the limit state equation that shrinks with cumulative damage in the generalized state space, and identify the intersection of the random load trajectory and the limit state surface. S50 generates a multidimensional vulnerable hypersurface. Based on parametric scanning and probabilistic statistical analysis, a vulnerability function of the guide embankment structure containing frequency migration coupling terms is obtained by fitting.

[0015] In this embodiment, the above method constitutes a closed-loop dynamic reliability assessment system. The system first converts the wave field and seismic field in the physical environment into stochastic process inputs required for numerical calculation through a probabilistic model. In step S10, unlike the traditional linear superposition of loads, this method simulates the specific encounter mode of the seismic peak and wave trough on the time axis through a phase coupling mechanism, thereby covering the most unfavorable working conditions caused by phase interference.

[0016] Subsequently, the data stream enters the dynamic response calculation module. In step S20, the guide embankment structure is characterized as a nonlinear system in which physical parameters evolve over time. This step identifies the moment when the pore flow state changes from laminar to turbulent by calculating the Reynolds number of the pore fluid in real time, and abruptly changes the coefficients in the damping matrix accordingly to reflect the nonlinear energy dissipation characteristics of the fluid inside the rockfill under strong earthquake action.

[0017] Step S30, as the core processing step, is performed alternately with step S20. As the structural stiffness degrades due to accumulated damage, the system's natural frequency drifts. This step monitors whether the frequency drift trajectory falls into the high-energy frequency band of the wave spectrum at the current moment. Once spectral overlap occurs, the system determines that the fluid-structure interaction resonance locking mechanism is triggered and nonlinearly amplifies the wave force acting on the structure. This amplification effect, in turn, further exacerbates the damage and stiffness degradation of the structure, forming a positive feedback coupling mechanism.

[0018] Step S40 abandons the single displacement threshold judgment criterion. In this step, the relative positional relationship between the load trajectory and the dynamic limit state surface is calculated in real time in the state space composed of generalized load effect and generalized resistance. If the load trajectory crosses this contraction surface at any time, the random simulation is judged as a failure event. This judgment method takes into account the temporal matching of resistance attenuation and random load fluctuation.

[0019] Finally, step S50 performs statistical regression on the calculation results of a large number of random samples. By performing gridded scanning within the parameter space composed of wave intensity index and seismic intensity index, the system outputs a three-dimensional vulnerable hypersurface containing a coupled resonance correction term, thereby realizing the probabilistic safety assessment of the guide embankment structure in a complex random environment.

[0020] In constructing a multidimensional non-stationary stochastic load model with a phase coupling mechanism (step S10), the determination of the evolving power spectral density function (sub-step S110) is specifically performed as follows: For wave loads, the JONSWAP spectral model, which characterizes the properties of incompletely grown wind waves at finite water depths, is selected as the frequency domain reference. This model uses peak-shape parameters to correct the peaks in the standard Pierson-Moskowitz spectrum, and its one-sided power spectral density function... Defined as: ; in, Indicates the angular frequency of the wave. Represents gravitational acceleration, spectral peak frequency Based on the design effective wave height within the recurrence period of the target sea area Spectral peak period The statistical relationship is determined; Phillips constant With peak shape parameters Values ​​were fitted based on empirical data on wind speed and wind zone length (Fetch) in the sea area, where... The value of is typically between 1.0 and 7.0, used to adjust the concentration of spectral energy in the frequency domain; peak width coefficient Take 0.07 to the left of the harmonic peak and 0.09 to the right.

[0021] For bedrock ground motion, the Clough Penzien spectral model, which includes low-frequency filtering correction, is selected. This model is based on the Kanai Tajimi spectrum with a high-pass filter in series to eliminate low-frequency integral drift distortion. Its power spectral density function... Expressed as: ; Among them, white noise spectral intensity It can be obtained directly from the peak ground acceleration (PGA) of the design earthquake through integration: ; excellent frequency of the site Damping ratio The filtering parameters are determined based on the shear wave velocity profile obtained from the geological survey of the engineering site. and Based on the baseline correction requirements of the seismic records, the settings ensure that the generated displacement time histories do not diverge.

[0022] A non-stationary power spectrum is constructed. Given the time-varying characteristics of natural disaster loads, an intensity modulation function is employed. Time-domain amplitude modulation is applied to the aforementioned stationary spectrum. For seismic motion, a three-segment envelope function is defined. The duration of its rising phase Duration of the steady phase and attenuation coefficient It is determined based on the attenuation relationship between the magnitude of the potential earthquake source and the epicentral distance.

[0023] Finally, the evolutionary power spectral density functions of waves and earthquakes are synthesized through the following product operations: and ; in The wave intensity envelope function represents the significant wave height in the above wave spectrum, especially when considering long-duration processes like storm surges. With spectral peak period Further set as time Slowly changing function and This leads to the evolution of the power spectral density function. It not only exhibits non-stationarity in intensity, but also includes the migration characteristics of spectral components over time.

[0024] In constructing a multidimensional non-stationary stochastic load model with a phase coupling mechanism (step S10), the phase coupling vector is introduced and random sampling is performed (sub-step S120), which is specifically executed as follows: Define independent stochastic phase and coupled phase variables. For two independent stochastic processes, wave and earthquake, generate two sets of independent phase variables within the interval... Uniformly distributed random phase angle sequence and ,in This corresponds to discrete frequency points. Based on this, a two-dimensional phase coupling vector is introduced. This vector does not change with frequency. The changes, as an overall offset, are added to their respective random phase sequences to shift the envelope peak position of the random process on the time axis.

[0025] Specifically, phase angle With time delay The relationship between them is satisfied ,in This corresponds to the dominant frequency (i.e., the peak frequencies) of the random process. By adjusting... and The difference can be used to simulate the time interval between the arrival times of different types of load peaks. A sampling space for the most unfavorable phase condition is constructed. To cover the various time series combinations that the guide embankment structure will encounter during its service life, phase coupling variables are set. and The range of values ​​covers a complete fluctuation cycle, that is... In numerical simulations, the focus is on phase difference. The working condition corresponds physically to the overlap of the peak ground acceleration moment of the earthquake and the trough moment of the wave (i.e. the moment of minimum water level).

[0026] At this point, the guide seawall experiences maximum buoyancy force, minimum effective self-weight, and maximum hydrostatic pressure difference towards the sea, constituting the most unfavorable boundary condition for anti-slip stability. Stratified random sampling is performed. The Latin hypercube sampling (LHS) method is used to sample the coupling vectors. Sampling is performed to generate Group phase sample pairs At the same time, corresponding generation A set of independent micro-random phase sequences By linearly combining the coupling phase and the microscopic phase, we obtain the first... The final phase parameters required for this simulation: ; Among them, the combined phase This will be directly substituted into the cosine function term of the subsequent spectral representation. Wherein, This ensures that the generated samples satisfy the statistical properties of the target power spectral density (i.e., second-order moment consistency), while This ensures that the generated sample set covers all temporal interference scenarios from "earthquake at the wave crest" to "earthquake at the wave trough". The specific algorithm for Latin hypercube sampling is well-known in this field and will not be elaborated here.

[0027] In constructing a multidimensional non-stationary stochastic load model with a phase coupling mechanism (step S10), the synthesis of joint load samples (sub-step S130) is specifically performed as follows: Determine the time-frequency discrete parameters. Set the cutoff frequency based on the energy distribution range of the evolving power spectral density. and divide it into Each frequency band, bandwidth The center frequency point is taken as At the same time, the time step is set according to the accuracy requirements of the dynamic analysis. Constructing discrete time series .

[0028] Calculate the time-varying amplitude modulation coefficients for each discrete frequency point. and time Using the evolution power spectral density function Calculate the corresponding physical amplitude: ; in, This formula represents discrete frequency points and is applicable to both wave and seismic spectra. Substituting the values ​​into the appropriate values... and The wave amplitude sequence can then be obtained. With earthquake acceleration amplitude sequence .

[0029] The stochastic process is synthesized and transformed into nodal load vectors. The combined phase parameters generated in step S120 are then used. Time history samples are generated by superimposing cosine series and mapped to spatial force vectors acting on the nodes of the structural finite element.

[0030] For wave loads, the time history of wave height is first synthesized. ; ; in, The wave amplitude sequence is represented, and then, based on this wavefront height time history, the elevation of the upstream face along the guide breakwater is calculated using the Godad wave pressure distribution formula or linear wave theory. Distributed hydrodynamic pressure Then, it is integraled over the area and transformed into an equivalent nodal force vector acting on each node of the guide dike on the water-facing side. For seismic action, first synthesize the time history of horizontal acceleration in the bedrock. ; ; in, The seismic acceleration amplitude sequence is represented, and then, according to d'Alembert's principle, the mass matrix of the structure is used. (Including additional water mass), calculate the equivalent seismic inertial force vector. ,in This is an indicator vector for the direction of seismic input.

[0031] Construct the final load input set. Superimpose the wave node force vectors and the seismic inertial force vectors to form the... Right-hand side terms of the nonlinear dynamic equation for the group of samples Repeat the process. Next, a load sample library containing different phase interference modes is established. The use of Fast Fourier Transform (FFT) to accelerate the trigonometric series summation process is a standard numerical algorithm in this field and will not be elaborated upon here.

[0032] In establishing the nonlinear dynamic equilibrium equations based on the abrupt change in flow regime (step S20), the discretization modeling of the porous medium of the guide embankment (sub-step S210) is specifically performed as follows: Discretization computational domain and element type selection: The cross-section of the guide embankment project is divided into three physical regions: rockfill, foundation soil, and external water. For the rockfill and foundation regions, a dynamic theory based on Biot porous media is adopted. Discretized using two-dimensional plane strain elements (solid displacement and pore water pressure). Each node contains two degrees of displacement freedom. and a pore water pressure degree of freedom This allows for the complete coupled simulation of skeleton deformation and pore water seepage.

[0033] Construct the solid-state skeleton mass matrix. For each discretized element, calculate the inertial contribution of its solid phase. Element skeleton mass matrix. Obtained through numerical integration: ; in, Porosity The density of the rubble particles, It is a displacement interpolation shape function.

[0034] The coupling term between the fluid's added mass and the hydrodynamic pressure is calculated. The inertial effect of the fluid on the structure is considered, and the process is divided into two parts: the fluid inside the pores and the external free water body.

[0035] For the interior of the rockfill, a pore-based tortuous courtyard is introduced. Additional mass density This is then superimposed onto the total mass matrix to characterize the inertial drag generated when the pore fluid accelerates relative to the skeleton.

[0036] For external water bodies, instead of directly dividing the fluid grid, the additional mass method is used to apply the hydrodynamic pressure equivalently to the nodes on the water-facing surface of the guide dike. Specifically, the additional mass per unit area distributed along the normal direction of the slope is calculated using the modified Westergaard parabolic distribution formula. : ; in Because of the water depth, It is the vertical distance from the water surface. Let be the slope angle of the guide embankment. The mass distribution is integrated and concentrated onto the diagonal of the mass matrix of the nodes near the water surface.

[0037] The total mass matrix of the system. This is formed by assembling the solid mass, the equivalent inertia term of the pore fluid, and the added mass of the external water body. This matrix serves as the dynamic equation. The inertia coefficients directly participate in the subsequent time-step integration. The specific size control of the mesh generation and the setting of artificial boundary conditions are standard preprocessing techniques in this field and will not be elaborated upon here.

[0038] In establishing the nonlinear dynamic equilibrium equation based on flow regime change (step S20), the construction of the flow regime-dependent damping matrix (sub-step S220) is specifically performed as follows: Construct the basic Rayleigh damping matrix. Establish a linear damping matrix to address the fluid viscosity effect under micro-deformation of the structural skeleton. : ; in, Let be the tangent stiffness matrix at the current time step. and Determined based on the fundamental frequency of the guide embankment and the target damping ratio. This part mainly reflects the energy dissipation due to internal friction of the medium.

[0039] Calculate the relative seepage velocity at the unit level. At each time step of the dynamic analysis, based on... The nodal solution of the format is used to calculate the seepage velocity vector of the pore fluid relative to the solid skeleton. Using the generalized form of Darcy's law, the relative velocity is determined by the current pore water pressure gradient. and skeletal acceleration Sure: ; in The permeability coefficient of the medium. It is the density of water.

[0040] A nonlinear turbulent drag term is constructed. To reflect the sharp increase in nonlinear drag caused by the transition from laminar to turbulent flow in pores under strong earthquakes and high flow velocities, a modified Ergun equation is used to describe the interaction forces between the fluid and the framework. This force includes linear laminar terms and nonlinear turbulent terms: ; Among them, laminar flow resistance coefficient turbulent drag coefficient .here For fluid dynamic viscosity, The equivalent diameter of the particle. This refers to the relative flow velocity.

[0041] Assemble the instantaneous equivalent damping matrix. To solve the finite element equations, the nonlinear resistance terms are linearized, and the equivalent flow is extracted.

[0042] Body damping coefficient Define the equivalent permeability coefficient. Decrease with flow rate: ; in, The initial Darcy permeability coefficient, This is the turbulence coefficient.

[0043] The permeability coefficient is converted into an equivalent damping coefficient according to the generalized Darcy's law, and the fluid-structure interaction damping matrix is ​​updated accordingly. The corresponding terms in the solution. During the solution process, as... As the value increases, the equivalent damping increases significantly, automatically realizing the physical evolution from linear Darcy flow (low Reynolds number) to nonlinear Forchheimer flow (high Reynolds number) without the need for manually setting a sudden change threshold. The final system damping matrix is... And it is updated in real time in each Newton-Raphson iteration step.

[0044] In establishing the nonlinear dynamic equilibrium equation based on abrupt changes in flow regime (step S20), the damage-driven time-varying stiffness mechanism is introduced (sub-step S230), which is specifically executed as follows: A damage evolution model based on energy dissipation is constructed. To describe the stiffness degradation of the rockfill under cyclic dynamic loading, a star-shaped damage variable is defined. Its evolution law is based on cumulative plastic shear strain. In each time increment step, the damage increment is calculated. : ; in, The equivalent plastic shear strain increment for the current step is determined by the plastic flow direction and plastic multiplier in the constitutive integration algorithm; For reference strain parameters; This represents the maximum shear strain amplitude experienced up to the present moment. For the damage threshold shear strain; and The material constant used to control the damage rate. When hour, This formula indicates that irreversible damage to material stiffness only occurs when significant plastic deformation occurs and the deformation amplitude exceeds historical extremes. The elastic predictive modulus of damage coupling is defined. In the prediction stage of the elastoplastic constitutive calculation, the damage-corrected instantaneous shear modulus is introduced. and bulk modulus .

[0045] For shear modulus, considering confining pressure dependence and damage effect: ; in, For reference shear modulus, For the cumulative damage variable, considering the slow degradation of the volumetric measure of the rockfill, a weaker damage correlation or only considering the confining pressure correlation is used for the bulk modulus: ; in Volumetric damage coefficient (value) The current mean effective principal stress, Atmospheric pressure.

[0046] The uniform tangent stiffness matrix is ​​derived. An implicit integration algorithm is used to solve the stress update process. This is done by constructing the element tangent stiffness matrix. At that time, it not only includes the traditional elastoplastic matrix Furthermore, it is necessary to introduce variables due to damage. Additional softening term resulting from strain variation. Specific uniform tangent modulus. Expressed as: ; in, It is an elastic-plastic matrix. For the current stress tensor, This is the derivative of the damage with respect to the strain in the pore, reflecting the softening term.

[0047] This correction term guarantees the second-order convergence rate of the Newton-Raphson iteration, avoiding numerical non-convergence caused by lag in stiffness updates or iteration calculations. Update the global dynamic characteristics. Update the stiffness matrix of each element. The instantaneous tangential stiffness matrix assembled into the whole system With damage variables The accumulation, The eigenvalues ​​gradually decrease. In solving the dynamic equations, this means that the instantaneous natural frequency of the structure... It will decrease monotonically over time. This time-varying frequency characteristic will serve as a key state parameter for subsequent resonance locking, used to determine whether the structure has entered the vulnerable low-frequency resonance region.

[0048] After constructing the time-varying stiffness and damping model, the real-time tracking of instantaneous natural frequencies (sub-step S310) is executed as follows: Extract and assemble the instantaneous system matrix. During the nonlinear dynamic time history analysis, the sampling interval for frequency tracking is set. (For example, every 10 or 50 integration steps). At each sampling point, the current system tangent stiffness matrix is ​​extracted from the solver. To avoid frequency calculation errors caused by singular or negative eigenvalues ​​in the stiffness matrix during numerical iteration, a modified secant stiffness matrix or a small amount of numerical damping is introduced into the tangent stiffness when plastic yielding is detected, ensuring the positive definiteness of the matrix. Simultaneously, the constant system mass matrix assembled in step S210 is invoked. Fast frequency estimation based on the Rayleigh quotient method. Given the high computational cost of complete eigenvalue decomposition, an iterative Rayleigh quotient method is used to quickly update the fundamental frequency. The first-order mode shape vector calculated at the previous sampling time (or the initial one) is used. As a trial vector, substitute it into the current stiffness and mass matrices to calculate the current Rayleigh quotient. : ; in, The first-order mode shape vector, For a moment The tangent stiffness matrix, Given the system mass matrix, the first-order instantaneous natural circular frequency of the structure is approximately: For higher precision, inverse lteration can be used to adjust the mode shape vectors. By performing one or two iterations, high-precision fundamental frequency evolution data can be obtained with extremely low computational cost.

[0049] Time-domain smoothing filtering is performed. Because the tangent modulus of the nonlinear constitutive model exhibits numerical jumps during the loading and unloading transitions, the directly calculated original frequency sequence... It contains high-frequency numerical noise. Construct a time window with a width of... The moving average filter smooths the original frequency sequence in real time. ; in The number of data points within the window. The value is typically taken as 0.5 to 1.0 times the initial period of the structure. This step filters out non-physical numerical oscillations and preserves the true trend of macroscopic degradation of structural stiffness as damage accumulates.

[0050] Output structure time-varying modal characteristics. The smoothed instantaneous fundamental frequency. It is stored as a time series vector. This series quantitatively describes the overall stiffness degradation process of the guide embankment structure under the coupled action of strong earthquakes and large waves, caused by the attenuation of the shear modulus of the rockfill and the softening of pore water pressure. It serves as the core state variable for subsequent judgment on whether "frequency encounter" and resonance locking have occurred.

[0051] After obtaining the time-varying frequency trajectory of the structure, the definition of the wave spectrum energy capture function (sub-step S320) is executed as follows: Estimate the instantaneous equivalent modal damping ratio. This is based on the structural damping matrix. It includes nonlinear turbulence terms that vary with flow velocity, and needs to be converted into the equivalent damping ratio in modal space. To determine the resonance bandwidth, the first-order mode shape at the current moment is used. The calculation is performed using the energy equivalence relationship between generalized modal mass and modal damping: ; in, For a moment The equivalent modal damping ratio, It is the instantaneous natural circular frequency. The system damping matrix, which includes nonlinear terms, projects the unevenly distributed fluid resistance in physical space onto the main vibration modes of the structure, resulting in a scalar index characterizing the intensity of energy dissipation.

[0052] Define the effective resonant frequency band. Based on the concept of the half-power point in linear vibration theory, define the frequency range in which the structure is highly sensitive to external excitation energy. Considering the frequency "softening" characteristic of nonlinear systems, the current instantaneous fundamental frequency is taken. The width on both sides is one times the damping ratio, which is used as the integration interval: ; ; Calculate the wave energy capture density function. To quantify the potential destructive power of the external wave field on the structure at each moment, calculate the wave spectrum energy falling within the aforementioned resonant frequency band. Evolve the power spectral density function. exist By performing definite integrals over the interval, the energy capture function is obtained. .

[0053] ; in, and These correspond to the upper and lower limits of the resonant frequency band, respectively. The larger the value of the energy capture function, the higher the overlap between the main energy frequency component of the wave and the instantaneous natural frequency of the structure, and the stronger the tendency of the structure to resonate and amplify.

[0054] A normalized resonance risk factor is constructed. To facilitate cross-sectional comparisons between samples from different sea states and to eliminate the influence of total wave energy intensity, a dimensionless resonance risk factor is defined. : ; Traverse the entire dynamic time history and extract The time history curve reflects the dynamic coupling between the structural damage process and the evolution of the random wave spectrum. When this factor exceeds a preset critical value (e.g., 0.4), the time period is marked as a strong resonance coupling zone, and the corresponding wave and earthquake phase combination is identified as the most unfavorable combination that can induce the maximum dynamic response of the structure.

[0055] After identifying the critical moment with high resonance risk, the amplification correction of the dynamic load (sub-step S330) is specifically executed as follows: Construct a frequency domain local amplification function. Based on the resonance risk factors identified in step S320. and the corresponding central resonant frequency Construct an amplitude amplification factor function for a specific frequency band. The amplification factor is described in a Gaussian distribution in the frequency domain to ensure that targeted strengthening is applied only to the weakest frequencies of the structure, while keeping the load characteristics far from the resonance region unchanged. ; in, To design a safety margin factor (with a value greater than 10). The equivalent damping ratio at resonance is used to control the width of the amplification bandwidth.

[0056] Correct the power spectral density of the design load. Obtain the power spectral density functions of the original design wave field and the input ground motion, respectively. The original spectral density is modulated using the aforementioned amplitude amplification function to generate the most unfavorable envelope power spectrum. : ; Note the following: The power spectral density is squared because it is proportional to the square of the signal amplitude. This step artificially constructs a virtual environmental load spectrum with energy concentrated at the natural frequency of the structure after the current damage.

[0057] Time history reconstruction based on Inverse Fast Fourier Transform (IFFT). To generate a load sequence usable for time-domain finite element analysis, the frequency domain energy needs to be converted back to the time domain. The key lies in the selection of the phase angle. To reproduce the resonance locking effect, the phase angle spectrum at the moment of maximum resonance risk in the original load time history is directly extracted. The Spectral Representation Method is used to combine the corrected amplitude spectrum with the original poor phase spectrum. ; ; Alternatively, it can be generated directly through inverse FFT transformation. The resulting corrected load time history It possesses both targeted resonant energy enhancement and retains the phase characteristics in the original random wave field that lead to the most unfavorable response of the structure.

[0058] Design verification was performed. The generated modified wave force time history and modified seismic acceleration time history were used as new boundary conditions and input into the nonlinear finite element model of the guide embankment for secondary dynamic analysis. It was then checked whether the permanent displacement of the embankment crest, the slippage of the revetment blocks, and the extent of foundation liquefaction under the "resonance enhancement" load still met the engineering allowable standards. This step represents a shift from "probabilistic design" to "design based on the most unfavorable resonance mechanism."

[0059] To comprehensively assess the multiple failure risks of the guide embankment under the coupled effects of strong earthquakes and large waves, the construction of the dynamic limit state surface (sub-step S410) is specifically executed as follows: Define a multidimensional performance state vector. Select independent physical quantities characterizing the key failure modes of deep-water guide dikes to construct the system state vector. For a typical rockfill guide embankment structure, the following three components are selected: 1. Normalized levee crest settlement Defined as the cumulative vertical displacement at the center point of the embankment crest. With dam height The ratio is used to assess overtopping risk.

[0060] 2. Maximum slippage of the protective layer The peak value of the tangential relative displacement of the facing block unit relative to the underlying cushion layer is extracted to assess the risk of structural integrity loss.

[0061] 3. Average pore pressure ratio in key areas of the foundation Select several monitoring points in areas prone to foundation liquefaction (such as below the toe of a slope) and calculate their excess pore water pressure. With initial vertical effective stress The average of the ratios is used to assess the risk of foundation bearing capacity loss.

[0062] ; Constructing a composite limit state function. Given that the three failure modes described above are characteristic of a series system where failure is determined by exceeding any one of the indices, the limit state function is defined using the envelope surface of multiple linear or nonlinear inequalities. To account for the potential interactions between different deformation modes (e.g., ground liquefaction exacerbating settlement and slippage), a convex generalized hyperellipsoidal failure surface is constructed: ; in, , , These are the individual performance thresholds set according to design specifications (e.g., a crest settlement limit of 1.5% and a critical pore pressure ratio of 0.85). Index The coupling coefficient is typically between 4.0 and 8.0. This equation defines a closed surface in three-dimensional state space that is approximately a flat-topped rectangle with rounded corners. When... When the structure is in a failure state, it is determined that the structure is in a failure state; when At that time, the structure is determined to be in a safe state.

[0063] A cumulative damage weakening factor is introduced. This considers that the cumulative damage effect over the duration of a seismic event not only increases the response but also reduces the critical resistance of the structure (e.g., the decrease in friction angle after particle breakage leads to a reduction in critical slip). A time-varying threshold attenuation function is defined. Adjust the static threshold: ; in, The cumulative plastic damage variable calculated in the preceding steps, The sensitivity coefficient is determined by materials experiments. Substituting the modified dynamic threshold into the limit state equation makes the failure surface... It shrinks dynamically as damage accumulates, thus more accurately reflecting the vulnerability of the structure in the later stages of a strong earthquake.

[0064] Numerical dispersion and decision-making implementation. At each time step of the time history analysis, the calculated instantaneous state vector is... Substituting the above Verification is performed within the equations. No complex mesh generation is required; the algebraic value is calculated directly. If several consecutive integration steps (e.g., consecutive 0.5 seconds) occur... If the time is recorded as "First Passage Time", then the specific failure mode combination (such as sedimentation and liquefaction combined failure) is marked as the basic sample data for reliability assessment.

[0065] After constructing the dynamic limit state surface, the first crossing determination of the random trajectory (sub-step S420) is executed as follows: Perform batch time-history simulations. The simulations generated in step S330, which include different phase and amplitude characteristics, will be used to perform batch time-history simulations. A set of corrected load sample sequences (e.g., 50 to 100 sets) is input one by one into the nonlinear finite element model of the guide embankment. Using parallel computing technology, multiple sets of dynamic time history analyses are performed simultaneously, and the corresponding results are output. Group structure state response vector trajectory Each trajectory records the evolution of settlement, slippage, and pore pressure ratio over time.

[0066] Calculate the algebraic safety margin time history. To avoid complex geometric distance calculations, the predefined limit state function is used directly. The value is used to measure the safety margin. Define the first... The sample trajectory at time Algebraic safety margin : ; According to this definition, when When, the structure is in a safe state; when At that time, the structure experiences functional failure or damage.

[0067] Calculate the probability of failure during the first crossing. Iterate through each calculated safety margin time history. To identify whether it is at any time An event occurred where the system crossed the zero axis downwards. The time-varying failure probability of the system was estimated using Direct Monte Carlo Counting. : ; in, This is an indicator function; it takes the value 1 if the condition is met, and 0 otherwise. This formula calculates the result up to time [time value missing]. The proportion of samples that have experienced at least one failed crossing. Compared to analytical approximate solutions based on the crossing rate assumption, this statistical method does not require the assumption of independence of failure events and is more suitable for strongly nonlinear and memory-based geotechnical dynamic systems.

[0068] Evaluate the reliability of the full-probability design. Obtain the cumulative failure probability after the entire seismic motion duration ends. By using the inverse function of the standard normal distribution function, it can be transformed into a global reliability index. : ; like If the reliability exceeds the set target (e.g., 3.7 for a primary guide embankment and 3.2 for a secondary guide embankment), then the current design scheme is deemed to meet the predetermined safety standards even considering the most unfavorable coupling effect of resonant waves and strong earthquakes. Otherwise, based on the statistical distribution of failure modes (e.g., failure mainly caused by liquefaction), targeted engineering reinforcement measures (such as widening the pressure foot or performing foundation treatment) should be taken, and the design should be re-verified.

[0069] After obtaining the initial crossing determination results of the batch samples, in order to more accurately estimate low failure probability events (i.e., rare events), the statistical calculation of failure probability (sub-step S430) is specifically performed as follows: Construct the maximum hazard extreme value sequence. To utilize the classic Generalized Extreme Value (GEV) distribution theory, the minimum safety margin problem is transformed into a maximum hazard problem. For each sample... Define the maximum risk index The negative of its minimum algebraic safety margin over the entire time span: ; like This means that the sample crossed the failure boundary at some point. By traversing all... Using simulated samples, construct a maximum risk sample sequence. .

[0070] Perform parameter fitting for the generalized extreme value distribution. Assume the maximum hazard sequence follows a generalized extreme value distribution, and its cumulative distribution function... It has the following forms: ; in, For position parameters, scale parameter For shape parameters. Probability Weighted Moments (PWM) or Maximum Likelihood Estimation (MLE) are used based on sample sequences. Estimate the three distribution parameters mentioned above. This step constructs a distribution model of the risk index across the entire probability space (including unsampled extreme regions) through statistical inference.

[0071] Calculate the tail-end failure probability. The failure probability of the structure is equivalent to the maximum hazard level. The probability of being greater than or equal to 0. Using the fitted cumulative distribution function, calculate the transcendence probability: ; This analytical calculation method effectively utilizes the tail characteristics of the sample distribution, even when the number of finite element simulations is limited (e.g., ...). Even when the number of direct failure samples is extremely small or even zero, it can still provide a failure probability estimate with a certain degree of confidence based on statistical trends. Reliability results output and confidence level test. The calculated failure probability... Converted into reliability index And output. To evaluate the stability of the statistical results, the Bootstrap resampling method is used: a large subset of samples is drawn with replacement from the original sequence, the above fitting process is repeated, and the results are calculated. The distribution range. If the coefficient of variation (COV) of the result is less than 0.15, the sample size is considered sufficient and the result is valid; otherwise, additional finite element analysis is required to increase the sample size.

[0072] See attached document Figure 2 To comprehensively assess the vulnerability of the guide embankment under different combinations of earthquakes and waves of varying intensities, the Intensity Index (IM) gridded scan (sub-step S510) is performed as follows: Define a multidimensional space for earthquake and wave intensity parameters. Select independent intensity indices that can characterize the potential for damage caused by external environmental excitation. For ground motion, select peak ground acceleration (PGA) or spectral acceleration (GFA). For wave loads, select the effective wave height ( and spectral peak period Construct a two-dimensional or multi-dimensional strength index space that includes all working conditions. Set the scan range to cover conditions ranging from common to extremely rare environments; for example, the PGA range could be set to... Significant wave height range .

[0073] Design a discretized grid sampling scheme. Discretize the continuous intensity space to generate an intensity combination node matrix for analysis. To balance computational cost with the fitting accuracy of the vulnerability curve, a stripes method or incremental dynamic analysis (IDA) point placement strategy is adopted. A denser sampling step size (e.g., [missing information]) is set in the sensitive intensity range where failure is expected (e.g., PGA between 0.3g and 0.6g). In low-intensity or extremely high-intensity sections, a sparse step size is used. This generates a sequence containing... A set of scenarios with discrete intensity combinations .

[0074] Construct the time history of the dynamic input for the target intensity. For each grid node... This generates the corresponding earthquake acceleration time history and wave force time history.

[0075] For ground motion: Samples are selected from a pre-selected set of benchmark strong motion records, and amplitude scaling is used to adjust the PGA of the original records to the target value. .

[0076] For waves: Considering wave nonlinearity and steepness limitations, directly linearly amplifying the wavefront elevation is not recommended. A target spectrum-based reconstruction method is used: based on the target's effective wave height... The target power spectral density function is generated using the JONSWAP or PM spectral formula, along with the corresponding spectral peak period. The corresponding wave time history is then synthesized using the random phase method. Multiple sets (e.g., 10 to 20 sets) of time history samples with different random phases are generated for each intensity node to account for the record-to-record variance at a given intensity.

[0077] To quantitatively describe the probabilistic characteristics of structural engineering requirement parameters (EDP) as a function of bivariate strength indices, the surface function fitting (sub-step S520) including coupling terms is performed as follows: A multivariate probabilistic demand model (PSDM) is constructed. To capture the nonlinear synergistic effect of the combined action of strong earthquakes and large waves, structural engineering demand parameters (EDP) are established, along with ground motion intensity (IM1, e.g., PGA) and wave intensity (IN2, e.g., ...). The analytical relationship between them. Unlike the traditional independent superposition assumption, this embodiment adopts a log-linear regression model that includes cross-coupling terms.

[0078] ; in, The median value of the response at a given intensity; The regression coefficients are to be determined. The coupling coefficient, whose sign and magnitude directly quantify the mutual enhancement or inhibition effects of the two disaster sources when causing structural damage; The model residuals follow a normal distribution with zero mean. The conditional standard deviation in logarithmic space. Perform parameter estimation and significance testing. Using the large-scale "intensity and response" sample dataset obtained in step S510, perform regression analysis on the above model parameters using least squares (OLS) or maximum likelihood estimation. Calculate the statistics and p-values ​​of each regression coefficient. If the coefficients of the coupling terms (such as...) If the p-value of the term is less than the significance level (e.g., 0.05), the term is retained, indicating that the coupling effect between earthquakes and waves is statistically significant; otherwise, the term can be removed, simplifying the model to an uncoupled one. This step ensures the simplicity and physical realism of the model structure.

[0079] Generate the three-dimensional vulnerable surface equation. Define the failure probability based on structural reliability theory. For a structural response exceeding a specified limit state threshold (e.g., the probability of ultimate settlement). Combining the median response model obtained from regression and the total uncertainty, the vulnerability surface function under each ultimate state is constructed: ; in, It is the standard normal cumulative distribution function; This is the uncertainty factor for the structural resistance capacity (usually taken as 0.2 to 0.3). To model the uncertainty coefficient. This equation describes a factor in three-dimensional space. A continuously changing surface.

[0080] Model validation and boundary checks. Plot the fitted vulnerability surface and check its physical plausibility. Focus on checking regions with extremely low values ​​for a single physical quantity (e.g., or The surface morphology of the curve was examined to confirm whether it could asymptotically degenerate into the classic vulnerability curve (S-shaped curve) under a single disaster. Simultaneously, the coefficient of determination was calculated. The root mean square error (RMSE) is used to quantify the goodness of fit of the model. If the fitting accuracy is insufficient (e.g.) If the problem persists, a piecewise regression model or a higher-order response surface model should be used for refitting.

[0081] After the vulnerable surface is constructed, it needs to be transformed into specific risk assessment indicators to guide engineering practice. The results output and risk assessment (sub-step S530) are executed as follows: Generate a multidimensional vulnerability map. Using the calculated vulnerability function, generate a three-dimensional surface plot or a two-dimensional isoprobability contour plot on the computer terminal. This map uses the horizontal axis to represent seismic ground motion intensity (e.g., PGA) and the vertical axis to represent wave intensity (e.g., wave strength ... The color map or Z-axis height represents the probability of structural failure. The map visually reveals the topological shape of the structural safety domain; for example, the degree to which contour lines bulge towards low-intensity regions directly reflects the magnitude of the risk increase caused by coupling effects.

[0082] The annual average failure frequency is calculated based on the total probability formula. To quantify the absolute risk level of the structure, risk integration is performed by combining the site's seismic hazard curve and a long-term statistical distribution model of waves. Considering that strong earthquakes are discrete Poisson processes while the wave environment is a continuous stochastic process, an annual failure frequency calculation model based on the total probability formula is adopted: ; in, The average annual failure frequency of the structure; This represents the long-term probability density function of waves; The absolute value of the derivative of the seismic hazard curve represents the annual probability density of earthquake intensity. The fragile surface function constructed in the preceding steps.

[0083] Implement risk-based design decisions. Calculate the average annual failure frequency. Compared with the permissible risk standards stipulated in the regulations (e.g.) Compare ( / year).

[0084] like If the current guide dike design meets the requirements for multi-hazard defense throughout the entire life cycle, then it is determined that the current design scheme meets the requirements for multi-hazard defense throughout the entire life cycle.

[0085] like If the design fails, optimization is required. By calculating the sensitivity coefficients of the vulnerability function to each design parameter (Sensitivity Analysis), the main controlling factors leading to excessive risk are identified. For example, if the failure probability is found to be most sensitive to the slip threshold value, measures such as increasing the mass of the facing blocks or using blocks with stronger interlocking (e.g., T-shaped blocks) are prioritized; if the failure probability is sensitive to the settlement threshold, foundation reinforcement measures are prioritized. After optimization, the design is re-entered into the above process for verification until the risk standards are met.

[0086] Output structured risk assessment data. The final assessment results will be compiled into a dataset containing key indicators, including: coefficients of the vulnerable surface parametric equations under various extreme states, a table of conditional failure probabilities under different return periods, and the determination of the dominant failure mode (e.g., whether the risk of instability caused by foundation liquefaction exceeds 80% in this marine environment). This data will serve as the core basis for engineering design documents, guiding construction drawing design or the formulation of reinforcement and maintenance strategies for existing structures.

Claims

1. A method for plotting the vulnerability curve of a guide embankment structure based on random load combinations, characterized in that, Includes the following steps: Based on fluid-structure interaction theory, the guide dike structure is discretized into a liquid-containing porous medium system, and a nonlinear finite element analysis model of the deep-water guide dike is constructed. A multidimensional non-stationary stochastic load model was constructed based on the phase coupling mechanism, generating a stochastic dynamic load sample sequence that includes seismic acceleration time history and wave force time history with time-domain coupling. The random dynamic load sample sequence is input into the nonlinear finite element analysis model of the deep-water guide embankment in batches to output the multidimensional dynamic response of the structure. The structural limit state threshold is set according to the cumulative damage evolution law, and a multidimensional dynamic limit state surface is constructed in the generalized state space. The structural failure state is determined based on the relative relationship between the multidimensional dynamic response of the structure and the multidimensional dynamic limit state surface. Based on the failure state of the structure, a regression analysis is performed on the multidimensional dynamic response of the structure to establish a probabilistic demand model containing coupling terms; Based on the probabilistic demand model containing coupling terms and the structural limit state threshold, a vulnerability surface describing the change of failure probability with the seismic intensity index and wave intensity index is plotted.

2. The method for plotting the vulnerability curve of a guide embankment structure based on random load combinations according to claim 1, characterized in that, The steps of constructing a multidimensional non-stationary stochastic load model based on the phase coupling mechanism and generating a stochastic dynamic load sample sequence containing time-domain coupled seismic acceleration time histories and wave force time histories specifically include: The JONSWAP spectral model and the CloughPenzien spectral model are selected as frequency domain references. Intensity modulation functions are introduced to modulate the time domain amplitude of the frequency domain references, thereby constructing the evolution power spectral density function to establish the multidimensional non-stationary random load model with phase coupling mechanism. Generate an independent random phase angle sequence that is uniformly distributed within the interval, define a phase coupling vector that characterizes the peak time delay between waves and earthquakes, and superimpose the phase coupling vector onto the independent random phase angle sequence to form a combined phase parameter that embodies the phase coupling mechanism; Substituting the evolved power spectral density function and the combined phase parameters into the cosine series superposition formula, multiple sets of seismic acceleration time histories and wave force time histories are synthesized and superimposed to form the random dynamic load sample sequence.

3. The method for plotting the vulnerability curve of a guide embankment structure based on random load combinations according to claim 1, characterized in that, The steps of discretizing the guide structure into a liquid-containing porous medium system based on fluid-structure interaction theory and constructing a nonlinear finite element analysis model of the deep-water guide structure specifically include: The rockfill and foundation of the guide embankment structure are discretized using up-format two-dimensional plane strain elements, and an additional mass density based on porosity tortuosity is introduced to construct a solid-phase skeleton mass matrix, thus characterizing the guide embankment structure as a liquid-containing porous media system. A flow-dependent nonlinear damping matrix is ​​defined, and the interaction force between the fluid and the skeleton is updated in real time according to the relative seepage velocity vector of the pore fluid to simulate the nonlinear energy dissipation characteristics caused by the transition of the pore flow from laminar to turbulent flow. Define a time-varying stiffness matrix driven by damage, calculate damage variables based on the cumulative plastic shear strain of the material, and use the damage variables to reduce the tangential stiffness matrix of the structure in real time. Based on the liquid-containing porous media system, a nonlinear finite element model of the guide embankment is constructed by assembling the solid skeleton mass matrix, the flow-dependent nonlinear damping matrix, and the damage-driven time-varying stiffness matrix.

4. The method for plotting the vulnerability curve of a guide embankment structure based on random load combinations according to claim 1, characterized in that, The steps for constructing a multidimensional dynamic limit state surface in the generalized state space specifically include: The settlement of the embankment crest was extracted from the multidimensional dynamic response of the structure and normalized; the maximum slip of the protective layer and the average pore pressure ratio of the key area of ​​the foundation were extracted as independent state components. Define a limit state function, the value of which is determined by the sum of the powers of the ratios of the independent state components to their respective corresponding structural limit state thresholds, minus a constant term. The multidimensional dynamic limit state surface is a convex closed surface formed by points in the three-dimensional state space where the limit state function value is zero.

5. The method for plotting the vulnerability curve of a guide embankment structure based on random load combinations according to claim 4, characterized in that, The steps for setting the structural limit state threshold based on the cumulative damage evolution law specifically include: The initial static threshold is determined according to the design specifications or functional requirements of the guide embankment structure; During dynamic time history analysis, the cumulative plastic damage variable of the material is calculated in real time; Define a decay function that monotonically decreases as the cumulative plastic damage variable increases; The initial static threshold is reduced by the decay function to calculate the structural limit state threshold at the current moment, so that the multidimensional dynamic limit state surface shrinks towards the origin as damage accumulates.

6. The method for plotting the vulnerability curve of a guide embankment structure based on random load combinations according to claim 4, characterized in that, The steps for determining the structural failure state specifically include: The algebraic safety margin is calculated and defined as the negative of the value obtained by substituting the multidimensional dynamic response of the structure into the limit state function expression; The instantaneous value of the algebraic safety margin is less than or equal to zero, which is determined as the structural failure state.

7. The method for plotting the vulnerability curve of a guide embankment structure based on random load combinations according to claim 1, characterized in that, Prior to the step of generating a stochastic dynamic load sample sequence containing time-domain coupled seismic acceleration time histories and wave force time histories, the process further includes determining the combination of values ​​for the seismic ground motion intensity index and the wave intensity index: Peak ground acceleration is selected as the seismic intensity index, and significant wave height is selected as the wave intensity index. The coverage range of the ground motion intensity index and the wave intensity index is set, and multiple intensity combination nodes are generated by using a discretized grid sampling scheme. A denser sampling step size is set in the expected vulnerable and sensitive intensity section. The target intensity corresponding to each of the intensity combination nodes is used to generate the random dynamic load sample sequence.

8. The method for plotting the vulnerability curve of a guide embankment structure based on random load combinations according to claim 1, characterized in that, The specific steps for establishing a probabilistic demand model that includes coupling terms include: The median value of the multidimensional dynamic response of the structure is defined to have a logarithmic linear relationship with the seismic intensity index and the wave intensity index. A cross-coupling term is introduced into the regression model, which is the product of the natural logarithm of the ground motion intensity index and the natural logarithm of the wave intensity index. The multidimensional dynamic response of the structure is analyzed by multiple regression using the least squares method to determine the regression coefficients and establish the probabilistic demand model containing coupling terms.

9. The method for plotting the vulnerability curve of a guide embankment structure based on random load combinations according to claim 8, characterized in that, The specific steps for drawing a vulnerable surface include: Based on structural reliability theory, a vulnerability function in the form of a standard normal cumulative distribution function is constructed; The independent variable of the vulnerability function is determined by dividing the difference between the median logarithm of the response predicted by the probabilistic demand model containing coupling terms and the logarithm of the structural limit state threshold by the combined logarithmic standard deviation. The overall logarithmic standard deviation is formed by combining the residual standard deviation of the regression model, the uncertainty coefficient of the structural resistance capability, and the modeling uncertainty coefficient through square root operations. Substitute the values ​​of the ground motion intensity index and the wave intensity index into the vulnerability function to calculate the corresponding failure probability, generate a spatial geometric shape describing the change of failure probability with the intensity index, and complete the drawing of the vulnerability surface.

10. The method for plotting the vulnerability curve of a guide embankment structure based on random load combinations according to claim 9, characterized in that, The method further includes a step of performing a full probability integral operation on the plotted vulnerable surface, specifically including: Obtain the seismic hazard curve for the target sea area based on the seismic ground intensity index, and the long-term wave probability density function based on the wave intensity index, as environmental input parameters; Using the total probability formula, the annual average failure frequency of the guide embankment structure is obtained by integrating the absolute value of the derivative of the vulnerability function, the seismic hazard curve in the environmental input parameters, and the long-term probability density function of the waves within the full strength space composed of the seismic intensity index and the wave intensity index.

Images