Paper deep dive
Learning Stiffness Dependent Fluid Structure Dynamics from Coarse Flow Representations
Chun-Jun Pu, Li-Wei Chen, Hai-Bo Huang
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 94%
Last extracted: 9/24/2026, 3:58:18 AM
Summary
This paper presents a data-driven framework for long-term prediction of fluid-structure interaction (FSI) dynamics, specifically focusing on the flow-induced vibration of a flexible plate. The proposed method utilizes a stiffness-conditioned neural evolution operator that jointly represents the Eulerian flow field and Lagrangian structural state using a hybrid CNN-Transformer architecture with bidirectional cross-attention. The model captures three distinct stiffness-dependent response regimes (deflected-flapping, deflected, and flapping) and enables accurate, differentiable aerodynamic force reconstruction via a Derivative-Moment Transformation (DMT) module, facilitating efficient parameter studies and optimization for flow-energy-harvesting applications.
Entities (8)
Relation Signals (5)
Derivative-Moment Transformation → enables → Aerodynamic Force Reconstruction
confidence 96% · we develop a differentiable aerodynamic-force module based on the derivative-moment transformation (DMT)... enabling accurate reconstruction of lift and drag.
Neural Evolution Operator → represents → Fluid-Structure Interaction
confidence 95% · A stiffness-conditioned neural evolution operator jointly represents the Eulerian flow field and Lagrangian structural state.
Neural Evolution Operator → captures → Flow-Induced Vibration
confidence 94% · a single operator captures three stiffness-dependent response regimes: deflected--flapping, deflected, and flapping.
Bending Stiffness → conditions → Neural Evolution Operator
confidence 93% · with nondimensional bending stiffness as a global conditioning variable.
CNN-Transformer → implements → Bidirectional Cross-Attention
confidence 90% · Bidirectional cross-attention couples fluid and structural representations within a hybrid CNN-Transformer architecture.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:This paper develops a data-driven framework for long-term prediction of fluid--structure interaction (FSI) dynamics, focusing on the flow-induced vibration (FIV) of a flexible plate. A stiffness-conditioned neural evolution operator jointly represents the Eulerian flow field and Lagrangian structural state. The plate is represented by 101 ordered structural tokens carrying nodal coordinates and velocities, with nondimensional bending stiffness as a global conditioning variable. Bidirectional cross-attention couples fluid and structural representations within a hybrid CNN-Transformer architecture. Trained with staged multi-step autoregressive rollouts and symmetry-reflected trajectories, a single operator captures three stiffness-dependent response regimes: deflected--flapping, deflected, and flapping. The predicted trajectories preserve the principal flow structures, structural oscillations, and dominant frequencies, while blind 1000-step rollouts remain bounded. The operator also interpolates to stiffness values excluded from training. To reduce sensitivity to under-resolved near-wall gradients in force reconstruction, we develop a differentiable aerodynamic-force module based on the derivative-moment transformation (DMT). Conventional wall-stress surface integrals are replaced by an enclosed 2D curve integral around the core vortex region, enabling accurate reconstruction of lift and drag. A Signed Distance Function (SDF) and smoothed Dirac-delta formulation make the integration fully differentiable while preserving gradient flow. The proposed framework provides an accurate and differentiable surrogate for stiffness-dependent FSI dynamics, enabling efficient parameter studies and future stiffness optimization for flow-energy-harvesting applications.
Tags
Links
- Source: https://arxiv.org/abs/2609.26816v1
- Canonical: https://arxiv.org/abs/2609.26816v1
Trouble viewing inline? Open PDF directly →
Full Text
97,636 characters extracted from source content.
Expand or collapse full text
Learning Stiffness-Dependent Fluid–Structure Dynamics from Coarse Flow Representations Chun-Jun Pu a , Li-Wei Chen b,∗ , Hai-Bo Huang a,∗ a University of Science and Technology of China, Hefei 230026, China b Beijing Institute of Aeronautical Systems Engineering, Beijing 100076, China Abstract This paper develops a data-driven framework for long-term prediction of fluid– structure interaction (FSI) dynamics, focusing on the flow-induced vibration (FIV) of a flexible plate. A stiffness-conditioned neural evolution operator is developed by jointly representing the Eulerian flow field and the Lagrangian structural state. The flexible plate is described by 101 ordered structural to- kens carrying nodal coordinates and velocities, while the nondimensional bend- ing stiffness is introduced as a global conditioning variable. Bidirectional cross- attention couples the fluid and structural representations within a hybrid CNN– Transformer architecture. Trained with staged multi-step autoregressive rollouts and symmetry-reflected trajectories, a single operator captures three distinct stiffness-dependent response regimes, including deflected–flapping, deflected, and flapping dynamics. The predicted trajectories preserve the principal flow organization, structural oscillation characteristics, and dominant frequencies, while blind 1000-step rollouts remain bounded. The same operator also inter- polates to stiffness values excluded from training. To reduce the sensitivity of force reconstruction to under-resolved near-wall gradients, we construct a differentiable aerodynamic-force reconstruction mod- ule based on the derivative-moment transformation (DMT). By replacing con- ventional wall-stress surface integrals with an enclosed 2D curve integral sur- ∗ Corresponding author Email addresses: liwei.chen.aero@gmail.com (Li-Wei Chen), huanghb@ustc.edu.cn (Hai-Bo Huang) arXiv:2609.26816v1 [physics.flu-dyn] 19 Sep 2026 rounding the core vortex region, the aerodynamic lift and drag forces can be accurately reconstructed. This differentiable integration is implemented through a Signed Distance Function (SDF) and a smoothed Dirac-delta formulation, pre- serving gradient flow throughout the force-reconstruction procedure. Overall, the proposed framework provides an accurate and differentiable surrogate for stiffness-dependent FSI dynamics, offering a foundation for efficient parameter studies and future stiffness optimization in flow-energy-harvesting applications. Keywords: Fluid–structure interaction; Flexible plates; Fluid-induced vibration; Neural evolution operator; Aerodynamic force reconstruction; Derivative-moment transformation 1. Introduction Fluid–structure interaction (FSI) involving slender and flexible bodies arises in a wide range of engineering and biological systems, including aeroelastic structures, energy-harvesting devices, aquatic propulsion, and flow-responsive materials. These problems are governed by a nonlinear feedback loop: the fluid loading deforms the structure, while the resulting structural motion continu- ously modifies the fluid boundary and reorganizes the surrounding vortical flow. A representative configuration is the inverted flexible plate, whose downstream end is clamped while its upstream end is free. Depending on its nondimensional bending stiffness and other fluid–structural parameters, the plate can exhibit straight, deflected, periodic flapping, and more complex oscillatory states [1–3]. The transitions between these regimes are accompanied by substantial changes in the mean plate configuration, oscillation amplitude, dominant frequency, wake topology, and aerodynamic loading. These stiffness-dependent changes are also directly relevant to flow-energy-harvesting systems, whose performance depends strongly on the amplitude, frequency, and loading characteristics of the structural motion. Resolving such stiffness-dependent dynamics therefore requires a model that simultaneously represents the Eulerian flow field, the moving Lagrangian structure, and their strongly coupled temporal evolution. 2 High-fidelity numerical methods can reproduce these phenomena but remain computationally demanding when long trajectories or densely sampled param- eter spaces are required. This cost becomes particularly restrictive in repeated- query applications such as parameter inference, uncertainty quantification, flow control, and aeroelastic design. Data-driven surrogate models provide a pos- sible alternative by learning the evolution or solution operator directly from high-fidelity simulations. Early operator-learning frameworks, represented by DeepONet and the Fourier neural operator (FNO), demonstrated that map- pings between infinite-dimensional function spaces can be approximated with- out solving each new PDE instance from the beginning [4, 5]. Graph-based simulators further extended learned physical modeling to irregular meshes and interacting discretized objects [6]. More recently, attention-based neural oper- ators, including OFormer and GNOT, have provided flexible representations of nonlocal spatial interactions, irregular sampling, multiple input functions, and multiscale PDE solutions [7, 8]. Applying these ideas to FSI is more difficult than predicting a fluid field over a fixed domain. The structural geometry is itself part of the evolving phys- ical state and determines where momentum is exchanged between the fluid and the solid. A fluid-only Eulerian surrogate must infer this moving interface im- plicitly from the velocity and pressure fields, so small interface errors can be reintroduced at every autoregressive step. Existing studies have explored CNN- based moving-boundary prediction, graph representations of coupled systems, and FNO-based hybrid FSI solvers [9–11]. Differentiable hybrid approaches have also embedded coarse numerical FSI models within trainable recurrent architec- tures, allowing gradients to propagate through the coupled rollout trajectory [12]. In parallel, codomain-attention operators have shown that tokenizing and exchanging information among different physical variables can improve the rep- resentation of multiphysics systems [13]. Nevertheless, explicit learned coupling between a dense Eulerian flow representation and an ordered Lagrangian flexible structure remains comparatively underexplored. In particular, directly concate- nating a geometry mask with the fluid channels does not preserve the distinct 3 physical identities and discretizations of the fluid and structural states. A second challenge concerns long-horizon prediction. Neural PDE models are often optimized using one-step or short-window errors, whereas practical FSI applications require repeated reuse of predicted states. Even a small one- step bias can accumulate during autoregressive rollout, leading to phase drift, blurred vortical structures, duplicated moving interfaces, or collapse toward an unphysical steady state. Multi-step message-passing solvers, refinement-based models, and denoising operator Transformers have been proposed to improve the stability of long PDE trajectories [14–16]. Recent systematic comparisons of one-step, differentiably unrolled, and gradient-truncated unrolled training further show that temporal unrolling improves autoregressive inference by ex- posing the model to its own rollout distribution. Full temporal differentiabil- ity provides additional benefits within matched training configurations, while non-differentiable unrolling can still retain much of the improvement associated with reducing the rollout distribution shift [17]. Benchmark studies have fur- ther shown that autoregressive performance depends strongly on the training horizon, noise treatment, temporal resolution, and underlying PDE dynamics [18]. For FSI, this difficulty is compounded by the need to preserve not only pointwise fields but also the coupled oscillatory attractor, including the mean deformation, dominant frequency, phase-space structure, and wake organization. Moreover, most existing surrogate models are demonstrated for a single physical condition or a limited response regime. Whether one conditional evolution op- erator can reproduce multiple stiffness-dependent FSI regimes while remaining stable under a fully autoregressive rollout is still insufficiently established. A third challenge is the recovery of aerodynamic forces from neural-predicted flow fields. Lift and drag are global functionals of the coupled solution and pro- vide a more integrated physical assessment than pointwise field errors. Conven- tional force evaluation uses pressure and viscous tractions on the instantaneous solid boundary. On a coarse neural grid, however, the thin boundary layer and near-wall gradients are generally under-resolved, making direct surface-stress 4 integration sensitive to local smoothing and differentiation errors. Control-surface formulations offer an alternative by transferring force evalua- tion from the body surface to a better-resolved outer-flow region. The derivative- moment transformation (DMT) expresses the unsteady force through vorticity moments and flux terms defined over a control domain and its enclosing contour [19, 20]. Control-surface force measurements and more recent weighted-integral developments have demonstrated the accuracy and diagnostic value of this class of formulations [21, 22]. For integration with neural models, force accuracy alone is not sufficient: the force-reconstruction operator should also retain derivatives with respect to the predicted velocity and pressure fields. Differentiable surrogate-based design has demonstrated the value of combining level-set geometry descriptions, auto- matic differentiation, and neural flow predictions [23]. Fully differentiable and GPU-accelerated CFD platforms have further enabled gradient-based inverse modeling and optimization through complete turbulence and FSI simulations [24]. Such approaches, however, generally require differentiating through the numerical solver itself. A complementary capability is to construct a differen- tiable force readout directly on coarse flow fields produced by a learned evolution model. This requires a continuous grid-compatible representation of the control contour and a regularized integration rule that preserves both the DMT force identity and a usable gradient path. In this work, we develop a stiffness-conditioned neural evolution framework for long-horizon prediction and force-based assessment of a two-dimensional in- verted flexible plate at Re = 200. The model evolves the fluid velocity and pressure on a fixed Eulerian grid while independently tracking the nodal coor- dinates and velocities of 101 ordered Lagrangian plate points. A convolutional encoder extracts spatial features from the flow field, which are then organized into a compact sequence of fluid tokens for Transformer attention. The position and velocity of each plate point are jointly embedded as a structural token, while the nondimensional bending stiffness is introduced through an additional global conditioning token. Bidirectional cross-attention transfers information from the 5 fluid to the structure and from the updated structure back to the fluid. The training data comprise 36 original stiffness-dependent trajectories, sup- plemented by eight reflected trajectories. To improve recursive prediction, the model is first trained using four-step rollouts and is subsequently fine-tuned using eight-step rollouts. Starting from a single coupled state, the complete fluid–structure state is advanced recursively without intermediate ground-truth correction. The resulting model provides a unified description of three stiffness- dependent response regimes, with particular attention to the deflected and flap- ping dynamics. The principal contributions of this work are summarized as follows. (1) Coupled Eulerian–Lagrangian representation for flexible-plate FSI. A neural evolution operator is developed to jointly represent the Eulerian flow field and the Lagrangian flexible structure, with bidirectional cross-attention coupling the two physical states. The structural state includes both nodal coor- dinates and velocities, while the nondimensional bending stiffness is introduced as a global conditioning variable. (2) Unified autoregressive prediction across multiple stiffness-dependent FSI regimes. A single stiffness-conditioned operator is trained over 36 bending stiffnesses and captures the deflected–flapping, deflected, and flapping regimes in a fully autoregressive manner. The predicted dynamics remain stable over long rollouts and interpolate to unseen stiffness conditions. (3) Differentiable aerodynamic-force reconstruction from coarse neu- ral flow fields. A differentiable force-reconstruction operator is developed by combining the DMT with a fixed control contour represented by a signed distance function and a smoothed Dirac-delta formulation. The resulting for- mulation recovers aerodynamic forces from coarse predicted flow fields while retaining gradients with respect to the flow variables. The remainder of this paper is organized as follows. Section 2 describes the physical FSI system, numerical solver, multi-stiffness dataset, data augmenta- tion, and autoregressive sequence construction. Section 3 presents the coupled CNN–Transformer evolution operator, the staged rollout-training strategy, and 6 the differentiable DMT force-reconstruction method. Section 4 evaluates the prediction accuracy across the stiffness parameter space, representative multi- regime dynamics, aerodynamic-force reconstruction, long-horizon rollout behav- ior, and interpolation to unseen stiffnesses. The main conclusions and limita- tions are summarized in Section 5. 2. Problem Description and FSI System We consider a two-dimensional flexible plate immersed in an incompressible viscous fluid flow. The structural boundary of the flexible plate is discretized into a set of N s = 101 Lagrangian points. Each Lagrangian point i∈1,...,N s is characterized by its instantaneous Cartesian coordinates x i (t) = [x i (t),y i (t)] T and structural velocity [u s,i (t),v s,i (t)] T . A single nondimensional bending stiff- ness κ is assigned uniformly to the plate in each simulation. The fluid state is defined on an Eulerian grid covering the computational domain Ωand con- sists of the velocity u = [u,v] T and pressurep, while the vorticity ω is derived from the velocity field. The coupling at the fluid-solid interface dictates that the fluid exerts aerodynamic forces on the plate, causing large-amplitude non- linear oscillations, which in turn modulate the vortex shedding patterns in the wake, exhibiting multiple distinct limit-cycle oscillation (LCO) modes as the nondimensional bending stiffness varies. The physical configuration and the computational domains are illustrated in Fig. 1. The flexible plate is aligned with the incoming uniform flow in its undeformed configuration, with its downstream end clamped and its upstream end free. The full non-uniform CFD domain is used to generate the reference trajectories, whereas a near-body window containing the plate and the principal wake structures is extracted for neural prediction. 2.1. Physical configuration and governing equations The present study focuses on the flow-induced vibration of an inverted flex- ible plate. The plate has a nondimensional length L = 1 and is initially aligned 7 -15.07025.15 x/L -15.07 0 15.07 y / L Uniform flow, U Near-body window Uniform far field: u = U, v = 0 (a) Full non-uniform CFD domain -1.0001.001.55 x/L -1.28 0 1.27 y / L Free end Clamped end 101 Lagrangian points Representative plate positions U Re = 200, M = 2.9, K S = 1000,varied 256 × 256 Eulerian samples, x/L =y/L = 0.01 L (b) Near-body computational window Figure 1: Physical configuration and computational domains of the flexible-plate FSI problem. (a) Full non-uniform CFD domain and the location of the near-body prediction window. (b) Near-body Eulerian window used by the neural evolution operator. with the incoming flow. Its downstream end is clamped at (x,y) = (0, 0), whereas its upstream end is free. Thus, the undeformed plate occupies −1 ≤ x≤ 0 and y = 0. A uniform flow with velocity u ∞ = (U ∞ , 0) is imposed in the positive x direction. This inverted configuration permits large-amplitude defor- mation of the upstream free end and produces a strongly coupled interaction between the plate motion and the wake dynamics. The fluid state is described by q f (x,t) = [u(x,t),v(x,t),p(x,t)] T ,(1) where u and v are the Cartesian velocity components and p denotes the pressure channel retained in the learning problem. Its normalization relative to the pressure appearing in the governing equations is specified below. ω = ∂v ∂x − ∂u ∂y .(2) The vorticity is not evolved as an independent variable by the neural model; it is derived from the velocity field and is used to characterize the wake regimes. Using the plate length L, reference velocity U ∞ , fluid density ρ f , and convec- tive time scale L/U ∞ , we define the pressure appearing in the governing equa- tions as p ⋆ = p dim /(ρ f U 2 ∞ ). The dimensionless incompressible Navier–Stokes 8 equations are ∇·u = 0,(3) ∂u ∂t +u·∇u =−∇p ⋆ + 1 Re ∇ 2 u,(4) where Re = U ∞ L ν = 200.(5) and ν is the kinematic viscosity. The dimensionless fluid stress tensor is written as σ f =−p ⋆ I + 1 Re ∇u +∇u T .(6) For the learning problem and the subsequent aerodynamic-force reconstruc- tion, the retained pressure channel p is normalized by the dynamic pressure q ∞ = 1 2 ρ f U 2 ∞ . Consequently, p = p dim /q ∞ = 2p ⋆ . This distinction is purely a normalization convention and does not change the underlying dimensional pressure field. The plate centreline is discretized by N s = 101 Lagrangian nodes connected by 100 geometrically nonlinear beam elements. The structural state associated with node i is defined as s i (t) = [x i (t), y i (t), u s,i (t), v s,i (t)] T , X s (t) =s i (t) N s i=1 ∈R N s ×4 . (7) Although the structural solver retains the complete finite-element degrees of freedom, only the in-plane nodal coordinates and velocities are retained in the learning problem. The principal control parameter is the nondimensional bending stiffness κ≡ K B = B ρ f U 2 ∞ L 3 ,(8) where B denotes the dimensional bending rigidity. The nondimensional axial or extensional stiffness is held fixed at K S = 1000.(9) Although the formulation permits a spatially varying local stiffness k i , all Lagrangian points in each simulation are assigned the same value, such that 9 k i = κ for i = 1,...,N s . The nondimensional mass ratio and Poisson ratio are fixed at m ∗ = ρ s h ρ f L = 2.9, ν s = 0.25,(10) respectively, while the structural damping coefficients are set to zero. The fluid traction is transferred to the Lagrangian structure, and the no-slip condi- tion requires the interpolated fluid velocity to match the local structural velocity on the moving plate. Consequently, the plate deformation modulates the vor- tex formation and shedding process, while the resulting fluid loading drives the structural motion. 2.2. Coupled numerical solver The high-fidelity trajectories are generated using a partitioned FSI solver. The incompressible flow is advanced by a D2Q9 multiple-relaxation-time lattice Boltzmann method on a nonuniform Cartesian grid. The full fluid grid contains 1468× 985 points and covers −15.073≤ x≤ 25.150, −15.073≤ y ≤ 15.073.(11) The mesh is refined in the vicinity of the plate and gradually stretched towards the far field. A uniform far-field state is imposed at the outer boundaries. The fluid and structure are coupled through an immersed-boundary forcing procedure. At every coupled time step, the fluid variables are first advanced, the immersed-boundary force is evaluated and transferred to the plate, and the structural state is subsequently updated. The plate dynamics are solved using a geometrically nonlinear finite-element formulation with a Newmark time-integration scheme. The fluid reference values used by the source solver are L = 1, U ∞ = 0.025, and ρ f = 1, yielding a convective reference time T ref = L/U ∞ = 40. The internal solver time step is ∆t = 0.01, corresponding to the dimensionless interval ∆t T ref = 2.5× 10 −4 .(12) 10 The flow and structural states used for data construction are written every 0.1 dimensionless time units. Before storage, the nodal coordinates are scaled by L, both the fluid and structural velocity components are scaled by U ∞ , and the retained pressure channel p is scaled by the dynamic pressure q ∞ = 1 2 ρ f U 2 ∞ . Thus, the stored pressure satisfies p = 2p ⋆ . Dimensional quantities are denoted explicitly only when needed hereafter. 2.3. Multi-stiffness FSI dataset Independent simulations are conducted by varying only the nondimensional bending-stiffness parameter κ, while all other fluid and structural parameters remain unchanged. The reference dataset consists of 36 independently simulated flexible-plate cases, with κ uniformly sampled from 0.05 to 0.40 at an interval of 0.01. The set of independently simulated stiffnesses is therefore defined as K train =0.05 + 0.01j|j = 0, 1,..., 35.(13) Collectively, these 36 stiffness conditions represent the characteristic plate vi- bration modes and coupled response regimes considered in this study. The reflection-based augmentation used to construct the final training set is de- scribed in Section 2.4. 2.3.1. Stiffness-dependent response regimes Previous studies have shown that the dynamical response of an isolated inverted flexible plate can be classified into several characteristic regimes, in- cluding the straight, flapping, deflected, and deflected–flapping modes. In the present study, the selected bending-stiffness range focuses on three dynamically active regimes: the deflected–flapping, deflected, and flapping modes. The structural response is characterized by the peak-to-peak transverse dis- placement A of the free tip, reported in normalized form as A L = max t∈T y tip (t)− min t∈T y tip (t) L ,(14) where L is the undeformed plate length and T denotes the statistically station- ary sampling interval. 11 0.050.100.150.200.250.300.350.40 Bending stiffness, κ 0.0 0.4 0.8 1.2 1.6 Flapping amplitude, A / L Deflected-flapping Deflected Flapping (a) Deflected-flapping κ=0.05 (b) Deflected κ=0.10 (c) Flapping κ=0.25 (d) Figure 2: Stiffness-dependent response regimes of the isolated inverted flexible plate. (a) Nor- malized peak-to-peak transverse displacement of the free tip. (b–d) Representative full-body envelopes of the deflected–flapping, deflected, and flapping modes, respectively. The dashed lines visually separate the response regimes represented by the discrete stiffness samples. As shown in figure 2, the densely sampled stiffness cases reveal three dis- tinct structural-response regimes. At κ = 0.05, the plate exhibits the deflected– flapping mode, characterized by a strongly deflected mean configuration accom- panied by finite-amplitude oscillations. For 0.06 ≤ κ ≤ 0.12, the response is classified as the deflected mode, in which the plate oscillates weakly around a bent mean configuration. When 0.13≤ κ≤ 0.40, the plate undergoes sustained large-amplitude flapping across the undeformed centerline. These three regimes collectively describe the characteristic stiffness-dependent structural responses covered by the present dataset. For each value of κ, a statistically developed segment over 250.1≤ t≤ 280.0(15) 12 is extracted from the original FSI simulation. The raw output interval is ∆t raw = 0.1, giving 300 synchronized fluid–structure snapshots for each con- dition. This late-time interval is selected to focus on the developed stiffness- dependent oscillatory responses rather than the initial transient from the un- deformed state. For every stiffness condition, the selected trajectory contains at least two complete oscillation cycles, thereby retaining the essential periodic evolution of both the plate motion and the associated wake dynamics. For neural prediction, a local region surrounding the moving plate and its near wake is extracted from the full numerical domain. The learning domain is Ω ML = [−1.00, 1.55]× [−1.28, 1.27](16) and is represented on a uniform 256× 256 Cartesian grid. Each stored CFD snapshot contains the velocity components u and v, the pressure p, and a binary geometry mask. Since the moving plate is represented independently by its La- grangian state, the geometry mask is not supplied to or predicted by the neural evolution operator. The fluid state used in the learning problem is therefore Q n f = [u n , v n , p n ]∈R 3×256×256 .(17) The corresponding structural state contains the nodal positions and velocities, X n s = x n i , y n i , u n s,i , v n s,i 101 i=1 ∈R 101×4 .(18) To reduce temporal redundancy, one snapshot is retained from every five raw outputs. The effective time interval of the learning problem is therefore ∆t ML = 5∆t raw = 0.5.(19) Each stiffness condition consequently contains 60 frames at t = 250.1, 250.6,..., 279.6.(20) The 36 original stiffness conditions contain 2160 synchronized fluid–structure states. Eight additional reflected trajectories, each containing 60 frames, con- tribute another 480 states. The augmented training set therefore contains 44 trajectories and 2640 synchronized fluid–structure states in total. 13 Table 1: Summary of the physical system and multi-stiffness dataset. QuantityValue Reynolds numberRe = 200 Plate lengthL = 1 Mass ratiom ∗ = 2.9 Poisson ratioν s = 0.25 Fixed extensional stiffnessK S = 1000 Training bending stiffnesses0.05–0.40, ∆κ = 0.01 (36 cases) Structural discretization101 nodes, 100 beam elements Full fluid grid1468× 985 Learning grid256× 256 Learning-domain extent[−1.00, 1.55]× [−1.28, 1.27] Raw snapshots per stiffness300 Raw output interval∆t raw = 0.1 Temporal subsampling factor5 Learning interval∆t ML = 0.5 Subsampled frames per stiffness60 Predicted fluid variablesu,v,p Predicted structural variablesx i ,y i ,u s,i ,v s,i , i = 1,..., 101 2.4. Preprocessing, reflection augmentation, and autoregressive sequence con- struction To reduce the directional imbalance associated with the deflected responses, reflection-based data augmentation is applied to the trajectories at K mirror =0.05, 0.06,..., 0.12. For a reflection about the plate centreline, the fluid variables are transformed according to u mir (x,y,t) = u(x,−y,t), v mir (x,y,t) =−v(x,−y,t), p mir (x,y,t) = p(x,−y,t), (21) 14 while the structural state is transformed as x mir i = x i , y mir i =−y i , u mir s,i = u s,i , v mir s,i =−v s,i .(22) The stiffness value is unchanged under reflection. This procedure introduces eight additional trajectories without adding new stiffness values, resulting in 36 original trajectories and eight reflected trajectories for training. The velocity and pressure channels have different numerical scales. Global minimum and maximum values are evaluated separately for q c ∈ u,v,p over all 44 original and reflected training trajectories. The fluid channels are then normalized according to eq c = q c − q c,min q c,max − q c,min + 10 −8 .(23) The structural coordinates and velocities are retained in their nondimen- sional physical scales without an additional dataset-wise min–max normaliza- tion. The stiffness condition is mapped to [−1, 1] using eκ = 2 κ− κ min κ max − κ min − 1, κ min = 0.05, κ max = 0.40.(24) Training samples are constructed using multi-step autoregressive rollout win- dows. Starting from an initial coupled state (Q n f ,X n s ), the model recursively predicts b Q n+j f , c X n+j s , j = 1,...,N r ,(25) where N r denotes the rollout horizon adopted for a given training config- uration, and the prediction at step j is used as the input to step j + 1. This construction allows the rollout horizon to be refined without changing the defi- nition of the underlying fluid–structure sequence. During training, zero-mean Gaussian perturbations with a standard devia- tion of 0.005 are added at each rollout step to the normalized fluid variables and all four components of the structural state. This exposes the network to off-trajectory states and improves its robustness during free autoregressive roll- out. 15 The model adopted in this study is trained using a staged rollout strategy. It is first optimized with a four-step autoregressive horizon (N r = 4). The resulting R4 checkpoint is then used to initialize an eight-step autoregressive fine-tuning stage (N r = 8). The selected R4-to-R8 checkpoint is used for all neural predictions and aerodynamic-force reconstructions reported in this paper. For each 60-frame stiffness trajectory, the R4 and R8 stages contain 56 and 52 admissible rollout windows, respectively, corresponding to 2464 and 2288 training windows over the 44 original and reflected training trajectories. The dataset spans stiffness-dependent changes in plate deformation, oscilla- tion amplitude, characteristic frequency, and wake topology. These qualitative response regimes are first identified from the plate envelopes and phase-aligned vorticity fields and are subsequently quantified using the free-tip trajectory, dominant frequency, mean drag, and root-mean-square lift. This organization provides the physical basis for evaluating whether a single conditional neural evolution operator can reproduce multiple FSI regimes and remain stable dur- ing long-term autoregressive prediction. 3. Methodology The core of the proposed framework consists of two main pillars: a hybrid CNN-Transformer neural evolution operator that models the autoregressive spa- tiotemporal FSI coupling, and a differentiable aerodynamic force reconstruction module based on the DMT force formulation, facilitated by an SDF-based soft integration method. An overview of the complete framework is shown in Fig. 3. 3.1. CNN-Transformer Hybrid Evolution Operator To resolve the multi-scale characteristics of the flow field and the discrete nature of the structure, we map the Eulerian fluid field and the Lagrangian solid points into a unified latent space using a hybrid tokenization strategy. At each discrete time step n, the normalized fluid state is represented by e Q n f = [eu n ,ev n ,ep n ]∈R 3×H×W , H = W = 256.(26) 16 Flatten + position prepend global stiffness token H c s θ :( ̃ Q n f ,X n s , ̃ κ)↦( ̂ Q n+1 f , ̂ X n+1 s ) ⋯ Predictions are recursively reused as inputs; all steps share θ and ̃ κ. Training: pure multi-step supervision against paired CFD and plate states ∂(C L ,C D )/∂( ̂ u, ̂ v, ̂ p) available by automatic differentiation Post-processing in this work; not included in the training loss. gradient path CNN encoder 4× downsample Fluid Transformer 4 self-attention layers Point embedding MLP + position ̃ κ MLP κ Fluid → Solid cross-attention Q s ; K f ,V f (2 layers) Solid → Fluid cross-attention Q f ; K s ,V s (2 layers) CNN decoder 4× upsample Solid decoder 128→64→4 ̂ Q n+1 f [ ̂ u, ̂ v, ̂ p]∈ℝ 3×256×256 ̂ X n+1 s [ ̂ x, ̂ y, ̂ u s , ̂ v s ]∈ℝ 101×4 q n =( ̃ Q n f ,X n s ) θ ̂ q n+1 θ ̂ q n+2 Predicted flow field [ ̂ u, ̂ v, ̂ p] Fixed contour φ, δ ε (φ) Differentiable DMT ∂ t ,∇,ω; H ε ,δ ε Eulerian area integration Lift / drag C L (t), C D (t) (a) Coupled neural evolution operator ̃ Q n f = [ ̃ u, ̃ v, ̃ p] X n s = [x,y,u s ,v s ] ∈ ℝ 101× 4 Bidirectional FSI coupling (b) Pure autoregressive multi-step rollout(c) Differentiable DMT force reconstruction 256×128 101×128102×128 Figure 3: Overview of the proposed framework. (a) The Eulerian flow field and the Lagrangian plate states, including nodal positions and velocities, are encoded as fluid and structural to- kens, respectively. The normalized bending stiffness is introduced as a global conditioning token, and bidirectional cross-attention exchanges information between the two physical rep- resentations. (b) The coupled state predicted at each time step is recursively reused as the input for subsequent autoregressive evolution. (c) Aerodynamic lift and drag are reconstructed from the predicted flow field using a differentiable DMT operator defined on a fixed level-set control contour. The instantaneous state of the flexible plate is represented by N s = 101 La- grangian points, X n s = s n 1 ,s n 2 ,...,s n N s T ∈R N s ×4 , s n i = x n i , y n i , u n s,i , v n s,i T .(27) The dimensionless bending stiffness is supplied to the network as the normalized scalar eκ. The one-step coupled evolution learned by the model is expressed as b e Q n+1 f , c X n+1 s =G θ e Q n f ,X n s ,eκ ,(28) where G θ denotes the neural evolution operator parameterized by θ, and the hat denotes a predicted quantity. Repeated application of equation (28) produces an autoregressive rollout of the coupled fluid–structure dynamics. 17 3.1.1. Fluid Tokenization via CNN At each time step n, the normalized two-dimensional fluid state e Q n f ∈ R 3×H×W , consisting of the velocity components eu and ev and the pressure ep, is processed by a convolutional encoder. The purpose of the encoder is to ex- tract hierarchical spatial features and convert the dense Eulerian flow field into a compact sequence of fluid tokens suitable for global attention. The convolutional encoder consists of four successive downsampling blocks. Each block employs a 4× 4 convolution with stride 2 and padding 1, followed by batch normalization and a GELU activation. The spatial resolution and channel dimension evolve as 3× 256× 256→ 32× 128× 128→ 64× 64× 64 → 128× 32× 32→ 128× 16× 16. (29) The resulting feature map can be expressed as F n f =E CNN e Q n f ∈R d m ×16×16 , d m = 128,(30) where E CNN denotes the convolutional encoder and d m is the latent embedding dimension. The feature map is subsequently flattened along its two spatial dimensions to produce N f = 16× 16 = 256 fluid tokens. A learnable positional embedding P f is added to retain the spatial identity of each token: H n f,0 = Flatten F n f +P f ∈R N f ×d m , N f = 256.(31) Before interacting with the structural representation, the fluid-token sequence is processed by a four-layer Transformer encoder. Each layer contains four-head self-attention, a pre-normalization structure, and a feed-forward network with hidden dimension 4d m = 512. The encoded fluid representation is therefore given by H n f =T f H n f,0 ∈R N f ×d m ,(32) where T f denotes the fluid Transformer encoder. 18 This hybrid representation is particularly suitable for the present FSI prob- lem. The convolutional encoder extracts local spatial features associated with velocity gradients and vortical structures while reducing the original 256× 256 Eulerian grid to only 256 latent tokens. Applying self-attention directly to all 256 2 grid points would be computationally prohibitive because its cost increases quadratically with the sequence length. The CNN-based tokenization reduces the attention sequence length by a factor of 256, making global interaction mod- eling computationally tractable. The subsequent Transformer encoder complements the local receptive fields of the convolutional layers by capturing nonlocal dependencies between the near-plate flow and the downstream wake. The resulting fluid representation therefore combines local flow features, long-range wake interactions, and com- putationally efficient global attention before coupling with the flexible-plate rep- resentation. 3.1.2. Solid Tokenization and Physical Embedding The structural state is represented by N s = 101 ordered Lagrangian points. Since this representation is already discrete, each point is directly embedded as a structural token without further spatial downsampling. For the ith point, the instantaneous nodal state s n i = x n i , y n i , u n s,i , v n s,i T is mapped into the latent space by an embedding MLP and combined with a learnable positional embedding: h n,0 s,i = MLP s (s n i ) +p s,i , i = 1,...,N s .(33) Here, p s,i distinguishes the fixed ordering of the Lagrangian points, while s n i describes both the instantaneous geometry and kinematics of the plate. The resulting structural-token sequence is H n s,0 = h h n,0 s,1 ,...,h n,0 s,N s i T ∈R N s ×d m , d m = 128.(34) 19 The normalized bending stiffness eκ is represented by an additional token, h 0 κ = MLP κ (eκ) +p κ ,(35) which is prepended to the structural-token sequence: H n sκ,0 = h 0 κ ;H n s,0 ∈R (N s +1)×d m =R 102×128 .(36) The stiffness token provides a global conditioning variable for the structural representation and participates together with the Lagrangian tokens in the sub- sequent cross-attention operations. 3.1.3. Bidirectional Cross-Attention Coupling Mechanism The fundamental FSI interaction is modeled via a bidirectional cross-attention mechanism, alternating between Fluid-to-Solid (f → s) and Solid-to-Fluid (s→ f) blocks. For the f → s coupling: b H n s = softmax Q s K T f √ d k ! V f + H n s .(37) where Q s comes from the solid tokens, while K f and V f come from the fluid tokens. b H n f = softmax Q f K T s √ d k V s + H n f .(38) Both directional coupling modules are implemented using two-layer Trans- former decoders. Each decoder layer consists of four-head attention, pre-normalization, residual connections, and a feed-forward network with hidden dimension 4d m = 512. For d m = 128, the feature dimension of each attention head is d k = 32. The fluid-to-solid attention operates between N s +1 = 102 structural tokens, including the stiffness token, and N f = 256 fluid tokens. Its attention matrix therefore has the dimension A f→s ∈R (N s +1)×N f =R 102×256 .(39) The subsequent solid-to-fluid attention uses the updated structural representa- tion as its memory, giving A s→f ∈R N f ×(N s +1) =R 256×102 .(40) 20 The two operations are evaluated sequentially, such that the structural tokens are first conditioned on the fluid representation and the fluid tokens are then up- dated using the coupled structural representation. This formulation couples the Eulerian and Lagrangian representations in the latent space without introducing an explicit pointwise interpolation operator. 3.1.4. Autoregressive Rollout and Pure Supervision The updated latent states are passed through the decoders to predict the next time step. The available trajectory for each stiffness contains T seq = 60 frames. The network is trained autoregressively using pure multi-step supervi- sion. Following the sequence construction described in Section 2.4, a staged au- toregressive training strategy is adopted. The network is first trained using a four-step rollout horizon (N r = 4). The resulting R4 checkpoint is sub- sequently used to initialize a second training stage with an eight-step rollout horizon (N r = 8). Within each training window, the predicted fluid and struc- tural states are recursively fed back to the model without replacement by the corresponding reference states. All prediction and force-reconstruction results reported in this study are obtained using the selected R4-to-R8 checkpoint. The training objective is defined as the sum of the fluid and structural mean- squared errors over the N r recursively predicted steps. Thus, N r = 4 in the initial training stage and N r = 8 in the subsequent fine-tuning stage: L =L f +L s ,(41) where L f = 1 N r N r X j=1 MSE b e Q n+j f , e Q n+j f , L s = 1 N r N r X j=1 MSE c X n+j s ,X n+j s . (42) The structural loss is evaluated jointly over the two nodal-coordinate com- ponents and the two structural-velocity components.The two loss components 21 are assigned equal weights in the present model. 3.2. Differentiable Aerodynamic Force Reconstruction via DMT and SDF Standard aerodynamic-force evaluation relies on direct integration of the pressure and viscous tractions over the moving solid boundary ∂B. For neural- predicted flow fields on a relatively coarse Eulerian grid, this evaluation can be sensitive to inaccuracies in the near-wall pressure and velocity gradients. We therefore employ the DMT formulation [19] on a fixed closed control contour Σ enclosing the plate in the better-resolved outer-flow region. In its general vorticity-moment form, the total aerodynamic force comprises the time derivative of the vorticity moment within the enclosed domain Ω c and a boundary flux integral along Σ: F =−ρ d dt Z Ω c r× (ωk)dΩ + I Σ T(u,ω, r, n)ds(43) where r is the position vector, k is the out-of-plane unit vector, and T collec- tively represents the pressure, viscous-vorticity, and momentum-flux contribu- tions on the control contour. A critical challenge in implementing Eq.43 within an end-to-end neural net- work is rendering the curve integral differentiable with respect to the Eulerian grid representation. Explicitly tracking the curve Σ breaks grid continuity. To achieve a fully differentiable line integral, we introduce a continuous Signed Dis- tance Function (SDF) following the implicit boundary treatment methodology thoroughly established in deep learning fluid surrogates [23]. Here, an SDF field φ(x) is constructed to implicitly represent the integration contour Σ: φ(x) = −d(x, Σ) if x∈ Ω c 0if x∈ Σ d(x, Σ) if x∈ Ω− Ω c (44) where d(x, Σ) denotes the Euclidean distance from a point x to the curve Σ. The outward unit normal vector is rigorously obtained as n =∇φ/∥∇φ∥. 22 We replace the hard line integral over Σ with an area integral over the Eulerian domain Ω using a regularized Heaviside function H ε and its derivative, the smoothed Dirac-delta function δ ε (φ). To remain consistent with the discrete implementation, a compactly supported cosine regularization is employed: H ε (φ) = 0,φ <−ε, 1 2 1 + φ ε + 1 π sin πφ ε , |φ|≤ ε, 1,φ > ε, δ ε (φ) = ∂H ε ∂φ = 1 2ε 1 + cos πφ ε , |φ|≤ ε, 0,|φ| > ε. (45) where ε controls the half-width of the regularized narrow band. For all aerodynamic- force results reported in this study, ε = 0.51 max(∆x, ∆y) = 0.0051. Consequently, a line integral of a generic integrand I over the curve Σ is approximated by a differentiable two-dimensional area integral weighted by the regularized Dirac-delta function: I Σ I ds≈ Z Ω Iδ ε (φ)∥∇φ∥ dΩ.(46) Applying this regularized SDF-based transformation to the general vorticity- moment equation gives the following approximate area-integral representation: F≈−ρ d dt Z Ω r× (ωk)H ε (−φ)dΩ + Z Ω Tδ ε (φ)∥∇φ∥dΩ(47) The explicit-contour and SDF-based formulations retain the same physical terms in the DMT force identity, differing only in how the control-contour inte- grals are represented numerically. This formulation seamlessly maps the curve integral onto the 2D Eulerian grid predicted by the CNN. Since spatial integration on a fixed grid and point- wise multiplication with δ ε (φ) are inherently tensor operations in modern deep learning frameworks, the entire force reconstruction process becomes rigorously differentiable. Error gradients from force-based loss functions can propagate 23 −1 0 1 y / L −1 0 1 y / L −101 −1 0 1 y / L 0.00.8 0.000.040.08 −101 −0.80.00.8 0.000.050.10 −101 −1.50.0 0.000.060.13 −101 0.00.8 0.000.020.05 −101 −0.80.00.8 0.000.040.08 −101 −3.0−1.50.0 0.000.050.09 −101 −101 0.000.020.04 −101 −0.80.00.8 0.000.020.04 −101 −20 0.000.050.10 (a) κ=0.05: Deflected-flapping uv p x/L (b) κ=0.10: Deflected uv p x/L (c) κ=0.30: Flapping uv p x/L Ground truth Prediction Absolute error Rollout step 59; t=279.6 Reference platePredicted plate Figure 4: Flow-field reconstruction at rollout step 59 for three representative response regimes: (a) κ = 0.05; (b) κ = 0.10; and (c) κ = 0.30. Within each group, the columns show u, v, and p, and the rows show the CFD ground truth, neural prediction, and pointwise absolute error. back through the softened Dirac-delta weighting seamlessly to both the CNN fluid tokens and the Transformer’s solid structural representation. 4. Results and Discussion 4.1. Prediction across stiffness-dependent response regimes 4.1.1. Representative regimes from one stiffness-conditioned operator Following the physical classification introduced in Section 2, three stiffnesses are selected to represent the distinct responses contained in the data set: κ = 0.05 (deflected-flapping), κ = 0.10 (deflected), κ = 0.30 (flapping). (48) All three trajectories are produced by the same trained operator; only the stiffness-conditioning token and the initial coupled state are changed. The first state is supplied from the corresponding CFD trajectory, after which the fluid fields and the positions and velocities of the 101 plate nodes are advanced re- cursively for 59 steps. Figure 4 shows that the predicted plate geometry and the principal high- and low-speed regions remain aligned with the CFD solution after repeated 24 autoregressive reuse. The distinction among the three regimes is not imposed by a discrete class label. It emerges from conditioning the common evolution map on κ and from recursively coupling the resulting structural motion back into the fluid representation. In particular, the model retains the biased mean configuration at κ = 0.05, the small-amplitude deflected response at κ = 0.10, and the large-amplitude nearly symmetric flapping at κ = 0.30. The field and structural errors at rollout step n are measured as e n q = ∥bq n − q n ∥ 2 ∥q n ∥ 2 , q ∈u,v,p,(49) and e n s = " 1 2N s N s X i=1 b x n s,i −x n s,i 2 2 # 1/2 .(50) Here, x n s,i = (x n i ,y n i ) T and b x n s,i = (bx n i ,by n i ) T denote the reference and predicted positions of the ith Lagrangian point, respectively. The quantity e n q is the relative L 2 error over the Eulerian grid, whereas e n s is the coordinate RMSE over the ordered Lagrangian points. As shown in Fig. 5, the errors remain bounded over the evaluated interval, although their magnitude depends on both the variable and the response regime. Averaged over the 59 predicted states, the relative errors are 1.89%, 5.39%, and 5.41% for u, v, and p, respectively, at κ = 0.05; 1.41%, 4.08%, and 3.36% at κ = 0.10; and 1.13%, 2.63%, and 3.94% at κ = 0.30. The corresponding mean plate-coordinate RMSE values are 2.91× 10 −3 , 9.17× 10 −4 , and 1.53× 10 −3 , respectively. The plate-coordinate errors remain of the same order across the three response regimes, indicating that the structural evolution is preserved even in the large-amplitude flapping case. 4.1.2. Accuracy across all 36 training stiffnesses The representative cases provide a qualitative view of the three regimes, whereas the complete stiffness sweep tests whether the same operator behaves consistently throughout the parameter set. Figure 6 summarizes the mean and final-step errors for all 36 training stiffnesses. 25 0.00 0.05 0.10 0.15 0.20 Relative L 2 error (a) = 0.05, Deflected-flapping: fluid uvp 0.000 0.001 0.002 0.003 0.004 0.005 0.006 0.007 Coordinate RMSE (b) = 0.05, Deflected-flapping: plate 0.00 0.05 0.10 0.15 0.20 Relative L 2 error (c) = 0.10, Deflected: fluid 0.000 0.001 0.002 0.003 0.004 0.005 0.006 0.007 Coordinate RMSE (d) = 0.10, Deflected: plate 01020304050 Rollout step 0.00 0.05 0.10 0.15 0.20 Relative L 2 error (e) = 0.30, Flapping: fluid 01020304050 Rollout step 0.000 0.001 0.002 0.003 0.004 0.005 0.006 0.007 Coordinate RMSE (f) = 0.30, Flapping: plate Figure 5: Stepwise rollout errors for the three representative response regimes over 59 predicted states. The left column shows the relative L 2 errors of u, v, and p, while the right column shows the coordinate RMSE of the plate. The rows correspond to κ = 0.05, 0.10, and 0.30, respectively. Hereafter, the plate-coordinate RMSE defined in Eq. (50) is denoted by e n xy . Since the present model also predicts the two structural velocity components, the corresponding plate-velocity RMSE is defined as e n uv = 1 2N s N s X i=1 bu n s,i bv n s,i − u n s,i v n s,i 2 2 1/2 .(51) Across the 36 training conditions, the maximum mean relative errors are approximately 2.19% for u, 6.72% for v, and 5.72% for p. The maxima for u and v occur at κ = 0.09, whereas the maximum mean pressure error occurs at 26 Table 2: Mean errors over the 59 predicted states for eight representative stiffnesses selected from the 36 training conditions. Fluid errors are relative L 2 errors in percent; e xy and e uv are the nondimensional RMSEs of the plate position and velocity, respectively. κ Regimee u (%) e v (%) e p (%)e xy e uv 0.05 Deflected–flapping 1.895.39 5.41 2.91× 10 −3 4.35× 10 −3 0.10 Deflected1.414.08 3.36 9.17× 10 −4 1.68× 10 −3 0.15 Flapping1.043.26 4.21 1.19× 10 −3 1.54× 10 −3 0.20 Flapping0.862.12 3.54 7.75× 10 −4 1.09× 10 −3 0.25 Flapping0.872.20 3.63 6.15× 10 −4 8.10× 10 −4 0.30 Flapping1.132.63 3.94 1.53× 10 −3 2.23× 10 −3 0.35 Flapping1.513.58 4.28 1.51× 10 −3 2.56× 10 −3 0.40 Flapping1.863.79 5.06 3.02× 10 −3 4.72× 10 −3 κ = 0.13. The maximum mean plate-position and plate-velocity RMSE values are 4.45 × 10 −3 and 6.22 × 10 −3 , respectively, both occurring at κ = 0.31. Because the reported errors compare time-aligned snapshots point by point, an accumulated phase difference between the predicted and reference oscillations can produce a relatively large instantaneous error without implying a qualitative failure of the predicted dynamics. 27 0 1 2 3 4 5 6 7 8 Relative L 2 error (\%) (a) Mean fluid-field error over 59-step rollout uvp 0 2 4 6 8 10 12 14 Relative L 2 error (\%) (b) Fluid-field error at rollout step 59 0.000 0.001 0.002 0.003 0.004 0.005 Position RMSE, e xy / L (c) Mean flexible-plate position error 0.000 0.001 0.002 0.003 0.004 0.005 0.006 Position RMSE, e xy / L (d) Plate position error at rollout step 59 0.050.100.150.200.250.300.350.40 Bending stiffness, 0.000 0.001 0.002 0.003 0.004 0.005 0.006 0.007 Velocity RMSE, e uv (e) Mean flexible-plate velocity error 0.050.100.150.200.250.300.350.40 Bending stiffness, 0.000 0.002 0.004 0.006 0.008 0.010 Velocity RMSE, e uv (f) Plate velocity error at rollout step 59 Figure 6: Prediction errors across the 36 training stiffnesses. (a) Mean fluid-field relative L 2 error over the 59 predicted states. (b) Fluid-field relative L 2 error at rollout step 59. (c) Mean flexible-plate position RMSE over the 59 predicted states. (d) Plate-position RMSE at rollout step 59. (e) Mean flexible-plate velocity RMSE over the 59 predicted states. (f) Plate-velocity RMSE at rollout step 59. The dashed separators indicate the deflected–flapping, deflected, and flapping regimes. 28 Despite these localized error variations, the neural operator reproduces the overall oscillatory evolution of the coupled system and preserves the principal spatial organization of the velocity and pressure fields over the evaluated roll- out. More importantly, the same stiffness-conditioned operator successfully cap- tures the deflected–flapping, deflected, and flapping response regimes, including their distinct mean plate configurations, oscillation amplitudes, and wake struc- tures. Overall, these results support the interpretation of the network as a single stiffness-conditioned evolution operator rather than a collection of inde- pendently fitted response curves. 4.2. Oscillation frequency, phase-space structure, and long rollout Pointwise field errors alone do not establish whether an autoregressive model preserves the underlying structural dynamics. A small phase drift can increase the instantaneous field error while leaving the limit-cycle topology and dominant frequency nearly unchanged. We therefore examine the free-tip displacement in both the time and frequency domains, as shown in Fig. 7. The dominant frequencies are recovered exactly at the resolution of the sam- pled sequence: f = 0.1333 for κ = 0.05, f = 0.4333 for κ = 0.10, and f = 0.1333 for κ = 0.30. The predicted spectral peak amplitudes are also close to the CFD values. In particular, the model distinguishes the higher-frequency, small- amplitude deflected response from the two lower-frequency flapping responses, rather than reproducing a common oscillation for all stiffnesses. The reconstructed limit cycles are further examined using the free-tip phase portraits shown in Fig. 8. The range-normalized displacement and velocity RMSE values are 1.27% and 1.44% for κ = 0.05, 2.71% and 1.82% for κ = 0.10, and 0.20% and 0.22% for κ = 0.30. The larger normalized displacement error at κ = 0.10 should be interpreted together with its very small displacement range: normalization by this range magnifies modest absolute deviations. The predicted orbit nevertheless retains the location, orientation, and periodic structure of the CFD cycle. 29 0.4 0.5 0.6 0.7 0.8 0.9 y tip / L (a) κ=0.05, Deflected-flapping: tip displacement Ground truthPrediction 0.000 0.025 0.050 0.075 0.100 0.125 0.150 Amplitude, A y / L (b) κ=0.05, Deflected-flapping: spectrum −0.87 −0.86 −0.85 −0.84 −0.83 −0.82 y tip / L (c) κ=0.10, Deflected: tip displacement 0.000 0.005 0.010 0.015 0.020 Amplitude, A y / L (d) κ=0.10, Deflected: spectrum 255260265270275 Time, t −1.0 −0.5 0.0 0.5 1.0 y tip / L (e) κ=0.30, Flapping: tip displacement 0.00.20.40.60.81.0 Frequency, f 0.0 0.2 0.4 0.6 0.8 1.0 Amplitude, A y / L (f) κ=0.30, Flapping: spectrum Figure 7: Free-tip displacement histories and frequency spectra for the three representative response regimes. The left column compares the predicted and CFD free-tip displacements, while the right column shows their corresponding amplitude spectra. The rows represent (a,b) κ = 0.05, (c,d) κ = 0.10, and (e,f) κ = 0.30. Finally, a blind autoregressive rollout over 1000 steps was performed to ex- amine long-horizon stability. The predicted coupled state remained bounded and preserved the regime-dependent periodic behavior over the full rollout. Since no synchronized CFD reference is available over this extended interval, this test is interpreted as a stability diagnostic rather than a direct measure of 1000-step prediction accuracy. Detailed stability diagnostics are provided in Appendix A. 30 0.40.60.8 y tip /L −0.4 −0.2 0.0 0.2 0.4 v tip =d( y tip / L )/d t NRMSE y = 1.27%, NRMSE v = 1.44% (a) κ=0.05, Deflected-flapping Ground truth Prediction −0.86−0.84−0.82 y tip /L −0.06 −0.04 −0.02 0.00 0.02 0.04 0.06 NRMSE y = 2.71%, NRMSE v = 1.82% (b) κ=0.10, Deflected −1.0−0.50.00.51.0 y tip /L −1.0 −0.5 0.0 0.5 1.0 NRMSE y = 0.20%, NRMSE v = 0.22% (c) κ=0.30, Flapping Figure 8: Free-tip phase portraits for three representative response regimes: (a) κ = 0.05, deflected–flapping; (b) κ = 0.10, deflected; and (c) κ = 0.30, flapping. The transverse tip velocity is obtained from the time derivative of the free-tip displacement. The normalized RMSE values shown above each panel quantify the displacement and velocity differences between the prediction and CFD reference. 4.3. Interpolation at two unseen stiffnesses The 36 training stiffnesses demonstrate that a single stiffness-conditioned operator can represent the three response regimes over the sampled parameter interval. To examine interpolation between the discrete training conditions, two additional CFD trajectories are considered at κ = 0.105 and κ = 0.335. These stiffnesses were excluded from training and lie midway between the neighboring training pairs (0.10, 0.11) and (0.33, 0.34), respectively. No model parameters are updated for either test. The two cases represent different interpolation conditions: κ = 0.335 lies within the established flapping regime, whereas κ = 0.105 lies in the deflected regime, where symmetry-related plate configurations can occur. 4.3.1. Interpolation within the flapping regime at κ = 0.335. Figure 9 compares the predicted and CFD flow fields at the final rollout step. The model preserves the large-scale flapping pattern, the instantaneous plate configuration, and the principal spatial organization of the velocity and pressure fields. Although localized errors remain near the moving plate and the convected wake structures, the predicted coupled state remains qualitatively consistent with the CFD solution after 59 recursive steps. 31 −1.0 −0.5 0.0 0.5 1.0 y / L (a) Ground truth u −1.0 −0.5 0.0 0.5 1.0 y / L (d) Prediction −1.0−0.50.00.51.01.5 x/L −1.0 −0.5 0.0 0.5 1.0 y / L (g) Absolute error Reference plate Predicted plate −1.0−0.50.00.51.01.5 0.000.010.020.030.040.050.060.07 (b) Ground truth v (e) Prediction −1.0−0.50.00.51.01.5 x/L (h) Absolute error −1.0−0.50.00.51.0 0.000.020.040.060.08 (c) Ground truth p (f) Prediction −1.0−0.50.00.51.01.5 x/L (i) Absolute error −3−2−101 0.0000.0250.0500.0750.1000.1250.150 κ=0.335, Unseen interpolation; rollout step 59, t=279.6 Figure 9: Flow-field reconstruction at rollout step 59 for the unseen stiffness κ = 0.335. The columns show u, v, and p, while the rows show the CFD ground truth, neural prediction, and pointwise absolute error. The temporal evolution of the prediction errors is shown in Fig. 10. Averaged over the 59 predicted states, the relative errors in u, v, and p are 2.20%, 5.20%, and 6.19%, respectively. The mean plate-position and plate- velocity RMSE values are 5.13× 10 −3 and 7.29× 10 −3 , respectively. At rollout step 59, the corresponding fluid-field errors are 1.97%, 3.79%, and 5.57%, while the plate-position and plate-velocity RMSE values are 7.34× 10 −3 and 8.77× 10 −3 . Both the fluid-field and structural errors remain bounded throughout the evaluated interval, indicating stable interpolation of the coupled Eulerian and Lagrangian states. 32 1020304050 Rollout step 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 Relative L 2 error (a) Unseen κ=0.335: fluid fields ΔuΔvΔp 1020304050 Rollout step 0.000 0.002 0.004 0.006 0.008 0.010 0.012 0.014 Coordinate RMSE (b) Unseen κ=0.335: flexible plate Figure 10: Evolution of the fluid-field relative L 2 errors and plate-position RMSE over the 59-step autoregressive rollout at the unseen stiffness κ = 0.335. The corresponding low-dimensional structural dynamics are presented in Fig. 11. 260270 Time, t −1.0 −0.5 0.0 0.5 1.0 y tip / L (a) Tip displacement Ground truth Prediction 0.00.20.40.60.81.0 Frequency, f 0.0 0.2 0.4 0.6 0.8 1.0 Amplitude, A y / L f GT =0.133 f pred =0.133 (b) Tip-displacement spectrum −1.0−0.50.00.51.0 y tip /L −1.0 −0.5 0.0 0.5 1.0 v tip =d( y tip / L )/d t NRMSE y = 0.75%, NRMSE v = 0.73% (c) Tip phase portrait Unseen stiffness κ=0.335: structural dynamics Figure 11: Structural dynamics at the unseen stiffness κ = 0.335: free-tip displacement, Fourier spectrum, and phase portrait. CFD values are shown by open circles and neural predictions by solid lines. The predicted and CFD free-tip signals have the same dominant frequency, f = 0.1333. The range-normalized displacement and velocity phase-space RMSE values are 0.75% and 0.73%, respectively. The predicted orbit therefore retains the amplitude, orientation, and periodic structure of the CFD limit cycle. 4.3.2. Interpolation within the deflected regime at κ = 0.105. Interpolation at κ = 0.105 is more sensitive because the plate response has a small oscillation amplitude and symmetry-related deflected configurations occur 33 in this stiffness range. Nevertheless, as shown in Fig. 12, the model recovers the upward-deflected mean configuration of the CFD trajectory and retains the principal near-plate and wake structures. 1.0 0.5 0.0 0.5 1.0 y / L (a) Ground truth u 1.0 0.5 0.0 0.5 1.0 y / L (d) Prediction 1.00.50.00.51.01.5 x/L 1.0 0.5 0.0 0.5 1.0 y / L (g) Absolute error Reference plate Predicted plate 0.50.00.51.0 0.000.020.040.060.080.100.120.14 (b) Ground truth v (e) Prediction 1.00.50.00.51.01.5 x/L (h) Absolute error 1.00.50.00.51.0 0.000.050.100.150.20 (c) Ground truth p (f) Prediction 1.00.50.00.51.01.5 x/L (i) Absolute error 3.02.52.01.51.00.50.00.5 0.000.050.100.150.200.25 = 0.105, Unseen interpolation; rollout step 59, t = 279.6 Figure 12: Flow-field reconstruction at rollout step 59 for the unseen stiffness κ = 0.105. The columns show u, v, and p, while the rows show the CFD ground truth, neural prediction, and pointwise absolute error. The pointwise discrepancies are larger than those at κ = 0.335, particularly in the transverse velocity and pressure fields. The predicted plate and flow fields nevertheless remain on the same deflected branch as the CFD reference. The error histories in Fig. 13 remain bounded, although their magnitudes are higher than those obtained at κ = 0.335. Averaged over the 59 predicted states, the relative errors in u, v, and p are 34 1020304050 Rollout step 0.00 0.05 0.10 0.15 0.20 Relative L 2 error (a) Unseen κ=0.105: fluid fields ΔuΔvΔp 1020304050 Rollout step 0.000 0.005 0.010 0.015 0.020 Coordinate RMSE (b) Unseen κ=0.105: flexible plate Figure 13: Evolution of the fluid-field relative L 2 errors and plate-position RMSE over the 59-step autoregressive rollout at the unseen stiffness κ = 0.105. The errors remain bounded, although their magnitudes are higher than those obtained within the flapping regime. 4.05%, 11.17%, and 10.28%, respectively. The mean plate-position and plate- velocity RMSE values are 1.22× 10 −2 and 7.98× 10 −3 , respectively. At rollout step 59, the corresponding fluid-field errors are 5.40%, 19.70%, and 11.05%, while the plate-position and plate-velocity RMSE values are 1.03× 10 −2 and 5.08×10 −3 . The absence of monotonic error growth indicates that the increased errors do not arise from numerical divergence of the autoregressive rollout. The structural response is examined in Fig. 14. 260270 Time, t 0.82 0.83 0.84 0.85 0.86 0.87 y tip / L (a) Tip displacement Ground truth Prediction 0.00.20.40.60.81.0 Frequency, f 0.000 0.005 0.010 0.015 Amplitude, A y / L f GT =0.433 f pred =0.433 (b) Tip-displacement spectrum 0.820.840.86 y tip /L −0.06 −0.04 −0.02 0.00 0.02 0.04 0.06 v tip =d( y tip / L )/d t NRMSE y = 17.18%, NRMSE v = 17.38% (c) Tip phase portrait Unseen stiffness κ=0.105: structural dynamics Figure 14: Structural dynamics at the unseen stiffness κ = 0.105: free-tip displacement, Fourier spectrum, and phase portrait. Although the dominant frequency and mean deflected branch are recovered, differences remain in the oscillation amplitude and phase-space trajec- tory. The model recovers the CFD dominant frequency of f = 0.4333 and retains 35 a repeatable oscillatory motion. The range-normalized displacement and ve- locity phase-space RMSE values are 17.18% and 17.38%, respectively. Thus, the model captures the mean deflected branch and the characteristic structural frequency, although differences remain in the oscillation amplitude and phase. The interpolation at κ = 0.105 is therefore less accurate than that at κ = 0.335, but it does not exhibit an incorrect branch orientation or an unstable rollout. Taken together, the two unseen-stiffness cases show that the same stiffness- conditioned operator can interpolate both within the flapping regime and within the branch-sensitive deflected regime. The interpolation at κ = 0.335 accu- rately preserves the flow-field organization and structural limit cycle, while the case at κ = 0.105 retains the correct mean branch and dominant frequency with larger instantaneous errors. These results support the use of a single con- ditional evolution operator across multiple stiffness-dependent FSI responses, while also identifying reduced instantaneous accuracy in the small-amplitude deflected regime. 4.4. Differentiable aerodynamic-force reconstruction from predicted flow fields The preceding results show that the coupled neural operator can maintain a coherent moving interface while evolving the surrounding velocity and pres- sure fields. Reconstructing the flow and structural states, however, does not by itself provide a complete description of the fluid–structure dynamics. The aero- dynamic loading is the quantity through which the fluid drives the structural motion, and the time-dependent drag and lift coefficients therefore provide a stricter physical assessment of the predicted coupled state. Unlike pointwise flow-field errors, aerodynamic forces are global function- als of the solution. They depend simultaneously on pressure, velocity gradients, vorticity, momentum transport, and their temporal variation. Small local errors that are barely visible in an instantaneous flow-field plot may accumulate co- herently during force integration. A predicted field may therefore reproduce the main wake structures while still yielding inaccurate aerodynamic coefficients. 36 The conventional evaluation of the aerodynamic force is based on the pres- sure and viscous tractions acting directly on the instantaneous body surface, F = Z ∂B (−pn +τ ·n) ds,(52) where ∂B denotes the deforming plate boundary. This formulation requires accurate near-wall pressure and velocity gradients. These quantities are partic- ularly sensitive to the spatial resolution of the neural prediction, because the thin boundary layer is only partially resolved on the coarse Eulerian grid. Direct traction integration consequently amplifies near-interface smoothing and local differentiation errors. To reduce this dependence on the under-resolved near-wall region, the DMT is employed to express the aerodynamic force using quantities evaluated over a control domain and its enclosing contour [19]. In the present two-dimensional setting, the control surface reduces to a closed control line surrounding the plate. The force evaluation is thereby transferred from the body surface to the better-resolved outer-flow region. 4.4.1. Fixed control contours for the flexible-plate motion The control contour must enclose the plate throughout the complete au- toregressive rollout. A contour that follows the plate at every time step would introduce a time-dependent integration geometry and complicate both the nu- merical evaluation and subsequent differentiation. A single fixed contour is therefore constructed separately for each stiffness. Denoting the i-th predicted Lagrangian point at rollout step n by (x n i ,y n i ), the contour is obtained from the envelope of all predicted plate configurations, x 0 = min n,i (x n i )− m, x 1 = max n,i (x n i ) + m, y 0 = min n,i (y n i )− m, y 1 = max n,i (y n i ) + m, m = 0.3,(53) after which the four sides are aligned with the Eulerian grid. The margin sep- arates the contour from the immediate vicinity of the plate, while the fixed 37 1.00.50.00.51.01.5 x/L 1.0 0.5 0.0 0.5 1.0 y / L (a) = 0.05, Deflected-flapping Sampled rollout plates Plate at step 59 Fixed Level-set DMT surface 1.00.50.00.51.01.5 x/L 1.0 0.5 0.0 0.5 1.0 (b) = 0.10, Deflected 1.00.50.00.51.01.5 x/L 1.0 0.5 0.0 0.5 1.0 (c) = 0.30, Flapping 0.50 0.25 0.00 0.25 0.50 0.75 1.00 1.25 1.50 Predicted streamwise velocity, u Figure 15: Fixed DMT control contours for the three representative flexible-plate regimes. The background shows the predicted streamwise velocity at rollout step 59. The white curves show sampled plate configurations during the rollout, the black curve denotes the plate at step 59, and the orange rectangle is the fixed control contour used for force reconstruction. geometry ensures that the same integration operator is used at every rollout step. Figure 15 illustrates the control contours for the deflected-flapping, deflected, and flapping regimes. The plate occupies substantially different regions of the domain in the three cases, so a separate contour is required for each stiffness. Within each case, however, the contour remains fixed and encloses the full struc- tural trajectory. This construction separates the force-integration geometry from the instantaneous plate motion and provides a common basis for the ex- plicit and Level-set DMT realizations. 4.4.2. Isolation of the DMT reconstruction error Before applying DMT to the neural-predicted fields, the force operator is evaluated using the original CFD flow fields. This intermediate validation is necessary because a discrepancy in the final force history may originate either from the DMT reconstruction or from the neural flow prediction. Supplying resolved CFD fields removes the second source and isolates the numerical error of the control-contour formulation. The validation is performed for the flapping case at κ = 0.30 and Re = 200, using 501 consecutive CFD fields at the original temporal resolution. The pressure-containing explicit Line DMT operator is evaluated on the fixed control 38 contour, and a cubic spline is used to differentiate the moment term in time. The resulting coefficients are compared with the aerodynamic forces reported directly by the CFD solver. The error metrics are computed using the 499 interior samples after excluding the two endpoint samples. The range-normalized root-mean-square error is defined as NRMSE = r 1 N P N n=1 b C n − C n 2 max(C)− min(C) × 100%,(54) where C denotes the direct CFD coefficient and b C its DMT reconstruction. This range-based normalization avoids the singular behavior of pointwise relative errors when the lift coefficient crosses zero. 1 2 3 4 C D NRMSE= 0.279%, R= 0.999991 (a) Drag coefficient Direct CFD forceDMT reconstruction 250260270280290300 Time, t −2 −1 0 1 2 C L NRMSE= 0.118%, R= 0.999991 (b) Lift coefficient Figure 16: Validation of the pressure-containing explicit Line DMT reconstruction using the original CFD flow fields at κ = 0.30. Open circles denote the direct CFD force and solid lines denote the DMT reconstruction: (a) drag coefficient and (b) lift coefficient. As shown in Fig. 16, the DMT reconstruction is nearly indistinguishable from the direct CFD force. The drag coefficient gives an NRMSE of 0.279% and a Pearson correlation coefficient of 0.999991, while the lift coefficient gives an NRMSE of 0.118% with the same correlation coefficient. The agreement 39 extends to both the principal extrema and the smaller secondary variations of the force histories. This result establishes that the explicit Line DMT formulation can recover the aerodynamic loading without direct access to the near-wall traction. More importantly, it provides a reference level for interpreting the neural-rollout re- sults: force errors substantially larger than the sub-percent CFD-field recon- struction error can be attributed primarily to inaccuracies in the predicted flow fields rather than to the underlying DMT formulation. 4.4.3. Aerodynamic loading recovered from neural rollouts Having isolated the DMT reconstruction error, the same force operator is applied to the autoregressively predicted u, v, and p fields. Three stiffnesses are selected to represent the distinct structural responses contained in the data set, κ = 0.05 (deflected-flapping), κ = 0.10 (deflected), κ = 0.30 (flapping). (55) For each case, only the initial coupled state is supplied by CFD; all subse- quent fluid and structural states are generated recursively by the neural opera- tor. The force reconstruction therefore evaluates whether the predicted sequence retains the global pressure and momentum organization required to reproduce the aerodynamic loading. Two realizations are considered. The explicit Line DMT samples and inte- grates the required quantities directly along the grid-aligned control contour. The Level-set DMT represents the same contour implicitly and replaces its line contribution with a regularized Eulerian-band integral. Both realizations receive the same neural-predicted flow sequence. Figure 17 shows that the force histories remain qualitatively different across the three structural regimes, while their principal temporal features are con- sistently recovered from the neural-predicted fields. The reconstructed curves preserve the dominant peak locations, oscillation periods, mean loading levels, and phase relationships of the direct CFD forces. 40 0.75 1.00 1.25 1.50 1.75 2.00 2.25 2.50 C D RMSE: line 0.058, level-set 0.076 (a) = 0.05, Deflected-flapping: C D Direct CFD force Line DMT Level-set DMT 1.50 1.25 1.00 0.75 0.50 0.25 0.00 C L RMSE: line 0.066, level-set 0.142 (b) = 0.05, Deflected-flapping: C L 1.6 1.8 2.0 2.2 2.4 C D RMSE: line 0.027, level-set 0.043 (c) = 0.10, Deflected: C D 0.4 0.2 0.0 0.2 C L RMSE: line 0.025, level-set 0.056 (d) = 0.10, Deflected: C L 250255260265270275280 Time, t 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 C D RMSE: line 0.052, level-set 0.088 (e) = 0.30, Flapping: C D 250255260265270275280 Time, t 2 1 0 1 2 C L RMSE: line 0.163, level-set 0.208 (f) = 0.30, Flapping: C L Figure 17: Aerodynamic-force reconstruction from the autoregressively predicted flow fields for three representative structural regimes. The rows correspond to κ = 0.05 (deflected-flapping), κ = 0.10 (deflected), and κ = 0.30 (flapping); the columns show C D and C L . Open circles denote the direct CFD force, dashed blue lines denote the explicit Line DMT reconstruction, and solid red lines denote the Level-set DMT reconstruction. For the explicit Line DMT, the NRMSE values range from 1.52% to 5.34% across the six aerodynamic-coefficient histories. This agreement indicates that the flow fields predicted by the stiffness-conditioned operator retain sufficient outer-flow information to recover the principal aerodynamic response. The deflected-flapping case produces a comparatively complex and asym- metric loading history. Although both DMT realizations recover its major force variations, the Level-set lift error reaches 11.49%, which is the largest error 41 Table 3: RMSE and range-normalized RMSE (NRMSE) of the aerodynamic coefficients re- constructed from the neural-predicted flow fields. κ RegimeCoefficientLine DMTLevel-set DMT RMSE NRMSE (%) RMSE NRMSE (%) 0.05 Deflected-flapping C D 0.05763.740.07654.96 0.05 Deflected-flappingC L 0.06635.340.142411.49 0.10 DeflectedC D 0.02663.270.04275.25 0.10 DeflectedC L 0.02493.650.05588.19 0.30 FlappingC D 0.05211.520.08842.58 0.30 FlappingC L 0.16263.450.20754.41 among the representative cases. The lift signal in this regime contains several closely spaced extrema and relatively rapid variations, making it sensitive to ac- cumulated flow-field errors and to the finite-width spatial averaging introduced by the Level-set band. For the deflected and flapping cases, the dominant loading periods and phases remain well preserved. In particular, the two DMT realizations give similar results for the large-amplitude flapping case at κ = 0.30, despite the substantial plate excursion. This agreement indicates that the coherent peri- odic loading is encoded robustly in the predicted outer-flow dynamics rather than depending exclusively on the under-resolved near-wall region. The CFD-field validation in Fig. 16 establishes a sub-percent reference error for the explicit Line DMT at κ = 0.30. The larger Line DMT errors obtained from the neural rollouts therefore primarily reflect accumulated inaccuracies in the predicted pressure, velocity, vorticity, and temporal momentum evolution. The additional difference between the Line and Level-set results also reflects the finite-width regularization used in the Level-set integration. The force compar- ison consequently provides a more integrated physical assessment than instan- taneous field errors because it tests whether the predicted coupled trajectory remains consistent with the global aerodynamic loading. 42 These results demonstrate that the neural operator does not merely re- produce visually similar flow structures. Across three distinct plate-response regimes, the predicted fields preserve sufficient global information to reconstruct the dominant drag and lift dynamics. The explicit Line DMT establishes the fea- sibility of force recovery from the predicted coarse-grid fields, while the Level-set realization provides a differentiable approximation of the same control-contour formulation. Aerodynamic loading at the unseen stiffnesses. The same force-reconstruction procedure is applied to the two unseen stiff- nesses examined in Section 4.3. No model parameters are updated, and both DMT realizations operate on the recursively predicted flow fields. At κ = 0.335, the reconstructed forces closely follow the direct CFD results throughout the evaluated interval. 0 1 2 3 4 C D RMSE: line 0.089, level-set 0.119 (a) Drag coefficient reconstructed from the predicted field Line DMTLevel-set DMTDirect CFD 250255260265270275280 Physical time, t −2 −1 0 1 2 C L RMSE: line 0.141, level-set 0.187 (b) Lift coefficient reconstructed from the predicted field Figure 18: Drag and lift coefficients reconstructed from the predicted flow fields at the unseen stiffness κ = 0.335, compared with the forces reported directly by the CFD solver. Open circles denote the direct CFD force, dashed blue lines denote the explicit Line DMT reconstruction, and solid red lines denote the Level-set DMT reconstruction. For the explicit Line DMT, the drag and lift NRMSE values are 2.38% and 43 3.00%, respectively. The corresponding Level-set DMT errors are 3.17% and 3.99%. Both realizations retain the principal extrema, oscillation period, and phase of the direct CFD forces. This agreement confirms that the interpolated flow fields preserve the global pressure and momentum organization required to recover the aerodynamic loading within the flapping regime. The force reconstruction at κ = 0.105 is more challenging, consistent with the larger flow-field and structural errors reported in Section 4.3. 1.6 1.8 2.0 2.2 2.4 C D RMSE: line 0.067, level-set 0.070 (a) Drag coefficient reconstructed from the predicted field Line DMTLevel-set DMTDirect CFD 250255260265270275280 Physical time, t −0.4 −0.2 0.0 0.2 0.4 0.6 C L RMSE: line 0.084, level-set 0.100 (b) Lift coefficient reconstructed from the predicted field Figure 19: Drag and lift coefficients reconstructed from the predicted flow fields at the unseen stiffness κ = 0.105, compared with the forces reported directly by the CFD solver. Open circles denote the direct CFD force, dashed blue lines denote the explicit Line DMT reconstruction, and solid red lines denote the Level-set DMT reconstruction. For the explicit Line DMT, the drag and lift NRMSE values are 9.00% and 11.73%, respectively, while the corresponding Level-set DMT values are 9.50% and 14.07%. Although these errors are larger than those at κ = 0.335, the reconstructed histories retain the principal periodic variations of the direct CFD forces. The larger discrepancies are consistent with the higher flow-field and structural errors observed in this branch-sensitive regime, while the bounded rollout errors provide no evidence of autoregressive divergence. 44 4.4.4. Level-set realization and differentiability Explicit control-line integration requires the extraction and ordered traversal of discrete contour points. Although suitable for post-processing, such index- ing operations are inconvenient within a tensor-based automatic-differentiation graph. The fixed contour is therefore represented implicitly as the zero level set of a signed-distance function, Σ =x∈ Ω : φ(x) = 0.(56) A generic contour integral is converted to an Eulerian area integral using a smoothed Dirac-delta function, I Σ I ds≈ Z Ω I δ ε (φ)|∇φ| dΩ,(57) where δ ε (φ) distributes the contour contribution over a narrow band of fi- nite width. Because the contour remains fixed, the signed-distance field, the smoothed delta function, and the normal-vector weights can be precomputed. The force reconstruction is then composed of spatial and temporal differences, pointwise tensor products, and weighted summations, all of which remain in the automatic-differentiation graph. The resulting discrete operator admits ∂(C D ,C L ) ∂(u,v,p) (58) for the complete predicted flow sequence. The differentiable formulation makes the aerodynamic-force reconstruction available as an active functional of the predicted flow field. A conventional non-differentiable post-processing procedure can evaluate lift and drag, but it cannot transmit information on how the upstream flow prediction or a prescribed design variable should change in order to modify these quantities. In contrast, the Level-set formulation retains the dependence of the reconstructed force on every velocity and pressure value participating in the Eulerian integration. When the force operator is composed with a differentiable neural rollout, a force-based objective J can, in principle, transmit sensitivity information 45 through the complete computational chain, dJ dα = ∂J ∂C ∂C ∂q ∂q ∂α , C = (C D ,C L ), q = (u,v,p),(59) where α denotes a differentiable model parameter, conditioning variable, con- trol variable, or design variable. The fixed Level-set contour is advantageous in this setting because its signed-distance field and quadrature weights remain unchanged during the rollout. Consequently, no time-dependent contour extrac- tion or discrete boundary-indexing operation interrupts the gradient path. The regularized Level-set integration requires a balance between smooth- ing and grid resolution. Analytical checks of the contour integral, geometric gradient, Taylor remainder, and translating-point-vortex force are presented in Appendix B. These tests verify the accuracy of the Level-set approximation and the consistency of its automatically differentiated gradients. In the present study, the force operator is applied only as a post-processing diagnostic, and no force-based objective is included in the training loss. The differentiable formulation nevertheless provides the computational connection required for future force-aware learning and gradient-based aeroelastic design. Taken together, the results establish a hierarchy of predictive capability. Across the 36 training stiffnesses, a single stiffness-conditioned operator recon- structs the coupled flow fields and plate motion over three physically distinct response regimes, with bounded 59-step prediction errors and stable long-term autoregressive rollouts. At the unseen stiffness κ = 0.335, the model accurately interpolates the flow-field organization, structural limit cycle, and aerodynamic loading within the flapping regime. The more challenging case at κ = 0.105 ex- hibits larger instantaneous field, structural, and force errors, but the prediction recovers the correct deflected branch and the dominant structural frequency. Together with the aerodynamic-force results, these findings show that the pro- posed model provides a unified representation of stiffness-dependent FSI without requiring a separate surrogate for each response regime. 46 5. Conclusions We have developed and systematically evaluated a stiffness-conditioned neu- ral evolution operator for fluid–structure interaction involving an inverted flex- ible plate. The framework jointly evolves Eulerian flow fields and an explicit Lagrangian structural state, using 101 ordered structural tokens carrying nodal coordinates and velocities together with a global stiffness-conditioning token. Bidirectional fluid–structure coupling enables a single operator to represent the coupled evolution of the moving plate and surrounding flow. Across 36 sampled stiffness values, the model reproduced three distinct re- sponse regimes, namely deflected–flapping, deflected, and flapping dynamics, without requiring a separate surrogate for each condition. The predicted tra- jectories preserved the principal flow organization, structural deformation, os- cillation frequency, and phase-space structure. Blind 1000-step autoregressive rollouts remained bounded and retained the characteristic dynamical behavior of the three regimes, although these tests assess long-term dynamical stability rather than pointwise prediction accuracy over the full horizon. The stiffness-conditioned representation also interpolated to two unseen stiff- ness values. The case at κ = 0.335 preserved the flow-field organization, struc- tural limit cycle, and aerodynamic response within the flapping regime. At the branch-sensitive κ = 0.105, the model recovered the same symmetry-related deflected branch as the CFD reference and retained the dominant oscillatory behavior, although with larger instantaneous errors. In addition, a differentiable aerodynamic-force readout was constructed by combining the derivative-moment transformation with a fixed control contour represented through a signed-distance function and a smoothed Dirac-delta for- mulation. The approach transfers force evaluation from the under-resolved mov- ing boundary to the better resolved outer-flow region while retaining gradients with respect to the predicted flow variables. Validation against resolved CFD fields confirmed the accuracy of the control-contour formulation, and applica- tion to neural rollouts demonstrated that the principal aerodynamic loading can 47 be recovered from coarse predicted flow fields. The differentiable force reconstruction developed here establishes a direct computational pathway from the learned fluid–structure state to aerodynamic loading and its gradients. This pathway can be extended naturally to inverse parameter identification from sparse flow, structural, or force measurements and to gradient-based optimization of stiffness and other design variables. The conditional neural operator also provides a reusable forward model for many- query parameter studies and design iterations, potentially avoiding repeated high-fidelity solutions for every candidate condition. Quantifying this compu- tational advantage and demonstrating end-to-end inverse and design tasks are important next steps, together with broader validation under unseen parameters and incomplete observations. 48 Appendix A. Long-horizon autoregressive stability diagnostics To examine the behavior of the learned evolution operator beyond the 59- step evaluation horizon considered in the main text, blind autoregressive rollouts are performed for 1000 recursive steps at the three representative stiffnesses, κ = 0.05, 0.10, and 0.30. Only the initial coupled state is supplied by CFD, and all subsequent fluid and structural states are generated recursively without intermediate reference-state correction. Figure A.1 presents the corresponding long-horizon stability diagnostics. The relative differences between consecutive predicted fluid states and the co- ordinate RMSE between consecutive predicted plate states are used to identify unbounded growth, abrupt state changes, or collapse of the autoregressive tra- jectory. The free-tip spectra obtained from the final 60 predicted states are compared with those of the corresponding 60-state, statistically developed CFD reference windows. The blind rollouts remain bounded and preserve regime-dependent periodic behavior over the complete 1000-step interval. The spectra obtained from the final 60 predicted states retain the dominant frequencies of the corresponding CFD reference trajectories. Repeated reuse of the predicted coupled state does not lead to immediate divergence, collapse to a steady solution, or unphysical growth of the plate motion. Because a synchronized CFD reference trajectory is not available over the complete extended interval, these results should be interpreted as evidence of long-horizon numerical and dynamical stability rather than as a direct measure- ment of 1000-step prediction accuracy. 49 0.00 0.25 0.50 0.75 1.00 1.25 Relative difference (a) = 0.05, Deflected-flapping: fluid stability UVP 0.00 0.02 0.04 0.06 Coordinate RMSE (b) = 0.05, Deflected-flapping: plate stability 0.0 0.2 0.4 0.6 0.8 1.0 Normalized FFT magnitude f model = 0.1333 f CFD = 0.1333 A model = 0.239 A CFD = 0.235 (c) = 0.05, Deflected-flapping: spectrum Model, final 60 steps CFD reference window 0.0 0.2 0.4 0.6 0.8 Relative difference (d) = 0.10, Deflected: fluid stability 0.000 0.005 0.010 0.015 0.020 0.025 Coordinate RMSE (e) = 0.10, Deflected: plate stability 0.0 0.2 0.4 0.6 0.8 1.0 Normalized FFT magnitude f model = 0.4333 f CFD = 0.4333 A model = 0.024 A CFD = 0.023 (f) = 0.10, Deflected: spectrum 2004006008001000 Rollout step 0.0 0.2 0.4 0.6 0.8 1.0 Relative difference (g) = 0.30, Flapping: fluid stability 2004006008001000 Rollout step 0.00 0.05 0.10 0.15 0.20 Coordinate RMSE (h) = 0.30, Flapping: plate stability 0.00.10.20.30.40.50.6 Frequency, f 0.0 0.2 0.4 0.6 0.8 1.0 Normalized FFT magnitude f model = 0.1333 f CFD = 0.1333 A model = 0.853 A CFD = 0.851 (i) = 0.30, Flapping: spectrum Figure A.1: Long-horizon stability diagnostics over 1000 fully autoregressive steps for κ = 0.05, 0.10, and 0.30. The left column shows the relative differences between consecutive predicted fluid states, and the middle column shows the coordinate RMSE between consecutive predicted plate states. The right column compares the normalized spectrum obtained from the final 60 predicted states with that of the corresponding 60-state CFD reference window. These quantities are used as stability diagnostics rather than errors against a synchronous 1000-step CFD trajectory. Appendix B. Analytical verification of the differentiable force recon- struction The Level-set quadrature and its geometric gradient are first tested on a circular contour with exact perimeter L = 2πR and derivative dL/dR = 2π. Fig. B.1(a,b) compares the corresponding errors under grid refinement and dif- ferent smoothing-band choices. Both quantities approach their analytical values when the regularized band is adequately resolved. Gradient consistency is fur- ther assessed using R 0 (h) =|J (R + h)− J (R)|, R 1 (h) =|J (R + h)− J (R)− hJ ′ (R)|. (B.1) 50 The fitted small-perturbation slopes are 1.01 and 2.08, consistent with the ex- pected first- and second-order behavior, respectively (Fig. B.1(c)). The vorticity-moment force is additionally tested using a point vortex of strength Γ translating with velocity (c x ,c y ), for which F = Γ(−c y ,c x ), ∂F ∂Γ = (−c y ,c x ).(B.2) The computed force and its automatic-differentiation gradient agree with the analytical expressions to approximately machine precision. Differences between automatic differentiation and centered finite differences are of order 10 −11 –10 −10 (Fig. B.1(d)). 64128256512 Grid size, N 10 5 10 4 10 3 Relative perimeter error (a) Level-set perimeter convergence Fixed = 0.06 = 0.5 x Narrow band: = 1.5x 64128256512 Grid size, N 10 4 10 3 10 2 10 1 Relative error in d L /d R (b) Geometric-gradient convergence 10 4 10 3 10 2 Perturbation magnitude, h 10 7 10 6 10 5 10 4 10 3 10 2 10 1 Taylor remainder Small-h slopes: 1.01, 2.08 (c) Taylor remainder test R 0 (h) R 1 (h) O(h) O(h 2 ) 0.010.020.040.08 Time step, t 10 14 10 13 10 12 10 11 10 10 10 9 Maximum relative error Maximum over N = 128, 256, 512 (d) Moving point-vortex verification Force vs. exact Autograd vs. exact Autograd vs. finite difference Figure B.1: Analytical verification of the Level-set integration and discrete gradients: (a) circular-perimeter convergence; (b) geometric-gradient convergence; (c) Taylor-remainder test; and (d) translating-point-vortex force and gradient verification. 51 Acknowledgments This work was supported by the National Natural Science Foundation of China under Grant Nos. 92470120 and 12672385. References [1] D. Kim, J. Cossé, C. Huertas Cerdeira, M. Gharib, Flapping dynam- ics of an inverted flag, Journal of Fluid Mechanics 736 (2013) R1. doi:10.1017/jfm.2013.555. URL https://doi.org/10.1017/jfm.2013.555 [2] J. E. Sader, J. Cossé, D. Kim, B. Fan, M. Gharib, Large-amplitude flapping of an inverted flag in a uniform steady flow—a vortex-induced vibration, Journal of Fluid Mechanics 793 (2016) 524–555. doi:10.1017/jfm.2016.139. URL https://doi.org/10.1017/jfm.2016.139 [3] A. Goza, T. Colonius, J. E. Sader, Global modes and nonlinear analysis of inverted-flag flapping, Journal of Fluid Mechanics 857 (2018) 312–344. doi:10.1017/jfm.2018.728. URL https://doi.org/10.1017/jfm.2018.728 [4] L. Lu, P. Jin, G. Pang, Z. Zhang, G. E. Karniadakis, Learning nonlinear op- erators via DeepONet based on the universal approximation theorem of op- erators, Nature Machine Intelligence 3 (2021) 218–229. doi:10.1038/s42256- 021-00302-5. URL https://doi.org/10.1038/s42256-021-00302-5 [5] Z. Li, N. B. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. M. Stuart, A. Anandkumar, Fourier neural operator for parametric partial differential equations, in: International Conference on Learning Represen- tations, 2021. URL https://openreview.net/forum?id=c8P9NQVtmnO 52 [6] T. Pfaff, M. Fortunato, A. Sánchez-González, P. W. Battaglia, Learning mesh-based simulation with graph networks, in: International Conference on Learning Representations, 2021. URL https://openreview.net/forum?id=roNqYL0\_XP [7] Z. Li, K. Meidani, A. Barati Farimani, Transformer for partial differential equations’ operator learning, Transactions on Machine Learning Research (2023). URL https://openreview.net/forum?id=EPPqt3uERT [8] Z. Hao, Z. Wang, H. Su, C. Ying, Y. Dong, S. Liu, Z. Cheng, J. Song, J. Zhu, GNOT: A general neural operator transformer for operator learning, in: Proceedings of the 40th International Conference on Machine Learning, Vol. 202 of Proceedings of Machine Learning Research, PMLR, 2023, p. 12556–12569. URL https://proceedings.mlr.press/v202/hao23c.html [9] O. Bublík, V. Heidler, A. Pecka, J. Vimmr, Neural-network-based fluid–structure interaction applied to vortex-induced vibration, Jour- nal of Computational and Applied Mathematics 428 (2023) 115170. doi:10.1016/j.cam.2023.115170. URL https://doi.org/10.1016/j.cam.2023.115170 [10] R. Gao, R. K. Jaiman, Predicting fluid–structure interaction with graph neural networks, Physics of Fluids 36 (1) (2024) 013622. doi:10.1063/5.0182801. URL https://doi.org/10.1063/5.0182801 [11] W. Xiao, T. Gao, K. Liu, J. Duan, M. Zhao, Fourier neural operator based fluid–structure interaction for predicting the vesicle dynamics, Physica D: Nonlinear Phenomena 463 (2024) 134145. doi:10.1016/j.physd.2024.134145. URL https://doi.org/10.1016/j.physd.2024.134145 [12] X. Fan, J.-X. Wang, Differentiable hybrid neural modeling for fluid– structure interaction, Journal of Computational Physics 496 (2024) 112584. 53 doi:10.1016/j.jcp.2023.112584. URL https://doi.org/10.1016/j.jcp.2023.112584 [13] A. Rahman, R. J. George, M. Elleithy, D. Leibovici, Z. Li, B. Bonev, C. White, J. Berner, R. A. Yeh, J. Kossaifi, K. Azizzadenesheli, A. Anandkumar, Pretraining codomain attention neural operators for solv- ing multiphysics PDEs, in: Advances in Neural Information Process- ing Systems, Vol. 37, Curran Associates, Inc., 2024, p. 104035–104064. doi:10.52202/079017-3306. URL https://proceedings.neurips.c/paper\_files/paper/2024/f ile/bc75fa9843a7905bbed9d83895a88f7f-Paper-Conference.pdf [14] J. Brandstetter, D. E. Worrall, M. Welling, Message passing neural PDE solvers, in: International Conference on Learning Representations, 2022. URL https://openreview.net/forum?id=vSix3HPYKSU [15] P. Lippe, B. Veeling, P. Perdikaris, R. Turner, J. Brandstetter, PDE-refiner: Achieving accurate long rollouts with neural PDE solvers, in: Advances in Neural Information Processing Systems, Vol. 36, Curran Associates, Inc., 2023, p. 67398–67433. URL https://proceedings.neurips.c/paper\_files/paper/2023/f ile/d529b943af3dba734f8a7d49efcb6d09-Paper-Conference.pdf [16] Z. Hao, C. Su, S. Liu, J. Berner, C. Ying, H. Su, A. Anandkumar, J. Song, J. Zhu, DPOT: Auto-regressive denoising operator transformer for large- scale PDE pre-training, in: Proceedings of the 41st International Confer- ence on Machine Learning, Vol. 235 of Proceedings of Machine Learning Research, PMLR, 2024, p. 17616–17635. URL https://proceedings.mlr.press/v235/hao24d.html [17] B. List, L.-W. Chen, K. Bali, N. Thuerey, Differentiability in un- rolled training of neural physics simulators on transient dynamics, Com- puter Methods in Applied Mechanics and Engineering 433 (2025) 117441. 54 doi:10.1016/j.cma.2024.117441. URL https://doi.org/10.1016/j.cma.2024.117441 [18] F. Koehler, S. Niedermayr, R. Westermann, N. Thuerey, APEBench: A benchmark for autoregressive neural emulators of PDEs, in: Advances in Neural Information Processing Systems, Vol. 37, Curran Associates, Inc., 2024, p. 120252–120310. doi:10.52202/079017-3822. URL https://proceedings.neurips.c/paper\_files/paper/2024/f ile/d9875ebcf74bccdc5076acab0dbee62c-Paper-Datasets\_and\_Be nchmarks\_Track.pdf [19] J.-Z. Wu, Z.-L. Pan, X.-Y. Lu, Unsteady fluid-dynamic force solely in terms of control-surface integral, Physics of Fluids 17 (9) (2005) 098102. doi:10.1063/1.2055528. URL https://doi.org/10.1063/1.2055528 [20] J.-Z. Wu, X.-Y. Lu, L.-X. Zhuang, Integral force acting on a body due to local flow structures, Journal of Fluid Mechanics 576 (2007) 265–286. doi:10.1017/S0022112006004551. URL https://doi.org/10.1017/S0022112006004551 [21] D. Lentink, Accurate fluid force measurement based on control surface integration, Experiments in Fluids 59 (1) (2018) 22. doi:10.1007/s00348- 017-2464-1. URL https://doi.org/10.1007/s00348-017-2464-1 [22] A.-K. Gao, C. Xie, X.-Y. Lu, Weighted integral methods for fluid force diagnostics in incompressible flows, Journal of Fluid Mechanics 1024 (2025) A57. doi:10.1017/jfm.2025.10854. URL https://doi.org/10.1017/jfm.2025.10854 [23] L.-W. Chen, B. A. Cakal, X. Hu, N. Thuerey, Numerical investigation of minimum drag profiles in laminar flow using deep learning surrogates, Journal of Fluid Mechanics 919 (2021) A34. doi:10.1017/jfm.2021.398. 55 [24] X. Fan, X.-Y. Liu, M. Wang, J.-X. Wang, Diff-FlowFSI: A GPU-optimized differentiable CFD platform for high-fidelity turbulence and FSI simula- tions, Computer Methods in Applied Mechanics and Engineering 448 (2026) 118455. doi:10.1016/j.cma.2025.118455. URL https://doi.org/10.1016/j.cma.2025.118455 56