© F. Sun et al., Published by EDP Sciences 2025

1 Introduction

Due to the disturbance caused by complex terrain, mountainous wind fields exhibit drastic wind speed gradients and strong unsteady vortex structures, resulting in highly random wind loads on bridge structures [1,2]. The aerodynamic characteristics of bridges are affected by local turbulence, stream shear, and terrain-induced vortices, making it difficult to precisely describe the lift, drag, and torque coefficients through fixed parameter models [3,4]. Wind-induced vibrations exhibit strong nonlinear coupling effects, making it difficult for traditional wind tunnel tests and numerical simulations to fully capture the dynamic response of bridges [5,6]. Establishing a modeling framework that integrates high-resolution wind field reconstruction and deep learning is of great value in analyzing the aerodynamic evolution laws and wind stability of bridges in complex wind fields.

Current bridge wind stability analysis methods face many challenges under complex mountain wind field conditions. Aerodynamic prediction [7,8] relies on wind tunnel tests and numerical simulations, but wind tunnel experiments are difficult to fully reproduce the characteristics of unsteady atmospheric turbulence. The scale effect affects practical applications; the test cost is high; and data acquisition is limited. Numerical simulation uses the Reynolds-Averaged Navier-Stokes (RANS) method [9,10] or large eddy simulation (LES) [11,12]. The former cannot accurately describe the transient turbulence effect under high Reynolds number conditions, and the latter is computationally intensive and is not suitable for long-term wind-induced vibration analysis of the entire bridge. Existing aerodynamic coefficient models are mostly based on empirical equations or low-dimensional parameterized expressions, lacking modeling of spatial correlation in complex wind fields and are difficult to adapt to the irregular changes in wind fields in mountainous environments [13,14]. Traditional wind load calculations usually use a one-way fluid-solid coupling method, ignoring the reaction effect of bridge structure vibration on the wind field, resulting in insufficient consideration of changes in aerodynamic damping and aerodynamic stiffness [15,16]. The solution to bridge wind resistance response often uses the finite element method [17,18] for structural dynamics analysis. However, in a strong nonlinear turbulent environment, the wind-induced vibration process involves multi-scale energy transfer and fluid-structure coupling characteristics. Existing methods make it difficult to precisely obtain the evolution of vibration modes and nonlinear wind-induced instability phenomena. The aerodynamic prediction model has limited generalization ability under different wind field conditions and cannot adapt to wind speed profiles, vortex shedding patterns, and flow field disturbance characteristics under different terrain conditions [19,20]. Although deep learning-based modeling methods have been applied to the field of wind engineering, most existing neural network architectures focus on local flow field prediction and is difficult to extract cross-scale correlation information of aerodynamic evolution from a global perspective, which limits its application in bridge wind resistance analysis.

For the wind-resistant stability of bridges, scholars have proposed a variety of methods, including large eddy simulation based on CFD (Computational Fluid Dynamics) [21,22], Reynolds time-averaged method, and fluid-solid coupling [23,24] calculation. Wan L conducted an in-depth analysis of the aerodynamic characteristics of grid-type high supports under typhoon conditions through wind tunnel force tests and particle image velocimetry technology, revealing the influence of wind direction, vortex core strength and turbulence on its stability, providing an important basis for the safe design and application of bridges in typhoon-prone areas [25]. Zhang Z analyzed the ice shock wave vibration caused by the collision of bridges, ice, water and air, and compared the semi-arbitrary Lagrangian-Euler and arbitrary Lagrangian-Euler methods, and found that the semi-arbitrary Lagrangian-Euler method was more suitable for complex fluid-solid coupling analysis. The results verified the applicability of the ice shock vibration theory and revealed the significant influence of velocity and thickness on vibration [26]. However, the LES method has a very high computational cost and is difficult to apply to long-term dynamic analysis of the entire bridge scale. The RANS method cannot accurately obtain transient wind load effects due to the limitations of the turbulent closed model. Although the fluid-solid coupling method can realize the calculation of wind-induced response, it has a high reliance on the spatiotemporal resolution of the wind field and is difficult to apply to unsteady complex wind fields. The application of deep learning in fluid mechanics modeling is increasing. Methods such as convolutional neural networks (CNNs) [27,28] and long short-term memory (LSTM) networks [29,30] have made certain progress in wind speed prediction and aerodynamic modeling. However, these methods are usually difficult to handle high-dimensional complex wind field data and lack cross-scale generalization capabilities, resulting in reduced prediction accuracy under mountain wind field conditions. Therefore, existing research still faces the problem of balancing computational efficiency and modeling accuracy in bridge wind stability analysis.

In recent years, fluid dynamics modeling methods based on neural operators such as FNO [31,32] and DeepONet [33,34] have emerged to address the problems of wind field complexity and computational cost. They achieve efficient flow field prediction by learning the mapping relationship between wind field variables and have strong generalization capabilities under unstructured data. Zhang K developed a deep learning model based on FNO to solve the two-dimensional oil/water two-phase flow partial differential equation problem. By improving the network structure and physical constraints, the model showed the potential of neural networks to replace traditional numerical simulation in terms of accuracy and generalization [35]. FNO extracts cross-scale features through frequency domain transformation and has shown superior performance in turbulence modeling. At the same time, the diffusion model is widely used in turbulence reconstruction and wind field simulation due to its advantages in generating high-fidelity data [36,37]. Combining the advantages of the diffusion model and FNO, three-dimensional wind field data containing complex terrain disturbances can be generated in a mountain wind field environment, and high-precision modeling of the three-force coefficients of the bridge can be achieved. At present, related research mainly focuses on single fluid dynamic characteristics and lacks bidirectional coupling modeling of the dynamic response of bridge structures. This paper constructs a method that combines Diffusion and FNO-ODE hybrid architecture to achieve a unified calculation framework for wind load prediction and bridge structure dynamic response.

This study constructs a three-force coefficient prediction framework based on deep learning to improve the accuracy and computational efficiency of wind stability analysis of bridges in mountain wind fields. Combining the diffusion model with FNO, end-to-end modeling from high-dimensional wind field data to aerodynamic characteristics is achieved. The diffusion model is used to reconstruct high-resolution turbulent wind fields. Physical constraints are applied on the basis of large eddy simulation data to improve the authenticity of data generation, so that the prediction of aerodynamic characteristics can cover the transient wind field structure under complex terrain disturbances. The FNO network uses frequency domain transformation to obtain the spatial correlation of the aerodynamic evolution process. Compared with conventional convolutional neural networks and long short-term memory networks, it has stronger cross-scale generalization capabilities and can effectively extract key features such as vortex separation and reattachment. Another contribution of the study is the construction of FNO-ODE, which embeds the wind load time series data predicted by deep learning into the Newmark-β solution process to form a high-precision calculation system that can describe the wind-structure coupling effect. Compared with the static aerodynamic fitting method, this framework can simulate the nonlinear dynamic response of bridges under different vibration modes and maintain stable calculation accuracy under high Reynolds number conditions. The Hilbert spectrum analysis strategy is applied to realize the automatic feature extraction of the bridge buffeting response time history, identify the wind-induced unstable state, and construct the probability distribution model of the unstable wind speed to improve the intelligent level of wind resistance safety assessment. This method can adapt to different wind field conditions and is suitable for the aerodynamic performance optimization design of large-span bridges in mountainous areas, providing new modeling ideas and calculation tools for wind resistance engineering.

2 Materials and methods

2.1 Construction of mountain wind field generation module and design of aerodynamic feature decoupling network

Figure 1 shows a wind stability evaluation framework for mountain bridges that integrates physical constraint deep learning and aeroelastic analysis. The core architecture consists of four coupling modules. (1) The mountain wind field generation module uses a physical constraint diffusion model to construct a three-dimensional fluctuating wind speed field based on WRF-LES data. The Navier-Stokes equation residual constraint is used to ensure mass conservation, and a 0.1D resolution wind field is generated by combining terrain curvature feature interpolation. (2) The aerodynamic feature decoupling module designs an eight-layer Fourier integral operator to decompose the multi-scale vorticity field in the wave number space, and uses the wall shear stress zero value tracking to quantify the separation bubble parameters. (3) The dynamic coupling solution module establishes an FNO-ODE hybrid model, maps the aerodynamic force to the Newmark-β time-history integral, compresses the structural freedom degree through the modal truncation technology, and realizes millisecond-level bidirectional coupling simulation. (4) The stability criterion module detects the spectral entropy mutation based on the Hilbert time-frequency spectrum attention mechanism to identify the divergence starting point, and constructs the critical wind speed probability distribution feedback to optimize the wind field generation. The architecture realizes multi-scale coupled analysis of bridge aerodynamic stability under complex terrain through Fourier operator frequency domain analysis, physical constraint embedding, and closed-loop feedback mechanism, providing a high-precision probabilistic evaluation tool for wind-resistant design.

Fig. 1 Evaluation framework.

2.1.1 Construction of physical constraint diffusion model

The turbulent characteristics of the WRF large eddy simulation dataset are embedded in the conditional diffusion framework to construct an implicit mapping relationship between terrain disturbance and wind speed field. The conditional diffusion process is defined as the reverse derivation of the Markov chain:

$$q\left(U_{t-1}|U_t,C\right)=\mathcal{N}\left(U_{t-1};{\mu }_{\theta }\left(U_t,t,C\right),{\sigma }_t^{2}I\right).$$(1)

Among them, C ∈ ℝ^H×W×3 is the three-dimensional conditional tensor of terrain elevation, roughness, and pressure gradient, and U_t is the wind speed field tensor of the diffusion state at time t (Unit: m/s), ${\sigma }_t^2$ is the diffusion noise variance, unit: m²/s², 𝒩 is a Gaussian distribution, I is the identity matrix. The optimization objective function of the diffusion model parameter θ integrates the physical conservation constraint term:

$$ℒ={\mathbb{E}}_{t,U_0,ϵ}\left[\parallel ϵ-{ϵ}_{\theta }{\left(U_t,t,C\right)\parallel }_2^2\right]+\lambda {\displaystyle {{\displaystyle \sum }}_{i=1}^3}\parallel \nabla {\cdot \left(\rho U_i\right)\parallel }_{L^2}.$$(2)

The second term forces the residual of the mass conservation equation to be minimized, and ρ is the air density, Unit: kg/m³, ϵ is the noise vector, unit: m/s, $\mathcal{L}$ is the total loss function, λ is the weight coefficient of the physical constraint term, $\parallel \nabla \cdot {\left(\rho U_i\right)\parallel }_{L^2}$ is the residual of the continuity equation. The U-Net backbone network is constructed using spectral normalization convolution blocks. The input layer embeds terrain curvature features, and the output layer correlates with the prediction of the eddy viscosity coefficient:

$${\nu }_t=\frac{{\left({\mu }_{\theta }^{\left(x\right)}\right)}^2+{\left({\mu }_{\theta }^{\left(y\right)}\right)}^2+{\left({\mu }_{\theta }^{\left(z\right)}\right)}^2}{3\left|S\right|}$$(3)

S is the strain rate tensor, and μ_θ is the mean velocity gradient predicted by the network. This design makes the generated wind speed field meet the boundary layer momentum transport characteristics. A phased optimization strategy is adopted in the diffusion model training: In the first stage, only the data reconstruction loss is optimized, the learning rate is set to 1 × 10^‒4, and the basic turbulence statistical characteristics are established; in the second stage, the mass conservation physical constraint is introduced, the learning rate is reduced to 5 × 10^‒5, and the terrain curvature characteristics are gradually integrated; in the third stage, the generator parameters are fixed, the decoder layer is fine-tuned, the learning rate is 1 × 10^‒5, and the spectral normalization technique is used to stabilize the training. The Adam optimizer is used in each stage, the batch size is kept at 64, and the learning rate is adjusted by cosine annealing. The network structure adopts a multi-scale feature extraction module, and applies jump connections between the encoder and the decoder to ensure the feature transfer of turbulence structures of different scales. In the loss function design, in addition to the physical constraint term, the spectral energy distribution matching term is also applied to ensure that the turbulence spectrum characteristics of the generated wind speed field are consistent with the WRF data. The spectral energy matching term is achieved by comparing the slope of the inertial sub-region of the generated data with that of the WRF data, using the Kolmogorov spectrum fitting error minimization constraint. Table 1 shows the comparison of the statistical characteristics of the generated wind speed field.

Table 1 shows the statistical characteristics of the wind speed field generated based on the physically constrained diffusion model and the original WRF large eddy simulation data. This comparison covers the key turbulence characteristic indicators of the wind speed field, including core parameters such as average wind speed, turbulence intensity, and Reynolds stress, which are used to verify the physical rationality of the generated data. The consistency between the generated data and the reference data was systematically evaluated through multi-dimensional indicators, with special attention paid to the turbulence characteristics unique to mountain wind fields. The comparison of turbulent kinetic energy and energy spectrum slope verifies the accuracy of the model in modeling the turbulent energy series process, while the comparison of integral length scales reflects the model's ability to capture large-scale turbulent structures. The error rates of various indicators are all controlled at a low level, indicating that the generated data can accurately reflect the impact of terrain disturbances on the wind field while maintaining statistical characteristics. This quantitative comparison provides a reliable data basis for subsequent aerodynamic calculations.

2.1.2 Generation of multi-resolution wind speed field

The latent space interpolation operation uses the non-uniform B-spline surface parameterization method to establish a nonlinear mapping between terrain features and turbulence scales. The generation process of the latent vector z ∈ ℝ²⁵⁶ is defined as:

$$z=f_{\varphi }\left(C\right)\oplus g_{\psi }\left(W_{LES}\right)$$(4)

f_ϕ is the terrain encoder; g_ψ is the WRF large eddy simulation data encoder; ⊕ which represents feature splicing. The interpolation kernel function uses the exponential radial basis function:

$$K\left(z_i,z_j\right)=\text{exp}\left(-\frac{\parallel z_i-{z_j\parallel }_2^2}{2{\sigma }_z^2}\right)$$(5)

K (z_i, z_j) is the dimensionless value of the interpolation kernel function, which measures the similarity or influence weight between spatial points z_i and z_j, σ is the width parameter of the radial basis function. Based on this, a three-dimensional wind speed field super-resolution operator is constructed:

$$U_{HR}=\mathcal{G}\left(U_{LR}\uparrow +{\displaystyle {{\displaystyle \sum }}_{k=1}^N}{\alpha }_{k}K\left(z,z_k\right)\cdot W_k\right)$$(6)

↑ Represents bilinear interpolation upsampling, and W_k is the turbulence base mode stored in the training phase, U_HR, U_LR are high/low resolution wind speed fields, in m/s, and α_k is the weight coefficient (dimensionless), $\mathcal{G}$ is the super-resolution operator, N is the number of fundamental modes of turbulence. This process makes the turbulent kinetic energy spectrum of the generated wind speed field satisfy the −5/3 scaling law in the inertial sub-region; the root mean square error of the pulsating wind speed is less than 0.15 m/s; the spatial resolution reaches 0.1D (beam height scale). The terrain disturbance effect is realized through the pressure gradient correction term ΔP = h (∇ z) ⊗ U_HR, and h is the mountain flow function modeled by the multi-layer perceptron. A progressive generation strategy is used in multi-resolution generation to generate a large-scale average wind speed field, which is then gradually refined to a small-scale turbulent structure. An adaptive convolution kernel is applied in the network design to dynamically adjust the size and shape of the convolution kernel according to the terrain characteristics to better obtain the flow field distortion caused by the terrain. A local weighted strategy is used in interpolation to dynamically adjust the interpolation weight according to the terrain characteristics and wind speed field characteristics.

Figure 2 shows the relationship between terrain height wind speed field and turbulent kinetic energy spectrum.

Figure 2a compares the distribution characteristics of low-resolution wind speed field and high-resolution wind speed field. The low-resolution wind speed field is generated by superimposing random noise with a sine function, showing obvious fluctuations and large noise, and cannot precisely obtain small-scale turbulent structures. The high-resolution wind speed field is generated by superimposing a finer sine function with smaller noise, and its fluctuations are smoother. The noise is significantly reduced, and it can better reflect the local influence of terrain on wind speed. The high-resolution wind speed field shows more detailed wind speed changes in areas with significant terrain changes (such as near z = 3 and z = 7), while the low-resolution wind speed field appears rough and discontinuous in these areas. This comparison verifies the effectiveness of the diffusion model in generating high-fidelity wind speed fields, especially in obtaining local wind speed changes caused by terrain disturbances.

Figure 2b shows the distribution characteristics of the turbulent kinetic energy spectrum of the high-resolution wind speed field with frequency. The energy spectrum is calculated by the Fourier transform of the high-resolution wind speed field. The energy in the low-frequency region is higher, and the energy gradually decays with increasing frequency. At a frequency of around 1 Hz, the energy spectrum shows an obvious −5/3 scaling law feature (marked by the gray dotted line), which is a typical feature of the turbulent inertial sub-region. The high energy value in the low-frequency region reflects the influence of terrain disturbances on large-scale turbulent structures, while the energy attenuation in the high-frequency region indicates that the energy distribution of small-scale turbulent pulsations is consistent with theoretical expectations. The −5/3 scaling law characteristics of the energy spectrum in the figure verify the physical consistency of the generated wind speed field in the turbulent statistical characteristics, indicating that the diffusion model can accurately capture the energy transfer process of turbulence and provide reliable input data for subsequent wind resistance stability analysis.

Fig. 2 Relationship between terrain height and wind speed field and turbulent kinetic energy spectrum. (a). Relationship between terrain height and wind speed field. (b). Turbulent kinetic energy spectrum.

2.1.3 Fourier multi-scale feature extraction

Eight-layer cascaded Fourier integral operator network is constructed to establish an implicit mapping between the flow field velocity tensor and the vortex evolution mode. Table 2 shows the configuration of the Fourier integral network layer parameters. The l th layer Fourier integral kernel operation is defined as:

$${\mathcal{F}}_l\left(u\right)\left(k\right)={\displaystyle \underset{D}{{\displaystyle \int }}}{e}^{-2\pi i⟨k,x⟩}{W}_l\left(k\right)\cdot u\left(x\right)dx$$(7)

W_l (k) is the learnable frequency domain parameter matrix, and k is the wave number vector, ${\mathcal{F}}_l$ (u) (k) is the output of the Fourier integral operator at wave number k in the lth layer (frequency domain characteristics), u (x) is the input three-dimensional velocity field (spatial domain), D is the spatial integration domain. The network input is the three-dimensional velocity field slice u in the range of 0.5D to 5D around the main beam, and multi-scale decomposition is achieved through cascaded truncated modes:

$$u_m={\displaystyle {{\displaystyle \sum }}_{\left|k\right|\le {\kappa }_m}}{\mathcal{F}}^{-1}\left({\mathcal{F}}_l\left(u\right)\odot M_m\left(k\right)\right)$$(8)

M_m (k) is the bandpass filter of the mth scale band; κ_m = 2^mκ₀ is the cutoff wave number; κ₀ corresponds to the minimum resolution scale of 0.5D, u_m is the velocity field after multi-scale decomposition of the mth scale band, and the unit is m/s, κ_m is the cutoff wave number of the scale band, in rad/m, ${\mathcal{F}}^{\mathrm{-1}}$ is the inverse Fourier transform operator, ⊙ is element-by-element multiplication, M_m (k) is a bandpass filter with a scale band. The network output layer h integrates the features of each scale:

$$h=\sigma \left({\displaystyle {{\displaystyle \sum }}_{m=1}^8}{\beta }_m\cdot {\text{Conv}}_{1\times 1}\left(u_m\right)\right)$$(9)

β_m is the adaptive weight coefficient, which is dynamically adjusted through the gating mechanism, Conv_1×1 is a 1 × 1 convolution operation used to adjust the channel dimension. The core of the adjustment is to calculate the energy norm of each scale feature in real-time as the input signal of the network; secondly, it performs nonlinear transformations on the features using learnable weight matrices, generating weight coefficients using the Sigmoid function for normalization. This mechanism enables the network to automatically enhance the contribution of key scale features and suppress noise interference based on the local characteristics of the input flow field. In a specific implementation, gating units receive features from all scales, achieving collaborative optimization through cross-scale information interaction. This architecture effectively separates the large-scale coherent structure of the boundary layer separation bubble and the small-scale turbulent pulsation in the wake region.

In the FNO network design, the selection of an eight-layer cascade structure is mainly based on the balance between the physical characteristics of multi-scale turbulence characteristics and computational efficiency. The theoretical basis is that the turbulent energy level series process in mountainous wind fields usually covers 3–5 orders of magnitude (from integral scale to dissipative scale), and the wave number range of the eight-layer network can cover the key spatial scale of 0.5D∼5D (corresponding to the separation bubble, tail vortex and other structures of the bridge flow). When there are less than six layers, small-scale vortices cannot be resolved, and when there are more than ten layers, the accuracy is slightly improved, but the amount of calculation increases significantly. This design avoids the problem of over-parameterization while ensuring the accuracy of the turbulence spectrum in the inertial sub-zone.

2.1.4 Quantitative modeling of vortex evolution

The characteristic constraint term of the vorticity transport equation is defined as:

$$\frac{\partial \omega }{\partial t}+\left(u\cdot \nabla \right)\omega =\left(\omega \cdot \nabla \right)u+\nu {\nabla }^2\omega$$(10)

ν is the dynamic viscosity coefficient, The network hidden state establishes a two-way coupling with the vortex field ω, and constructs the vortex center detection operator:

$$\mathcal{C}\left(x\right)=\text{arg}{\displaystyle \underset{x^{\prime }\in \mathrm{\Omega }}{\text{max}}}\parallel {{\displaystyle {\displaystyle \int }{B}_r\left(x^{\prime }\right)}}_{}\omega \left(x\right)dx{\parallel }_2$$(11)

B_r (x^′) is a spherical domain with a radius of r = 0.2D, $\mathcal{C}\left(x\right)$ is the vortex center detection operator, which is used to locate the spatial coordinates of the vortex core, Ω is the computational domain. The reattachment length of the detached bubble is tracked by the zero point of the wall shear stress:

$$L_r=\text{max}\left\{x\right|{\tau }_w\left(x\right)=0,x\in {\mathrm{\Gamma }}_{\text{wall}}\}$$(12)

L_r is the separation bubble reattachment length, that is, the maximum distance the flow reattaches to the wall after separation, in meters, τ_w (x) is the wall shear stress, Γ_wall is the wall boundary, The network output layer is designed with a dual-branch structure: branch 1 uses a three-dimensional deformable convolution kernel to extract the vortex topological feature ϕ_v, and Branch 2 compresses the wall pressure pulsation signal ϕ_p through the spectral pooling layer. The final prediction value ${{\displaystyle \hat{L}}}_r$ is generated by the feature fusion layer:

$${{\displaystyle \hat{L}}}_r=\text{MLP}\left(\left[{\varphi }_v\oplus {\varphi }_p\right]\right).$$(13)

The vorticity conservation regularization term is applied in the training stage:

$${\mathcal{L}}_{\omega }={\displaystyle {{\displaystyle \sum }}_{t=1}^T}\frac{\partial {\omega }_t}{\partial t}-{\mathcal{N}}_{\omega }{\left(u_t\right)}_{H^1}$$(14)

${\mathcal{N}}_{\omega }$ is the right-hand side term of the vorticity transport equation, ${\mathcal{L}}_{\omega }$ is the vorticity conservation regularization loss term, which is used to constrain the vorticity evolution predicted by the network to conform to physical laws, and the H¹ Sobolev norm enhances the smoothness of the solution. The vorticity conservation regularization term, by enforcing the physical constraints of the vorticity transport equation on the network, can significantly enhance training stability and prevent non-physical phenomena. This constraint ensures that the predicted vortex evolution process maintains spatiotemporal continuity, improving the prediction accuracy of key features such as the reattachment length of separated bubbles.

The network parameter update adopts the quasi-Newton algorithm accelerated by Jacobi preconditioning, and the Hessian matrix is approximated as H_k:

$$H_k=J_k^T{J}_k+\mu {\displaystyle {{\displaystyle \sum }}_{i=1}^8}\frac{{\partial }^2{ℒ}_{\omega }}{\partial W_l^{\left(i\right)2}}$$(15)

J_k is the Jacobian matrix of the loss function ${\mathcal{L}}_{\omega }$ with respect to the network parameters, $W_l^{\left(i\right)2}$ is the trainable parameter of the Fourier integral operator network, This design ensures the spatiotemporal continuity of the vortex evolution law and realizes the sub-grid scale analysis of the separation bubble parameters.

2.2 Dynamic response coupling solvern and intelligent calculation of stability criteria

2.2.1 FNO-ODE hybrid architecture design

A hybrid solution framework that couples Fourier neural operators with ordinary differential equations is constructed to realize real-time two-way coupling of aerodynamic forces and structural vibrations, and the state space equation of the bridge dynamic system is defined as:

$$\frac{dX\left(t\right)}{\mathrm{d}t}=f\left(X\left(t\right),F_a\left(t\right)\right)$$(16)

$X\left(t\right)={\left[u,{\displaystyle \dot{u}}\right]}^T\in {ℝ}^{2n}$ is the state vector of the structural displacement u and velocity ${\displaystyle \dot{u}}$, and F_a (t) is the three components of aerodynamic force . The FNO module maps the flow field characteristics to the aerodynamic force:(lift, drag, and torque)

$$F_a\left(t\right)={\mathcal{G}}_{\theta }\left(U\left(x,t\right)\right)$$(17)

${\mathcal{G}}_{\theta }$ is the parameterized Fourier integral operator, and U (x, t) is the velocity distribution of the flow field around the bridge. The ODE solver uses the implicit Newmark-β format:

$$M{{\displaystyle \ddot{u}}}_{n+1}+C{{\displaystyle \dot{u}}}_{n+1}+K{u}_{n+1}=F_{a,n+1}$$(18)

u_n+1 is the node displacement vector (at time n + 1), unit: m, $u_{n+1}$ is the node velocity vector, unit: m/s, ${{\displaystyle \ddot{u}}}_{n+1}$ is the node acceleration vector, unit: m/s². The time discretization format is:

$$\{\begin{array}{c}\hfill u_{n+1}=u_n+\mathrm{\Delta }t{{\displaystyle \dot{u}}}_n+\frac{\mathrm{\Delta }{t}^2}{2}\left[\left(1-2\beta \right){{\displaystyle \ddot{u}}}_n+2\beta {{\displaystyle \ddot{u}}}_{n+1}\right]\hfill \\ \hfill {{\displaystyle \dot{u}}}_{n+1}={{\displaystyle \dot{u}}}_n+\mathrm{\Delta }t\left[\left(1-\gamma \right){{\displaystyle \ddot{u}}}_n+\gamma {{\displaystyle \ddot{u}}}_{n+1}\right]\hfill \end{array}$$(19)

β = 0.25 and γ = 0.5 ensure unconditional stability. The FNO-ODE coupling is implemented by the Jacobi-Free Newton Krylov method:

$$\mathcal{J}\left(X_{n+1}\right)\mathrm{\Delta }{X}_{n+1}=-\mathcal{R}\left(X_{n+1}\right)$$(20)

$\mathcal{J}$ is the approximate Jacobian matrix, and $\mathcal{R}$ is the residual function. Compared to the traditional Newton-Raphson method, the Jacobi-Free Newton-Krylov method significantly enhances computational efficiency in solving FNO-ODE coupling problems. This method does not require the explicit construction of a Jacobian matrix; instead, it uses vector products to approximate the directional derivatives, which reduces storage and computational costs. It is particularly suitable for high-dimensional nonlinear systems. The convergence rate is comparable to traditional methods, but each iteration's computational load is significantly reduced, thus improving overall solution efficiency. The design compresses the time step to 0.001 s to meet the accuracy requirements of aeroelastic coupling. The time step is set to 0.001 s mainly based on the numerical calculation requirements of the fluid-structure coupling system. This step first meets the CFL (Courant-Friedrichs-Lewy) condition to ensure the stability of the fluid calculation, where the minimum grid size corresponds to one-tenth of the bridge beam height; secondly, it meets the unconditional stability requirements of the Newmark-β method to ensure the convergence of the structural dynamic response calculation. This time step can accurately capture the high-frequency vortex shedding phenomenon and the low-frequency bridge vibration characteristics at the same time, and realize the accurate coupling calculation of wind load and structural vibration.

In the FNO-ODE coupling solution, the Jacobi-Free Newton Krylov method has significant efficiency advantages over the traditional Newton-Raphson iteration. This method does not explicitly construct the Jacobian matrix, but uses vector products to approximate directional derivatives, avoiding the high cost of calculating and storing large dense Jacobian matrices in traditional methods. Especially for high-dimensional nonlinear coupling systems, Jacobi-Free Newton Krylov can greatly reduce the amount of calculation for each iteration while maintaining a similar convergence speed, which significantly improves the overall solution efficiency while maintaining the calculation accuracy.

2.2.2 Bidirectional coupling iterative solution

The bidirectional coupling iterative format $F_{a,n+1}^{\left(k\right)}$and $X_{n+1}^{\left(k\right)}$ of aerodynamic force and structural vibration is established:

$$\{\begin{array}{c}\hfill F_{a,n+1}^{\left(k\right)}={𝒢}_{\theta }\left(U\left(x,t_n+\mathrm{\Delta }t;X_{n+1}^{\left(k-1\right)}\right)\right)\hfill \\ \hfill X_{n+1}^{\left(k\right)}=\mathrm{Newmark}-\beta \left(F_{a,n+1}^{\left(k\right)}\right)\hfill \end{array}$$(21)

k is the number of iteration steps, and the convergence condition is: $\frac{{\parallel X}_{n+1}^{\left(k\right)}-{X_{n+1}^{\left(k-1\right)}\parallel }_2}{\parallel X_{n+1}^{\left(k\right)}\parallel 2}<ϵ$, ϵ = 10^‒6. The FNO module adopts an adaptive time-step strategy:

$$\mathrm{\Delta }{t}_{\text{FNO}}=\text{min}\left(0.001,\frac{C_{FL}}{\text{max}\left(\parallel u{\parallel }_{\infty }\right)}\right)$$(22)

Δt_FNO is the module adaptive time step, The structural solver applies the modal truncation technology to reduce the system's degrees of freedom from n to m:

$${\mathrm{\Phi }}^{T}M\mathrm{\Phi }{\displaystyle \ddot{q}}+{\mathrm{\Phi }}^{T}C\mathrm{\Phi }{\displaystyle \dot{q}}+{\mathrm{\Phi }}^{T}K\mathrm{\Phi }q={\mathrm{\Phi }}^T{F}_a$$(23)

Φ is the first m-order modal matrix, and q is the generalized coordinate. This reduced-order model maintains the calculation accuracy while significantly improving the efficiency.

2.2.3 Nonlinear coupling stability

The Lyapunov index λ_L of the coupled system is defined as:

$${\lambda }_L={\displaystyle \underset{t\to \infty }{\text{lim}}}\frac{1}{t}\text{log}\frac{\parallel \delta X\left(t\right)\parallel }{\parallel \delta X\left(0\right)\parallel }$$(24)

δX is the state disturbance vector. The QR decomposition method is used to calculate λ_L in real time to monitor the system stability. The aerodynamic damping ratio ζ_a is extracted by the Hilbert transform:

$${\zeta }_a=\frac{1}{2{\omega }_0}\frac{d}{\mathrm{d}t}\text{log}A\left(t\right)$$(25)

A (t) is the vibration amplitude envelope, and ω₀ is the natural frequency. This index is used to evaluate the influence of the aeroelastic effect on structural stability. The solver adopts a mixed precision calculation strategy; FP16 is used for aerodynamic prediction; FP32 is used for structural response solution, which ensures accuracy while improving calculation efficiency.

2.2.4 Hilbert spectrum feature attention modeling

A time-frequency joint analysis framework is constructed, and the Hilbert-Huang transform of the buffeting response time history signal ℋ (t, ω) is defined as:

$$\mathscr{H}\left(t,\omega \right)=Re[\begin{array}{c}\hfill {\displaystyle {\int }_{-\infty }^{\infty }}y\left(\tau \right)h\left(\tau -t\right)e^{-i\omega \left(\tau -t\right)}d\tau \hfill \end{array}]$$(26)

In the equation, h is the adaptive window function, and the window length T_w (t) is determined by the local feature scale:

$$T_w\left(t\right)=\frac{2\pi }{\text{max}\left(\frac{d}{\mathrm{d}t}\varphi \left(t\right)\right)}$$(27)

ϕ (t) is the instantaneous phase function. A multi-head attention mechanism is designed to extract the correlation of spectral features, and the kth head attention weight ${\alpha }_{ij}^{\left(k\right)}$ is calculated as:

$${\alpha }_{ij}^{\left(k\right)}=\frac{\text{exp}\left(Q_k{\left(t_i\right)}^T{K}_k\left(t_j\right)/\sqrt{d_k}\right)}{{{\displaystyle \sum }}_{m=1}^N\text{exp}\left(Q_k{\left(t_i\right)}^T{K}_k\left(t_m\right)/\sqrt{d_k}\right)}$$(28)

Q_k and K_k are the query vector and the key vector, and d_k is the feature dimension. The output of the spectral feature fusion layer S (t_i) is:

$$S\left(t_i\right)={\displaystyle {{\displaystyle \sum }}_{k=1}^H}{\displaystyle {{\displaystyle \sum }}_{j=1}^N}{\alpha }_{ij}^{\left(k\right)}{V}_k\left(t_j\right)\odot W_{\omega }\left({\omega }_j\right)$$(29)

V_k is the value vector, and ω is the frequency weight function, satisfying ${\displaystyle \underset{\mathrm{\Omega }}{{\displaystyle \int }}}{W}_{\omega }\left(\omega \right)d\omega =1$. This architecture effectively captures the migration law of the energy concentration frequency band.

2.2.5 Divergent vibration starting point detection

The spectral entropy mutation detection operator is defined as:

$$\mathrm{\Delta }ℰ\left(t\right)=\mid \frac{\partial }{\partial t}\begin{array}{c}\hfill \left(-{\displaystyle {{\displaystyle \sum }}_{\omega }}{\mathscr{H}}^2\left(t,\omega \right)\text{ln}{\mathscr{H}}^2\left(t,\omega \right)\right)\hfill \end{array}\mid$$(30)

Δε (t) is the time rate of change of spectral entropy, which is used to detect sudden changes in energy distribution (the starting point of divergent vibration), ${\mathcal{H}}^2$ (t, ω) is the Hilbert time spectrum, The threshold for spectral entropy mutation detection is set based on the statistical characteristics analysis of Hilbert spectra. The equilibrium point is determined through the optimization of receiver operating characteristic curves, ensuring that the true positive rate in system stability assessment exceeds 95% while keeping the false positive rate below 8%. The threshold setting method of the spectral entropy mutation detection operator can determine the equilibrium point by optimizing the ROC curve, reducing the dependence on empirical parameters. This method combines the theoretical distribution of spectral entropy with the requirements of engineering safety margin, and only needs a small number of calibration tests to verify its adaptability, which has low empirical dependence and high stability evaluation accuracy. This method integrates the theoretical distribution characteristics of spectral entropy with actual engineering safety margin requirements, not entirely relying on empirical parameters, but requiring a few calibration tests for different bridge types to verify adaptability.

A bidirectional LSTM network is constructed to predict the probability of the divergent point:

$$p_d\left(t\right)=\sigma \left(W_d\left[S\left(t\right)\oplus \mathrm{\Delta }ℰ\left(t\right)\right]+b_d\right)$$(31)

p_d (t) is the probability of divergent vibration, W_d is the weight matrix, b_d is the bias vector, S (t) is the spectral feature vector output by the attention mechanism. The interval loss function ${\mathcal{L}}_d$ is used in the training phase:

$${ℒ}_d={\displaystyle {{\displaystyle \sum }}_{t=t_c-\mathrm{\Delta }t}^{t_c+\mathrm{\Delta }t}}\left[{\lambda }_1\text{max}\left(0,\gamma -p_d\left(t\right)\right)+{\lambda }_2\text{max}\left(0,p_d\left(t\right)-\gamma \right)\right]$$(32)

t_c is the actual divergence moment, and γ = 0.5 is the classification threshold, p_d (t) is the divergence probability predicted by LSTM. The detection algorithm simultaneously monitors the instantaneous energy growth rate η_e (t):

$${\eta }_e\left(t\right)=\frac{1}{T}{\displaystyle {{\displaystyle \int }}_{t-T}^t}\frac{\partial \parallel \mathscr{H}\left(\tau ,\omega \right){\parallel }_2^2}{\partial \tau }d\tau$$(33)

When η_e (t)>κ_e lasts for more than 3 cycles, it is determined that the system enters the divergent state.

2.2.6 Estimation of the probability distribution of unstable wind speed

A non-parametric conditional probability model is established:

$$P\left(V_{cr}\le v\mid \mathrm{\Theta }\right)=\frac{1}{Z}{\displaystyle {{\displaystyle \sum }}_{i=1}^M}K\left(\frac{v-V_i}{h}\right)\cdot {\displaystyle {{\displaystyle \prod }}_{j=1}^J}\delta \left({\theta }_j-{\mathrm{\Theta }}_j\right)$$(34)

Θ is the environmental parameters such as turbulence intensity and incoming flow angle of attack, P (V_cr ≤ v∣Θ) is the probability that the critical wind speed does not exceed v under condition Θ; K is the Epanechnikov kernel function; the bandwidth is 0.1${\displaystyle \stackrel{\leftharpoonup }{V}}$. The sequential Monte Carlo method is used to update the probability distribution $w_t^{\left(i\right)}$:

$$w_t^{\left(i\right)}=w_{t-1}^{\left(i\right)}\cdot \frac{p\left(y_t\mid V_{cr}^{\left(i\right)}\right)p\left(V_{cr}^{\left(i\right)}\mid {\mathrm{\Theta }}_t\right)}{\pi \left(V_{cr}^{\left(i\right)}\mid {\mathrm{\Theta }}_t\right)}$$(35)

$p\left(y_t|V_{cr}^{\left(i\right)}\right)$ is the observation likelihood (how well the actual vibration response matches the prediction), $p\left(V_{cr}^{\left(i\right)}|{\mathrm{\Theta }}_t\right)$ is the prior probability (historical data distribution),The number of valid particles $N_{\text{eff}}=1/{\displaystyle {{\displaystyle \sum }}_{i=1}^N}{\left(w_t^{\left(i\right)}\right)}^2>N_{\text{th}}$ is retained in the resampling phase. The final output is the confidence interval (CI) of the unstable wind speed:

$$C{I}_{95\%}=\left[{{\displaystyle \hat{V}}}_{cr}-1.96{{\displaystyle \hat{\sigma }}}_V,{{\displaystyle \hat{V}}}_{cr}+1.96{{\displaystyle \hat{\sigma }}}_V\right]$$(36)

${{\displaystyle \hat{\sigma }}}_V$ is the standard deviation of the kernel density estimate, ${{\displaystyle \hat{V}}}_{cr}$ is the estimated critical wind speed. The model dynamically adjusts the probability distribution through an online learning mechanism to adapt to the non-stationary characteristics of the mountain wind field.

3 Method effect evaluation

3.1 Experimental setting

The experimental object is a large-span steel box girder bridge in a mountain area. The main beam cross-sectional geometric parameters are shown in Table 3. The total length of the bridge is L = 1200; the main span is 800 m; the beam height is 4.5 m; the aspect ratio is 12.5; the mass density is 7850 kg/m³; the first-order torsion frequency is 0.25 Hz; the first-order vertical bending frequency is 0.12 Hz; the structural damping ratio is 0.5%.

The experimental data is a fusion of WRF large eddy simulation and wind tunnel measured data. The WRF simulation uses a triple nested grid, with the innermost grid resolution of 10 m, a time step of 0.1 s, and a simulation time of 72 h, covering typical mountain wind conditions. The terrain elevation data is derived from the 30 m resolution DEM (Digital Elevation Model), and the surface roughness length is randomly perturbed according to the Cauchy distribution. The wind tunnel test model has a scale ratio of 1:200; the passive turbulence generator generates a boundary layer flow field with a turbulence intensity of [5%, 25%]; the wind speed in the test section is [5, 80] m/s. The flow field data acquisition system deploys 64 high-frequency pressure sensors with a sampling frequency of 2000 Hz, covering the circumferential angle of 0°–180° on the main beam surface. The three-dimensional particle image velocimetry system measures the 0.5D-5D flow field area around the main beam, with a spatial resolution of 0.02D × 0.02D × 0.04D and a temporal resolution of 500 Hz. The structural response monitoring uses a laser Doppler vibrometer to measure the torsional and vertical bending modal displacements. The comparison model parameters follow the industry standards. The DES (Detached Eddy Simulation) uses the Spalart-Allmaras turbulence model, and the near-wall grid scale is less than 1; the RANS model uses the Realizable k-ε and SST (Shear Stress Transport) closure scheme. The deep learning benchmark model is compared with PINN (Physics-Informed Neural Network), with a hidden layer dimension of 256, a training cycle of 500 epochs, a batch size of 64, and a cosine annealing strategy. The non-stationary wind field input signal is generated by the random phase reconstruction method; the hardware platform is configured with NVIDIA A100 GPU to accelerate training; the structural solver is based on the Intel Xeon Platinum 8360Y processor for parallel calculation. The CFD simulation uses a 4096-core CPU (Central Processing Unit) cluster. The experimental environment temperature is controlled at 23 ± 1 °C, and the humidity is 50 ± 5% to eliminate the influence of thermal expansion on the modal parameters.

3.2 Correlation between terrain characteristics and wind speed gradient, multi-scale eddy field decomposition

Table 4 quantifies the nonlinear relationship between terrain curvature and wind speed gradient in the generated wind field. For gentle terrain (Slope < 15°), the correlation coefficient is low, indicating that the terrain has a weak disturbance on the flow. In moderate slope terrain (15° ≤ Slope < 30°), the correlation is significantly enhanced, reflecting a stronger interaction between terrain and flow. Steep terrain (Slope ≥ 30°) shows the highest correlation, with an average wind speed gradient of 1.56 ± 0.25 m/(s · m), indicating that the model can accurately capture the flow separation and acceleration effects caused by terrain. These results verify the model's ability to simulate wind field changes in complex mountainous terrain and provide an important quantitative basis for bridge wind resistance design.

Figure 3 shows the vortex field structure based on Fourier multiscale decomposition, including three-dimensional velocity field slices at 8 scales, with a spatial resolution of 64 × 64 × 64 and a scale range of k = 2¹ to k = 2⁸. Each sub-graph corresponds to a vortex field of a specific scale. The color mapping represents the velocity amplitude. The horizontal and vertical coordinates represent the spatial position X and Y, respectively, and the unit is the beam height D. The small-scale vortex field is manifested as a large-scale low-frequency fluctuation, and the velocity amplitude distribution is relatively uniform, reflecting the overall movement trend of the flow field. As the scale increases, the vortex field gradually presents a high-frequency small-scale turbulent structure, and the velocity amplitude distribution is more localized, showing strong spatial inhomogeneity. This multi-scale decomposition can effectively separate the large-scale coherent structure of the boundary layer separation bubble from the small-scale turbulent pulsation in the wake region, verifying the effectiveness of the Fourier integral operator in capturing the law of vortex evolution. It is worth noting that the spatial resolution of the simulation data is limited. For higher scales, the wave number k is close to the Nyquist frequency, and in the Fourier transform, the energy of the high-frequency component is usually low. Affected by the numerical truncation error, there is no obvious difference in the visualization of the higher-scale vortex field. The characteristics of the vortex field at different scales are consistent with theoretical expectations, indicating that this method can precisely analyze the multi-scale characteristics of the flow field and provide a reliable flow field data basis for the wind stability analysis of bridges.

Fig. 3 Multiscale vortex field decomposition.

3.3 Aerodynamic force and structural response time history, lyapunov index and damping ratio

Figure 4 shows the relationship between the aerodynamic force (lift and drag) and structural response (vertical displacement) of the bridge under wind load over time. The lift curve shows a decaying sinusoidal wave characteristic, reflecting the energy dissipation caused by the aerodynamic damping effect; the drag curve contains high-frequency harmonic components, with relatively fixed mean and amplitude, reflecting the unsteady characteristics of the turbulent wind field. The displacement curve phase lags the lift, with an amplitude of 0.2 m, and the frequency is consistent with the lift, both of which are 0.5 Hz, indicating that there is a time delay between the structural response and the aerodynamic force, verifying the aeroelastic coupling effect. The frequency consistency of the aerodynamic force and displacement curves shows that the vibration of the bridge is mainly caused by the periodic excitation of the wind load, while the attenuation of the lift amplitude over time reflects the inhibitory effect of aerodynamic damping on the structural vibration. The data changes intuitively reveal the dynamic characteristics of the wind-bridge coupling system and provide an important basis for wind-resistant stability analysis.

Figure 5 shows the relationship between the Lyapunov index λ_max and the aerodynamic damping ratio ζa with wind speed, which is used to evaluate the stability of the bridge under wind load. The wind speed range is 20–60 m/s. The Lyapunov index curve (orange) transitions from negative to positive as the wind speed increases. The critical wind speed occurs when the Lyapunov index is equal to 0 (≈45 m/s). This critical point has a clear physical significance: it marks the boundary where the bridge system transitions from a stable state to a divergent instability state. When the Lyapunov index is negative, minor disturbances in the system decay over time (stable state); when the index is positive, disturbances are exponentially amplified (divergent state). The critical value of 45 m/s reflects the threshold at which aerodynamic stiffness and structural stiffness reach dynamic equilibrium, at which point the energy input from wind loads exactly offsets the energy dissipated by structural damping. This phenomenon aligns with classical flutter theory, confirming that the system exhibits dynamic instability near the critical wind speed. This result can provide a clear stability criterion for wind-resistant bridge design. The aerodynamic damping ratio curve (cyan) gradually decreases from 0.8% with increasing wind speed. When the wind speed is greater than 60 m/s, the aerodynamic damping ratio is less than 0, indicating that the energy input exceeds the dissipation and triggers flutter. The gray dotted line marks the stability threshold of the Lyapunov index = 0, which clearly distinguishes the stable and unstable areas. This analysis verifies the accuracy of the critical wind speed and provides a theoretical basis for the wind-resistant design of bridges.

Fig. 4 Aerodynamic force and structural response time history.

Fig. 5 Lyapunov index and damping ratio analysis.

3.4 Critical wind speed prediction deviation rate verification

The evaluation process is based on wind tunnel test data, and a multi-model comparison framework is constructed to verify the prediction accuracy. For different working conditions (covering different turbulence intensities [5%, 25%] and incoming flow angles), the measured values of critical wind speed under each working condition are collected, and the prediction error calculation is defined as the relative error${ϵ}_r=\left|V_{\text{critica} \text{pred}}-V_{{\text{critical}}_{ }\text{exp}}\right|/V_{\text{critica} \text{exp}}$. The relative error values of each model are counted for each working condition. To eliminate random interference, the prediction is repeated 10 times for each working condition, and the average is taken. The evaluation uses non-parametric statistical methods, and the error distribution of each model is visualized through box plots. The median, 25%/75% quantiles, and extreme value ranges are quantified, and the error concentration and discrete characteristics are analyzed.

Figure 6 shows the distribution of critical wind speed prediction errors, and Table 5 shows the statistical significance analysis.

Figure 6 compares the relative error distribution of critical wind speed predictions of five models (PINN, DES, k-ε, SST, and Diffusion+FNO). Each box represents the 25%–75% interquartile range of the error. The whiskers extend to 1.5 times the interquartile range, and the external discrete points are outliers. From the data distribution, the Diffusion+FNO method shows the best prediction performance, with a median error of 3.9% and an interquartile range of 2.9%–4.3%, indicating that it has extremely high prediction consistency and robustness in complex mountain wind fields. The errors of traditional CFD models (k-ε and DES) are significantly higher: the median error of k-ε is 17.5%, and the interquartile range is 14.8%–19.1%. There are multiple outliers, mainly because the turbulent closed model cannot accurately capture the transient evolution of the separation bubble. The median error of DES is 12.3%, but the interquartile range is wide (10.1%–14.9%), and the errors of some working conditions are high due to insufficient grid resolution. The PINN and SST models perform in the middle, with a median error of 7.8% for PINN, but some working conditions have high errors (10.4%) due to insufficient physical constraints. The median error of SST is 10.6%, which is limited by the modeling ability of the RANS framework for unsteady flow. The statistical significance analysis in Table 5 shows that the prediction error distribution of the Diffusion+FNO framework is significantly different from that of the traditional method and PINN. The p-values of the comparison between Diffusion+FNO and PINN, DES, k-ε, and SST are all <0.001, which verifies the statistical advantage of its prediction performance. The effect size analysis further quantifies the degree of difference: the Cliff's Delta values are all >0.89 (Large effect), indicating that the error difference has practical engineering significance. The Cohen's d confidence interval does not contain 0, proving the stability of the difference. Diffusion+FNO significantly exceeds the prediction limit of traditional methods through cross-scale modeling and physical embedding, providing a high-confidence solution for the wind-resistant design of bridges in complex wind fields.

Fig. 6 Critical wind speed prediction error distribution.

3.5 Flow field pressure distribution and vortex shedding effect

The flow field pressure distribution on the bridge surface is extracted, and the average pressure coefficient and turbulent pulsation characteristics under wind load are calculated. Combined with experimental data, the changes in the vortex shedding pattern on the bridge surface are analyzed to determine its impact on wind-resistant stability.

Figure 7a shows the average pressure coefficient and the RMS fluctuating pressure within the circumferential angle range of 0° to 180° on the main beam surface. The average pressure coefficient reaches a minimum value of −1.5 on the windward side (θ = 0°), then gradually rises and recovers to 0.25 on the leeward side (θ = 180°). The RMS fluctuating pressure reaches a peak value of 0.7 near θ = 90°, indicating that the turbulent pulsation in this area is most significant. It directly reflects the typical physical phenomenon of flow separation. When the incoming flow separates on the windward side, a shear layer and a recirculation zone will be formed on the leeward side. The position of θ = 90° is in the transition zone of separation bubble reattachment. The flow state in this area is extremely unstable. The periodic generation and fragmentation of the vortex structure in the shear layer leads to strong velocity gradient changes, thereby causing significant local pressure pulsations. This flow characteristic is directly related to the vortex shedding frequency in the process of boundary layer separation-reattachment, and is also an important excitation source for the bridge buffeting response. The data changes show that the flow separation starts at θ≈60° (the average pressure coefficient slope suddenly changes), and the reattachment point is located at θ≈120° (the average pressure coefficient slope slows down). The high fluctuating pressure in this area may cause a sudden change in local aerodynamic loads and aggravate the structural buffeting response.

Figure 7b shows the vibration amplitude (blue bar) corresponding to different vortex shedding frequencies (0.9 Hz to 1.5 Hz), and the red dashed line marks the vertical bending natural frequency of the bridge, 1.2 Hz. The data shows that when the vortex shedding frequency is close to the vertical bending natural frequency of the bridge (1.1 Hz to 1.3 Hz), the vibration amplitude increases significantly, and the peak value appears at 1.2 Hz (A = 0.8). This phenomenon indicates that when the dominant frequency of vortex shedding is coupled with the structural mode, it may trigger limited vibration or even divergent instability, verifying the potential risk of vortex-induced resonance. In actual engineering, the risk of vortex-induced resonance can be avoided through the following measures: optimizing the aerodynamic shape of the main beam, installing guide vanes or nozzles to disrupt periodic vortex shedding patterns; secondly, adjusting structural dynamic characteristics by optimizing mass distribution or stiffness to shift the natural frequency away from common vortex shedding bands; additionally, installing tuned mass dampers and other vibration reduction devices to suppress resonant responses. These methods should be combined with wind tunnel tests for validation to ensure that hazardous frequency locking does not occur within the operating wind speed range.

Fig. 7 Surface pressure characteristics and vortex shedding frequency analysis. (a). Surface pressure characteristics. (b). Vortex shedding frequency analysis.

4 Conclusions

This study proposes a deep learning framework that combines a physical constrained diffusion model with a Fourier neural operator, which successfully solves the problem of predicting the three force coefficients and evaluating the wind stability of bridges in mountainous wind fields. By constructing a conditional diffusion model based on WRF large eddy simulation data, a three-dimensional pulsating wind field with a resolution of 0.1D is generated, and the physical consistency of the generated wind speed field is ensured by embedding the residual constraints of the mass conservation equation and the terrain curvature characteristics. The eight-layer Fourier integral network realizes the decoupling of the multi-scale vortex evolution characteristics in the range of 0.5D∼5D around the main beam, significantly improving the accuracy of cross-scale flow field modeling. Through experimental verification, the median relative error of the Diffusion+FNO framework in the prediction of critical wind speed is 3.9%, which is much lower than the traditional PINN, DES, k-ε and SST models, proving the advantage of frequency domain transformation in flow field modeling. In addition, the FNO-ODE hybrid architecture compresses the time step to 0.001 s by embedding the aerodynamic force prediction into the Newmark-β solution process, ensuring efficient bidirectional coupled iterative calculation and achieving extremely high prediction accuracy. In each step of the research plan, through optimization and experimental verification at different stages, a numerical tool with high accuracy and high efficiency was successfully constructed, which can quantitatively evaluate the wind resistance risk of bridges in complex mountain wind fields, such as aerodynamic instability phenomena such as vortex-induced resonance and flutter. The research results not only improve the accuracy of bridge wind resistance design, but also provide the potential for near-real-time early warning for bridge wind resistance safety monitoring. Future research will extend the model to non-steady rainfall coupling conditions and explore the deployment of edge computing technology to further improve the real-time early warning capability. The practical significance of this study is that by introducing the combination of advanced deep learning technology and physical constraints, a new bridge wind resistance assessment method is provided, which provides a scientific basis for bridge design and wind resistance monitoring in mountain wind fields. In the future, we expect this framework to be widely used in the field of wind engineering, especially for the wind resistance design and monitoring of critical infrastructure, to provide strong support for more intelligent and accurate wind disaster prevention and control.

Acknowledgments

Not applicable.

Funding

This work was supported by Guangxi Natural Science Foundation (2023GXNSFAA026418), National Natural Science Foundation of China (52178468; 52268023) Basic scientific research business expenses in Heilongjiang Province (YWF10236230212).

Conflicts of interest

The authors declare no conflict of interest

Data availability statement

This article has no associated data generated and/or analyzed / Data associated with this article cannot be disclosed due to legal/ethical/other reasons.

Author contribution statement

Fangjin Sun, Jun Peng designed the research study. Yufei Li and Daming Zhang analyzed the data. Fangjin Sun wrote the manuscript. All authors contributed to editorial changes in the manuscript. All authors read and approved the final manuscript.

References

All Tables

All Figures

	Fig. 1 Evaluation framework.
In the text

	Fig. 2 Relationship between terrain height and wind speed field and turbulent kinetic energy spectrum. (a). Relationship between terrain height and wind speed field. (b). Turbulent kinetic energy spectrum.
In the text

	Fig. 3 Multiscale vortex field decomposition.
In the text

	Fig. 4 Aerodynamic force and structural response time history.
In the text

	Fig. 5 Lyapunov index and damping ratio analysis.
In the text

	Fig. 6 Critical wind speed prediction error distribution.
In the text

	Fig. 7 Surface pressure characteristics and vortex shedding frequency analysis. (a). Surface pressure characteristics. (b). Vortex shedding frequency analysis.
In the text