PaperlayerView exact revision on arXiv

Source-linked AI summary

A full Eulerian finite difference approach for solving fluid-structure coupling problems

Kazuyasu Sugiyama, Satoshi Ii, Shintaro Takeuchi, Shu Takagi, Yoichiro Matsumoto

arXiv:1009.3609v1physics.comp-phphysics.flu-dyn

TL;DR

FSI simulations must couple Eulerian fluid dynamics with Lagrangian structural deformation, while moving or complex geometries can make mesh-based methods costly. This paper develops a full Eulerian finite-difference method using VOF geometry representation and Eulerian transport of a left Cauchy-Green tensor for nonlinear hyperelasticity. The method is validated against analytical and previously available simulation results, and captures nonlinear deformation and shape recovery in the reported tests.

  • Problem

    The central problem is coupling fluid and structural dynamics when conventional formulations use different Eulerian and Lagrangian numerical frameworks.

  • Method

    The method solves fluid-structure equations on a fixed Cartesian grid, represents geometry with VOF, and updates the left Cauchy-Green tensor in the Eulerian frame.

  • Results

    The method is validated against an analytical oscillatory-layer solution and available simulations of lid-driven-cavity and Couette-flow interactions.

  • Takeaways & Limitations

    The fixed-grid formulation can use voxel-based geometry and perform FSI simulations without a mesh-generation procedure once the initial solid volume fraction is supplied.

Abstract

from arXiv · show

A new simulation method for solving fluid-structure coupling problems has been developed. All the basic equations are numerically solved on a fixed Cartesian grid using a finite difference scheme. A volume-of-fluid formulation (Hirt and Nichols (1981, J. Comput. Phys., 39, 201)), which has been widely used for multiphase flow simulations, is applied to describing the multi-component geometry. The temporal change in the solid deformation is described in the Eulerian frame by updating a left Cauchy-Green deformation tensor, which is used to express constitutive equations for nonlinear Mooney-Rivlin materials. In this paper, various verifications and validations of the present full Eulerian method, which solves the fluid and solid motions on a fixed grid, are demonstrated, and the numerical accuracy involved in the fluid-structure coupling problems is examined.

I. INTRODUCTION

FSI couples fluid and structural dynamics across different numerical frameworks, motivating a full Eulerian method that represents both phases on a fixed Cartesian grid. The approach uses volume fraction and one-continuum equations to describe geometry and interactions without relying on a body-fitted moving mesh.

  • Motivation: FSI coupling is challenging because fluid dynamics are conventionally Eulerian while structural dynamics are usually Lagrangian.The difficulty arises from the differing numerical frameworks and the multiphysics nature of the problem.
  • Existing approaches: Body-fitted interface-tracking methods automatically satisfy interface kinematics but require meshes that conform to moving fluid-solid boundaries.ALE and DSD/ST approaches use a shared surface mesh between fluid and solid phases.
  • Existing approaches: Eulerian-Lagrangian methods keep fluid and solid phases on separate fixed Eulerian and Lagrangian grids, respectively.Immersed Boundary, Fictitious Domain, and related methods communicate between the two formulations.
  • Research gap: Lagrangian hyperelastic simulations require temporal displacement updates tied to reference and current configurations, while remeshing can be computationally intensive for complicated geometries or many objects.This motivates describing deformation without repeatedly reconstructing the mesh.
  • Present approach: The proposed method writes governing equations on a fixed Cartesian grid, uses VOF geometry representation, and incorporates nonlinear Mooney-Rivlin constitutive behavior through the Eulerian deformation description.The method targets flexible hyperelastic materials and implementation within standard incompressible fluid solvers.
  • Governing representation: A one-continuum formulation uses one velocity field over the entire domain, making velocity continuity across the fluid-structure interface immediate.The formulation is constructed through volume averaging and uses the solid volume fraction to represent phase occupancy.

B. Constitutive equations for solid

The constitutive formulation expresses incompressible visco-hyperelastic solid stress in the Eulerian frame, using a nonlinear Mooney-Rivlin law based on deformation-tensor invariants. The framework also covers linear Mooney-Rivlin, neo-Hookean, and incompressible Saint Venant-Kirchhoff special cases.

  • Eulerian constitutive framework: The Eulerian formulation writes constitutive equations over the current configuration, unlike finite-element hyperelastic laws commonly written over the reference configuration.This distinction determines how stress is differentiated and incorporated into the momentum equation.
  • Material assumptions: The solid is modeled as an incompressible visco-hyperelastic material undergoing only isochoric motion.The deviatoric Cauchy stress is used for the solid constitutive description.
  • Mooney-Rivlin law: The Mooney-Rivlin formulation retains nonlinear contributions through O(B^2) in the deviatoric Cauchy stress.Its strain-energy representation uses reduced invariants of the right Cauchy-Green tensor.
  • Eulerian closure: The constitutive law can be expressed as a function of the left Cauchy-Green deformation tensor B once coefficients c1, c2, and c3 are specified.This enables the constitutive equation to close entirely in Eulerian variables when B is determined in the Eulerian frame.
  • Special cases: The linear Mooney-Rivlin, neo-Hookean, and incompressible Saint Venant-Kirchhoff models are treated as specific material cases within the formulation.The neo-Hookean model is identified as a particular case of the linear Mooney-Rivlin model.
  • Resulting law: The implemented constitutive law falls within the class of nonlinear Mooney-Rivlin laws.This follows after substituting the deformation-tensor relations into the constitutive expression.

2. Transport of left Cauchy-Green deformation tensor field in the Eulerian frame

The method transports a modified left Cauchy-Green tensor on the fixed Eulerian grid to represent solid deformation while avoiding numerical growth in fluid and interfacial cells. Finite differences, a staggered MAC arrangement, and an exact pressure-Poisson solve complete the fixed-grid procedure.

  • Tensor transport: The left Cauchy-Green tensor B is transported in Eulerian form using the velocity gradient and the upper-convected evolution equation.The initial unstressed solid uses B = I.
  • Stabilization: Updating modified tensor ˜B = φs^αB suppresses deformation information in fluid cells and avoids instability caused by exponential growth of B.Because ˜B vanishes where φs = 0, the fluid domain does not accumulate the same deformation growth.
  • Stabilization: The exponent is restricted to α ≤ 1/2, and the method chooses α = 1/2 while applying a threshold φmin near the interface.The threshold sets ˜B to zero where φs < φmin to balance numerical consistency and stability.
  • Stress coupling: The updated modified tensor ˜B is used to evaluate the phase-averaged solid deviatoric stress together with the constitutive coefficients.The resulting stress expression is evaluated alongside the transport equation.
  • Discretization: All governing quantities are updated by finite differences on a uniform fixed Cartesian grid using a staggered MAC arrangement.Velocity components are placed on cell faces and pressure at cell centers.
  • Pressure solution: The method solves the pressure Poisson equation over the entire domain to satisfy the solenoidal condition, using the fast Fourier transform for accuracy and efficiency.This connects the FSI procedure to standard incompressible-flow algorithms.

B. Time-stepping algorithm

The method advances fluid-structure variables through a three-step projection-based algorithm on a fixed grid, combining explicit and implicit time discretizations with structure-preserving spatial operators.

  • Time-stepping stages: The time-stepping algorithm uses three stages: deformation and volume-fraction update, unprojected velocity and stress calculation, then pressure projection.The first stage updates φ and the modified left Cauchy-Green tensor; the second computes v* and σ̃*; the final stage updates pressure, solenoidal velocity, and stress.
  • Time-stepping stages: The first stage applies a second-order Adams-Bashforth scheme to update the volume fraction and modified left Cauchy-Green deformation tensor.
  • Time-step restriction: The time increment is selected so the CFL number based on advection and elastic-wave speeds is 0.1 or less in all computations.
  • Time-stepping stages: The second stage combines second-order Adams-Bashforth and Crank-Nicolson schemes to calculate the unprojected velocity and stress tensor.
  • Projection stage: The projection stage determines incremental pressure by solving a Poisson equation and updates the solenoidal velocity and stress tensor.
  • Spatial discretization: Spatial derivatives use second-order central differences, fifth-order WENO advection discretization, and an energy-conserving discrete momentum identity.The discretization also chooses interpolation and finite-difference formulas to preserve a stated integral relation in discrete form.
  • Spatial discretization: The spatial treatment accounts for staggered locations of diagonal and off-diagonal deformation components to accurately represent isochoric solid deformation.The derivation uses component-wise deformation-tensor relations, divergence-theorem manipulations, wall kinematics, and the solenoidal condition.

IV. VALIDATION TESTS

Validation uses oscillatory parallel fluid-solid layers as an analytically tractable test, then compares profiles and convergence against sharp-interface solutions. The method captures nonlinear constitutive behavior, while coupling errors show near-first-order convergence.

  • Validation tests: Three validation tests examine an analytically solvable oscillatory layered problem, deformable-solid motion in a lid-driven cavity, and particle interaction in Couette flow.A third validation also checks reversibility of hyperelastic shape under external shear strain.
  • Oscillatory parallel layers: The parallel-layer benchmark contains fluid-solid-fluid layers between two oscillating walls, with comparisons made against sharp-interface solutions.The upper and lower walls move with opposite sinusoidal velocities, and the numerical simulation is periodic in the x direction.
  • Material response: The simulations compare neo-Hookean and Saint Venant-Kirchhoff materials, whose constitutive responses differ between linear and nonlinear cases.The Saint Venant-Kirchhoff formulation has c3 ≠ 0 and is nonlinear in displacement, whereas the linear Mooney-Rivlin case has c3 = 0.
  • Material response: Increasing spatial resolution makes the computed velocity profiles converge toward sharp-interface solutions and captures the nonlinear solid response.For the Saint Venant-Kirchhoff case, the profile approaches the nonlinear reference solution as the number of grid points increases.
  • Grid convergence: L2 and L∞ velocity errors are nearly proportional to N_y^-1, indicating near-first-order accuracy in fluid-structure coupling.The L2 norm is treated as a global indicator and the L∞ norm as a local maximum indicator.
  • Grid convergence: The near-first-order convergence is attributed to locally first-order mixture-stress accuracy at the interface, which dominates global accuracy.The shear-stress grid-convergence study likewise confirms first-order interface accuracy with respect to grid size.
  • Parameter sensitivity: Wall-friction sensitivity is evaluated across linear Mooney-Rivlin, neo-Hookean, and Saint Venant-Kirchhoff models as a function of 2(c1 + c2).The wall-friction root-mean-square amplitude is sampled over the final cycle after a 20-cycle computation.

B. Comparison with independently conducted FSI analyses

The full Eulerian method is compared with independently conducted FSI analyses and examined for grid convergence in particle motion and deformation. The comparisons show agreement with reference results, while convergence studies indicate first-order accuracy.

  • Comparison with Zhao et al.: The present full Eulerian simulation agrees closely with Zhao et al.’s Lagrangian-tracking result for highly deformed particle shapes in lid-driven flow.The comparison uses a 1024 × 1024 mesh, and no special particle-wall overlap-avoidance artifact is implemented.
  • Grid convergence: Centroid trajectories converge toward the highest-resolution curve as the grid is refined from 64 × 64 to 1024 × 1024.The trajectories are evaluated over t ∈[0, 20].
  • Grid convergence: Both centroid errors are nearly proportional to N_x^-1, indicating first-order accuracy in the fluid-structure coupling.The errors are measured in L2 and L∞ norms relative to the highest-resolution result.
  • Modal deformation: The n = 2 deformation mode reaches its largest elongation at about t = 5, while modal profiles settle to convergent curves with grid refinement.The synchronized increases across modes provide a deformation-based verification of the method.
  • Modal deformation: Modal-amplitude errors are nearly proportional to N_x^-1, indicating first-order accuracy for capturing particle deformation.The errors are reported in both L2 and L∞ norms for modes n = 0, 1, 2, 3, 4 and 5.

2. Two particles interaction in a Couette flow

The full Eulerian method is tested on interactions between two deformable particles in Couette flow against a body-fit Lagrangian-mesh analysis. It reproduces the reported interaction modes and approaches the reference centroid motion as resolution increases.

  • Comparison with Gao & Hu: The present simulation reproduces the complicated “roll over” and “bounce back” particle interactions reported by Gao & Hu.Particle shapes are compared at five time instants using a 1024 × 512 mesh.
  • Comparison with Gao & Hu: The simulated particle shapes agree with Gao & Hu’s well-validated body-fit Lagrangian-mesh result.The dotted markers visualize deformation but are not used to compute solid stress or strain.
  • Centroid motion: Higher spatial resolution brings the present centroid trajectories closer to the full Lagrangian result, especially when the particle gap is narrow.Peaks occur around t = 3.0, t = 16.0 and t = 20.0, corresponding to narrow particle gaps.

C. Reversibility in shape of hyperelastic material

The method is examined for reversibility when a hyperelastic particle is released from shear. Refinement improves shape recovery, while energy-budget diagnostics show that fluid-solid energy exchange is well conserved.

  • Reversibility setup: The fixed Eulerian grid is tested for recovering a circular particle after shear is released by stopping both plates at t = 4.The particle is initially circular with radius 0.75 and is simulated between plates at y = 1 and y = −1.
  • Reversibility dynamics: After release, the particle shape recovers while transverse elastic waves propagate inside the solid and reflect at the fluid-structure interface.The waves are damped through repeated reflections, and the particle exhibits tank-treading-like motion before release.
  • Energy conservation: The kinetic-energy budget closes with an absolute numerical error below 10^-5, much smaller than the individual energy-term variations.The result indicates that the system is well conserved during the simulation.
  • Grid convergence: Increasing spatial resolution produces convergent interface shapes at t = 4 and better recovery of the initial curve at t = 200.The interfaces are identified by the isoline φ_s = 1/2 for both linear Mooney and incompressible Saint Venant-Kirchhoff materials.
  • Grid convergence: Residual higher-order deformation amplitudes decrease nearly as N_x^-1, indicating first-order accuracy in recovering the unstressed shape.The residual modes examined are |R_n| for n = 2, 4 and 6 at t = 200.

2. Shape reversibility of a rectangular particle

The rectangular particle undergoes large deformation during imposed shear and gradually recovers its unstressed shape after release. Grid refinement improves shape reversibility, while energy exchange between fluid and solid remains reasonably accurate.

  • Shape recovery: The rectangular particle gradually recovers its unstressed shape after the shearing force is released at t = 4.At t = 6, the deformed particle is shown recovering toward its initial shape.
  • Energy exchange: The kinetic-energy budget indicates that numerical error is much smaller than the variation of individual contributions.This indicates that energy exchange between fluid and solid phases through solid deformation is reasonably guaranteed.
  • Material-point reversibility: Tracer distributions demonstrate recovery of material points toward the initial configuration after the particle deforms at t = 4.The strongly stretched left-top and right-bottom corners gradually resume their original shapes.
  • Grid convergence: Increasing the number of grid points makes the particle shapes at t = 4 and t = 200 converge and improves shape reversibility.The rectangle edges remain smeared because numerical dissipation affects the fifth-order WENO update of φs.
  • Overall accuracy: The method captures shape reversibility when grid resolution is sufficiently high and has first-order accuracy with respect to grid size in fluid-structure coupling.The method was also validated against analytical and available simulation results for several fluid-structure interaction problems.
  • Accuracy limitations: Sharper interface capture and more robust time advancement are identified as ongoing requirements for improving accuracy toward practical applications.The present method smooths strain-rate discontinuities across the fluid-structure interface at grid scale.

Appendix A: Finite difference descriptions

The appendix introduces finite-difference operators used to describe the discretized equations.

  • The finite-difference formulation begins by introducing operators used to discretize the governing equations.

2. For momentum conservation equation (26)

The momentum equation discretization is organized term by term, with interpolation and tensor-related operations specified for the Cartesian grid.

  • Interpolation operators are introduced for discretizing the momentum equation and related quantities.
  • The discretization separately addresses advection terms, pressure gradients, and divergences of deviatoric stress tensors.
  • The trace of the modified deformation tensor is written under the condition Bzz = 1.

3. For the advection terms in (13) and (22)

The advection terms for the solid volume fraction and deformation tensor are evaluated with fifth-order WENO discretization and interpolated velocities on the staggered grid.

  • A fifth-order WENO scheme is applied to advection terms for φs, B̃xx, and B̃yy defined at cell centroids.
  • The advection term vx∂xq is discretized for centroid-defined quantities, with a positive ε preventing division by zero.
  • For B̃xy defined at a cell apex, the advection term v · ∇B̃xy is evaluated using interpolated velocities.

Appendix B: Spectral algorithm to find sharp interface solution for the parallel layers problem

The sharp-interface solution uses Fourier expansions on the symmetric upper-half domain, enforcing interface, wall, and symmetry conditions. Modal equations are truncated for fast transforms, with separate linear and nonlinear time-integration treatments.

  • Sharp interface formulation: The method exploits symmetry about y = 0 by solving in the upper half-region and expanding fluid and solid velocities in Fourier series.The expansions are constructed to satisfy velocity continuity at the interface, no slip at the upper wall, and symmetry in the solid.
  • Sharp interface formulation: The modal formulation derives solid displacement, nonlinear solid stress, momentum relations, and shear-stress continuity at the fluid-structure interface.The nonlinear stress contribution is related to its expansion coefficients before applying sine-function orthogonality.
  • Spectral discretization: The infinite modal summation is truncated at k = K − 1, and powers-of-two K enable fast Fourier sine transforms for evaluating fluid velocity and solid displacement.The spectral algorithm solves for vf,k, us,k, VI, and UI using the truncated system.
  • Time integration: For the linear Mooney-Rivlin case c3 = 0, σNL vanishes and harmonic separation of variables reduces the time-dependent differential equations to algebraic equations.The velocity and displacement fields are represented with complex harmonic amplitudes.
  • Time integration: For c3 ≠ 0, the nonlinear system requires numerical time integration using second-order Adams-Bashforth and Crank-Nicolson schemes, producing temporally periodic solutions after sufficient computation.The algorithm updates UI, vf,k, us,k, and VI at the next time level from the two preceding time levels and prescribed wall velocity.
  • Validation: The truncated-mode convergence check identifies K = 2048 as sufficiently accurate to serve as the reference solution within the parameter range shown in figure 9.This reference is used for comparison and validation of the full Eulerian model.
Loading 1009.3609v1…