Open-source simulation of strongly-coupled fluid-structure interaction between non-conformal interfaces

Design code-based “life-safety” requirements for structural earthquake and tsunami design offer reasonable guidelines to construct buildings that will remain standing during a tsunami or seismic event. Much less consideration has been given to assessing structural resilience during sequential earthquake and tsunami multi-hazard events. Such events present a series of extreme loading scenarios, where damage sustained during the earthquake influences structural performance during the subsequent inundation. Similar difficulties exist with respect to damage sustained during tropical events, as wind and fluid loading may vary with structural response or accumulated damage. To help ensure critical structures meet a “life-safety” level of performance during such multi-hazard events, analysis software capable of simulating simultaneous structural and fluid dynamics must be developed. To address this gap in understanding of non-linear fluid-structure-interaction (FSI), an open-source tool (FOAMySees) was developed for simulation of tsunami and wave impact analysis of post-earthquake non-linear structural response of buildings. The tool is comprised of the Open-source Field Operation And Manipulation software package and OpenSeesPy, a Python 3 interpreter of OpenSees. The programs are coupled via preCICE, a coupling library for partitioned multi-physics simulation. FOAMySees has been written to work in a Linux OS environment with HPC clusters in mind. The FOAMySees program offers a partitioned conventional-serial-staggered coupling scheme, with optional implicit iteration techniques to ensure a strongly-coupled two-way FSI solution. While FOAMySees was developed specifically for tsunami-resilience analysis, it may be utilized for other FSI applications with ease. With this coupled Computational Fluid Dynamics (CFD) and Finite Element Analysis (FEA) program, tsunami and earthquake simulations may be run sequentially or simultaneously, allowing for the evaluation of non-linear structural response to multi-hazard excitation.


Introduction
Design code-based "life-safety" requirements for structural earthquake and tsunami design offer reasonable guidelines to construct buildings that will remain standing during a tsunami or seismic event. Much less consideration has been given to assessing structural resilience during sequential earthquake and tsunami multi-hazard events. Such events present a series of extreme loading scenarios, where damage sustained during the earthquake influences structural performance during the subsequent inundation. Similar difficulties exist with respect to damage sustained during tropical events, as wind and fluid loading may vary with structural response or accumulated damage. To help ensure critical structures meet a "life-safety" level of performance during such multi-hazard events, analysis software capable of simulating simultaneous structural and fluid dynamics must be developed. Additionally, such a program must not only account for initial earthquake damage, but also the interplay of fluid-induced forces and the motion of structures on which they act. This integrated approach of analysis of cascading hazards will allow engineers to better represent the physical conditions and expected conditions, enabling more resilient designs for critical structures.
To address this gap in understanding of non-linear fluidstructure-interaction (FSI), an open-source tool (FOAMySees) was developed for simulation of tsunami and wave impact analysis of post-earthquake non-linear structural response of buildings. The tool is comprised of the Open-source Field Operation And Manipulation (OpenFOAM, Weller et al., 1998) software package and OpenSeesPy (Zhu et al., 2018), a Python 3 interpreter of OpenSees (McKenna et al., 2010). The programs are coupled via preCICE (Bungartz et al., 2016), a coupling library for partitioned multi-physics simulation. FOAMySees has been written to work in a Linux OS environment with HPC clusters in mind. The FOAMySees program offers a partitioned conventional-serialstaggered coupling scheme, with optional implicit iteration techniques to ensure a strongly-coupled two-way FSI solution. While FOAMySees was developed specifically for tsunamiresilience analysis, it may be utilized for other FSI applications with ease. With this coupled Computational Fluid Dynamics (CFD) and Finite Element Analysis (FEA) program, tsunami and earthquake simulations may be run sequentially or simultaneously, allowing for the evaluation of non-linear structural response to multi-hazard excitation.
1.1 Research motivation 1.1.1 Numerical methods for fluid-structure interaction simulation Several commercial packages with non-linear FSI capabilities are currently available, with a long-standing history of validation and community support (see LS-Dyna, STAR-CCM+, ABAQUS, COMSOL, etc.). Though these programs have been shown to be capable of solving complicated multiphysics problems with various coupled-numerical methodologies, the large barriers to entry associated with these commercial packages make these programs undesirable for research purposes. These barriers to entry include costs of licensing of the commercial codes and overhead for preprocessing and post-processing software packages, inflexibility of the sub-routines making solution of specialized problems complicated, as well as lack of transparency of the source codes. While suitable for solution of a broad range of fluid-structure interaction problems and have been used extensively in numerical and experimental studies of civil engineering structures, these programs each have drawbacks and none of the programs listed have been specifically developed for the purposes of civil/structural engineering, making earthquake analysis a challenge with many of these codes. A list of commonly-utilized parallelized software programs for Civil/Structural engineering along with their capabilities is shown in Table 1.
Methods currently exist to couple open-source finite element codes (Deal-II, CalculiX, FENiCs, etc.) with computational fluid dynamics codes (see preCICE, associated adapters), but none of the existing commercial structural mechanics solvers that have been coupled with a fluid dynamics solver allow for representation of three-dimensional geometries within the fluid domain with onedimensional, 6 degree of freedom structural beam elements. Such methods are currently in development which utilize an unstructured polyhedral finite volume method for solution of the structural deformation of modelled structures within Solids4Foam, an open-source fluid-structure interaction toolkit available for OpenFOAM. Despite the finite volume method offering high accuracy in solution of structural mechanics problems, this is a developing field of research and the finite volume method (unlike the finite element method) has not been extensively validated for earthquake engineering purposes. While strongly-coupled software programs capable of performing non-linear fluid-structure interaction of deformable bodies along with turbulence modelling within the fluid domain exist currently, few are open-source, and of those that are open-source, all utilize structural analysis software programs, which are not commonly utilized for the purposes of nonlinear time history analysis of civil engineering applications and/or are geared toward conjugate heat transfer analysis rather than nonlinear fluid-structure interaction.

Impetus for research
It is assumed that for very small deformations without periodicity a CFD model with non-slip velocity boundary conditions will suffice in determining applied force from impinging fluids, but for cases where displacements involve resonant vibration of a structure or where a structure possesses insufficient stiffness in particular degrees of freedom to resist incoming forces, the change in structural position and velocity within a flow can have a much more noticeable effect on fluid behavior and forces. In the case of larger deflections, or in the case of fluid-excited structural vibrations, fluid forces and their direction of application can result in changes of the pitch or position of the Frontiers in Built Environment frontiersin.org surfaces upon which they impinge, possibly leading to changes in loading. For tall, slender structures or structures with long spans such as bridges or highways, small changes in rotation at connections can result in large displacements of components of the structure. For example, a wide bridge with a long span could experience twisting of the bridge deck along the axis of its span, changing the position of the leading edge of the bridge deck structure within a flow transverse to the bridge span and the angle of inclination of the deck with respect to the flow. Furthermore, frequency based interactions between the structure and loading from surrounding flows can result in resonance between the structure and the vortexes shed in the wake of the structure, leading to vortex induced vibration. In order to determine actual forces and impulses imparted from tsunami inundation events upon structures post-earthquake, expanding the extents of numerical analyses to realistic length scales of design structures subjected to tsunamis is necessary to accurately simulate structural deformations due to wave impingement and periodic fluid forcing from post-impact inundation flow. Additionally, since the feasibility and costs of experimental methods are prohibitive at full scale, utilizing numerical methods such as CFD for the determination of tsunami demands on structures at full scales presents the most economical alternative as well. By validating numerical models of experiments conducted at a reduced scale, validated numerical models may be used to simulate events at full scale in order to assess the accuracy of design equations and make recommendations for improvements. Thus identifying compatible combinations of CFD and FEM software is crucial to realizing full-scale simulations of consecutive earthquake and tsunami multi-hazard events.
Bearing these considerations in mind, OpenSees and OpenFOAM were selected for this research as the FEM and CFD simulation software, respectively. OpenSees is commonly used for non-linear FEM analysis of structures, particularly for earthquake simulation and hysteretic analysis of structures. OpenSees has a long-standing history in the civil engineering field and has been validated against experimental results countless times, providing reasonable confidence in the accuracy of the program when utilizing it for non-linear structural mechanics and dynamics. Similarly, OpenFOAM is used within the civil engineering community for CFD simulation of wind and water hazards, and has been validated against experiments involving hydrodynamic impact with great accuracy (Douglas and Nistor, 2014;Douglas et al., 2015;Motley et al., 2016;Wong, 2015;Sarjamee et al., 2017a;Sarjamee et al., 2017b;Qin et al., 2018;Qin, 2019;Winter 2019;Croquer et al., 2022;Elsheikh et al., 2022a;Elsheikh et al., 2022b;Liu et al., 2022).
OpenFOAM and OpenSees are both free, scalable, and plenty of example cases exist for both codes, which allows engineers to utilize the tools handily and learn through trial and error. The proposed application programming interface that couples OpenFOAM and OpenSees will assist engineers in analysis of structures subject to cascading hazards by offering a methodology for determining earthquake and tsunami forces at realistic scales while accounting for accumulated damage and non-linear fluid-structure interaction; this is particularly of importance when testing prototypical structures is not feasible or when experimental results from small-scale tests cannot be extrapolated accurately beyond their test-scales. Thus, the intention of this research is to expand the current capabilities of these open-source modelling methodologies such that they may be utilized to simulate full-scale damaged-state structural responses to non-linear fluid loading, and to validate the program against analytical and experimental test cases.

Scope of research
FOAMySees, an application programming interface (API) to allow for strongly-coupled FSI between OpenFOAM and OpenSeesPy, was developed to address the gaps in simulation capability described in the previous section. FOAMySees shows strong correlation with other FSI simulation software when benchmarked against analytical test cases. FOAMySees allows for any element formulation within OpenSeesPy to be utilized and coupled with a CFD simulation for the purposes of two-way strong fluid-structure interaction. Furthermore, the capabilities of FOAMySees allow for either preliminary or concurrent simulation of seismic excitation of structures with the OpenFOAM simulation, opening the possibility of simulation of non-linear structural response to multi-hazard events. The API is written such that users may build and implement any OpenSeesPy model, utilizing it for FSI analysis, before, during, or after additional forcing time histories as defined by the user. The program capabilities are intended to be of use to the civil engineering community to assist in the understanding of multi-hazard analysis and design of resilient structures. In particular, the capabilities of the program support high-fidelity CFD analyses of tsunami bore impact and inundation loading applied to non-linear post-yield multiple degree of freedom finite element models, allowing for the investigation of non-linear structural system responses to extreme hydrodynamic loading.

Computational mechanics methodologies 1.2.1 Computational fluid dynamics-OpenFOAM
Fluid modelling was completed using computational fluid dynamics (CFD) with the Open-source Field Operation And Manipulation (OpenFOAM) (Weller et al., 1998) software package, namely, OpenFOAM Foundation's OpenFOAM.org-Version 8. OpenFOAM is a collection of C++ libraries that may be compiled to create individual applications, which are broadly categorized as either solvers or utilities. Cases examined here were simulated using the olaFlow solvers (Higuera, 2018), an open-source project developed within the OpenFOAM framework which solves the three-dimensional Volume Averaged Reynolds Averaged Navier Stokes equations (VARANS) using the finite volume discretization, allowing for the simulation of physically correct two-phase incompressible fluids. Additionally, olaFlow provides wave generation and active absorption functionality, which is useful for coastal engineering applications including tsunami-like wave generation. In cases when pure fluid phases are modeled where porosity is neglected, the VARANS equations utilized within olaFlow are reduced to the classical RANS equations utilized in the standard OpenFOAM application interFoam, from which the solvers for olaFlow were derived. Since porosity effects were not included in this study, similar CFD modelling procedures to those conducted by Motley et al. (2016) were performed for the following analytical correlation studies. OpenFOAM patch boundary conditions utilized in this study are listed in Table 2.

Modelling of multi-phase air-water mixture
The two incompressible phases (water and air) are tracked using the Volume of Fluid (VOF) technique to represent complex free surface configurations Hirt and Nichols (1981). An indicator function α is defined for the volume fraction of the two-phase fluid, which has a value of 1.0 corresponding to regions occupied by one phase, in this case water (ρ = 1, 000 kg/m 3 , ] = 1.0 × 10 −6 m/s 2 ), and a value of 0.0 for the other, in this case air (ρ = 1.22 kg/m 3 , ] = 1.48 × 10 −5 m/s 2 ), where ρ = mass density of the fluid; and ] = kinematic viscosity. Intermediate values indicate cells contain a mixture of water and air, where the freesurface between the fluids is not resolved explicitly. Where free-surface tracking was necessary in this study, a volume fraction of 0.5 was used to identify an approximate free-surface. The cell's fluid phase fraction is represented by a scalar field with the variable, α. Furthermore, the VOF method assumes the fluid phases are immiscible, where each phase remains largely separate from one another, which is reflected by the fact that a single set of continuity and momentum governing equations are solved for both fluids, as opposed to Eulerian Multiphase (EMP) methods, which fully resolve phase interactions by incorporating interaction models and solving separate sets of continuity and momentum equations for each fluid phase.

Mesh motion in openFOAM
Motion of FSI interfaces for the present study is handled within OpenFOAM by the displacementLaplacianFvMotionSolver Class, which calculates the near field cell displacements required to satisfy a Laplacian diffusivity scheme specified by the user. This allows for smooth motion and deformation of the entire mesh during movement of boundary patches. Many other options for solution of cell deformation exist within OpenFOAM; however, the current coupling methodology only utilizes displacementLaplacian based mesh motion. As such, a discussion of alternative mesh motion techniques within OpenFOAM are omitted here for brevity. The displacementLaplacian finite volume motion solver only requires the displacement field at the boundary patch be specified either in memory or within the case files. This is handled dynamically in memory by the proposed coupling method, allowing real-time updates of the patch deformation during coupled fluid-structural analysis. Explanation of the displacement Laplacian methodology for finite volume method mesh motion is available in the original paper by Jasak and Tukovic (2006). For more detailed information on dynamic mesh capabilities within OpenFOAM, see Jasak (2009).

Computational solid mechanics-openSeesPy
OpenSeesPy (Zhu et al., 2018 Due to the flexibility of OpenSees, sequential analysis of earthquake loading and tsunami loading to structures may be implemented readily. For implementation of a fluid-structure interaction analysis, particle finite element methods (PFEM) have been developed for OpenSees; however, PFEM in OpenSees is designed to solve fluid flow problems without considering turbulence modeling. Fluidstructure interaction within the present study is handled within OpenSees with an additional solution loop after initial analysis is complete, obtaining fluid forces from an OpenFOAM simulation running in parallel to the OpenSees simulation. This loop functions as a listener from the coupling driver, controlling the communication of nodal forces and displacements to and from the coupling interface software. Python module commands not included within the OpenSeesPy library are utilized for this loop, and are listed and explained in detail in Section "Usage of preCICE Python Language Bindings."

Coupling methodology-preCICE
preCICE (Bungartz et al., 2016) is a library developed for the purpose of coupling multiphysics simulations in a partitioned method, utilizing iteration to ensure interface acceleration convergence between coupled surfaces and maintain stability. The approach allows for recovery of the coupling solution as if the problem were to be solved wihtin a single program through use of a partitioned "black box" approach and running coupled applications in parallel. The library functions by transferring data arrays between coupling partners for serial or parallel computation of fully-coupled solutions between codes. Convergence measures, data transfer, timestep iteration, and surface-to-surface mapping between coupling partners are managed by the preCICE library. These surface-to-surface mappings are most generally handled with radial basis function (RBF) mesh interpolation from the source mesh to a destination mesh.

preCICE adapters
preCICE is structured in a modular manner, with applications interfacing with the preCICE library via API function calls, which allow access to the routines of the library. Applications may be "plugged-in" to the preCICE library coupling routines by means of a coupling "adapter." These adapters interface with the coupled program, taking data from the program and passing it to the preCICE library routines to be utilized in a manner specified by the user through configuration files. preCICE offers many adapters for open-source finite element and finite volume codes; however, many of these codes have not been validated for civil/structural engineering applications, and are not optimized to allow for a structural engineer to pick up and readily utilize them for engineering design purposes.
The preCICE API function calls are available in many languages via bindings. These bindings allow for usage of the preCICE library within code written in the language shared by the bindings, and offer flexibility for users to implement their own scripts and code for dynamically or temporally assigning boundary conditions within a desired numerical simulation.

Usage of preCICE python language bindings
Calls to the preCICE API are managed through a precice.Interface() object. This object has various classes that control data communication and solver progression. The classes of precice.Interface utilized in FOAMySees include: The latter two call upon the following preCICE action objects to return Boolean values: action_write_initial_data () action_write_iteration_checkpoint () action_read_iteration_checkpoint () A simplified version of initialization routines for FOAMySees along with their associated preCICE python language binding commands are shown in Figure 1A. The precice.Interface class and sub-classes are utilized within the main solution loop of FOAMySees to control the progression of the solution and manage iteration of OpenSeesPy until preCICE determines the coupling residual tolerance has been met. A simplified version of the main solution loop is shown schematically in Figure 1B. Calls to preCICE API functions are highlighted in orange, whereas the functions themselves are shown in blue.

A new framework for open-source fluidstructure interaction-FOAMySees
The proposed program, "FOAMySees", a portmanteau of "OpenFOAM' and "OpenSees," may be qualified as a coupling adapter for OpenSeesPy finite element models, structured within Python using preCICE's Python language API bindings. While it is specialized to work with OpenFOAM, other CFD and FEM solvers could be coupled with OpenSeesPy via FOAMySees. In addition to this, the possibility exists within preCICE to couple multiple model models of similar types together. This would allow multiple OpenSeesPy models to be coupled to a single OpenFOAM model, allowing for simultaneous FSI of ensembles of structures. It would also allow users to couple another solver to the OpenFOAM/OpenSeesPy coupling, such as a shallow water solver to apply CFD inlet boundary conditions.

Theoretical background
The solution of complex turbulent fluid dynamics problems such as a tsunami inundation must be completed in three Frontiers in Built Environment frontiersin.org dimensions, due to the three-dimensional nature of turbulence. Since many structural finite element models utilize onedimensional center-line and spring elements to represent three-dimensional structural components, additional steps must be taken to convert data such as rotations of structural nodes to displacements that accurately represent the rotation of the three-dimensional structural element modelled by elements along its center-line or mid-plane.
To account for this, volume must be generated in the coupled modelling space by extruding cross sectional shapes along the center-line elements. In FOAMySees, the cross sectional shape of the center-line elements is constructed by defining "data communication branch nodes" which exist solely outside the fluid solution and structural solution. "Data communication branch nodes" refer to the additional non-structural nodes which are utilized by FOAMySees and inaccessible to both solid and fluid solutions. These data communication branch nodes are radially centered relative to the FEM nodes and rotate rigidly about them.
Displacements of the finite element model are interpolated to these data communication branch nodes by the FEM node/branch node relationship described in Eq. 1 for each FEM-branch node group. The branch displacements from the FEM model at the completion of the current coupling timestep (t n ), are utilized for progression of the finite volume domain displacement and evolution of the fluid model during the next coupling timestep (t n+1 ). In Eq. 1, FIGURE 1 (A) preCICE commands utilized within FOAMySees initialization, (B) Diagram of preCICE subroutine calls within FOAMySees. This diagram is solely intended for informing about preCICE API function calls and where these calls are necessary, and is not indicative of the full coupling procedures for FOAMySees. For specific location of preCICE API function calls within FOAMySees, refer to the FOAMySees source code.

Frontiers in Built Environment frontiersin.org
BranchRotation refers to rigid-body rotation of branch nodes about FEM nodes using rotation matrices.
x i,branchj t n+1 The mesh created by the collection of data communication branch nodes for all FEM nodes, the data communication mesh, is coupled by preCICE to the fluid mesh patch surfaces which represent the structure within the CFD model. Forces from the fluid model are conservatively applied to these branch nodes using radial basis function mapping from the CFD mesh patch nodes to the branch nodes. For all FEM nodes, the force vectors of the branch nodes associated with a given FEM node are summed to construct the force vectorF i (t n ) at each timestep. More specifically, for a standard FOAMySees model FEM node with N branches, the force relationship between an FEM node and its branch nodes at any particular time is described in Eq. 2. Summarily, each FEM nodes' forces are a combination of the forces applied to all of its branches.
Equilibrium at the i th node of a finite element model is defined as the balance of structural resistance and kinetic motion forces with external applied forces. For equilibrium at an interface between a structural dynamics model and an external dynamic force field to be satisfied, input and output work for a period of time must be equal or within satisfaction of a tolerance for residual in force between two sources. This means externally applied work for each node i must be roughly balanced by the product of the displacement of the structure with the resistance of the motion of the structure, as given in Eq. 3. For a single node within the structural dynamics model, acceleration in six degrees of freedom, € x, is integrated in time for the duration of a single timestep and converted to kinetic motion in the form of velocities, _ x, and displacements, resulting in shape changes to the undeformed structural configuration,x. The forces associated with each of these components of motion -inertial, damping, and internal structural forces -are calculated by finding the product of these components of motion with their respective coefficients within the characteristic equation of motion, mass (m i ), damping (c i ), and tangent stiffness (k i ), respectively. The effective resisting force of a given node is determined by summing these three contributions to internal force together as specified in Eq. 3.
Unbalanced forces as defined by Eq. 4 are calculated by subtracting the internal dynamic resistance force components of the structural motion and deformed configuration at the current time from the applied forces, denoted asF i (t n ) in Eq. 2.
These unbalanced forces are used at the end of the fluid model's progression to the next coupling checkpoint at t n+1 to unbalanced forces between the beginning (t n ) and end (t n+1 ) of the current coupling timestep. This relationship is summarized in Eq. 5. The unbalanced forces,F i,APPLIED (t n ), between the current iteration of the progressed fluid model state (at t = t n+1 ) and the current iteration of the un-progressed structural model state (at t = t n ) are applied to the structural dynamics model to be resolved into kinetic motion, damping resistance, and elastic and inelastic structural resistance over the FEM model's progression from the beginning (t n ) to the end (t n+1 ) of the current coupling timestep.
Depending on the coupling scheme, the solution process for each solver is completed serially or in parallel, and implicit or explicit; as such, these processes may occur in a different order than listed here, or may occur several times prior to progressing to the next timestep.

FOAMySees solution procedure
Prior to entering the solution procedure loop, FOAMySees initializes the OpenSeesPy model and FEM/branch node interpolator and sends FOAMySees data communication mesh (es) and initial values to preCICE for coupling to OpenFOAM. The following operations are then performed each timestep in sequence by FOAMySees components preCICE, OpenFOAM, and OpenSeesPy.

OpenFOAM: Reads data communication mesh displacements
and integrates fluid dynamics model over time with mesh motion 2. preCICE: Calculates branch nodal forces at end of coupling timestep, and reads displacements and forces from both OpenFOAM and OpenSeesPy at end of OpenFOAM timestep to ensure residual convergence for implicit coupling schemes. 3. FOAMySees: Reads forces from preCICE, sums branch node forces from CFD analysis to FE model node locations, and calculates net nodal forces for OpenSeesPy timestep by subtracting nodal unbalanced forces and applies net nodal forces to OpenSeesPy nodes 4. OpenSeesPy: Integrates structural dynamics model over time and returns nodal displacements and unbalanced forces to FOAMySees 5. FOAMySees: Projects displacements from each FEM node to its associated branch nodes via rotation matrices 6. preCICE: Passes FOAMySees calculated data communication mesh displacement array to OpenFOAM for application of displacement boundary conditions during the next coupling timestep A full description of the participants in the coupling process is shown in Figure 2. Software participants are listed along the left hand side of the figure. The order of operations for both FOAMySees initialization and timestepping is diagrammatically shown as an initialization phase, two interior loops, one solution loop for fluid solution internal residuals and one for those of the structural solution internal residuals, and one exterior loop per timestep, ensuring satisfaction of tolerances for residuals of coupling variables between participants.

Surface-to-structure coupling using preCICE and FOAMySees
After initialization of the structural model along with optional branch nodes, the data communication mesh is passed to preCICE

Structural modeling and branch interpolation
During each timestep, after forces are mapped by preCICE via RBF to the FOAMySees branch nodes associated with each FEM structural node, branch node forces are summed and applied to the FEM nodes for evolution of the structural model during the next timestep or coupling iteration. The displacements of the FEM nodes are calculated using a user-specified numerical integration scheme at a maximum specified timestep, for a given number of sub-cycles. The branch node displacements are passed to preCICE at the end of the converged timestep. The displacements from FEM nodes are applied to the branch nodes associated with each FEM node via translation such that each of the n branch nodes per FEM node moves rigidly with the FEM node. The process of rigid-body rotation of branches about their FEM nodes is shown schematically in Figure 6. The branch nodes are then rotated rigidly about the

Frontiers in Built Environment
frontiersin.org FEM nodes via rotation matrix transformation, according to each nodes' calculated local slope. Using the locations of these branch nodes at the previous timestep with respect to the deformed position of the nodes along the surface patch within the OpenFOAM model, preCICE interpolates the displacements from the branch nodes to the CFD patch surface nodes using the RBF mapping used previously. This process is shown schematically in Figure 7.

Time integration and time stepping
Time integration of the fluid model is completed with an Euler scheme. Stability of the fluid model requires a timestep small enough for satisfaction of Courant-Frederichs-Lewy cell fluid velocity criteria. The Newmark method with coefficients of γ and β valued 0.5 and 0.25 respectively was utilized for time integration of the structural model. As such, the maximum timestep increment for the structural dynamics simulation is contingent on structural material wave propagation timestep requirements and minimum element eigenvalue. A timestep sufficiently small enough to meet all stability criteria for the numerical integration routines of both structural and fluid solvers is chosen as the timestep for both the exchange of coupling data and individual solver time progression.

Coupling scheme
An explicit time progression scheme for coupling data is used in the following validation cases. Sub-cycling for CFD and FEM solvers is possible; alternatively, iterative techniques for interface manifold acceleration convergence such as Aitken Under-Relaxation and Iterative-Quasi-Newton Inverse-Least-Squares (IQN-ILS, Degroote and Vierendeels, 2011) may be implemented. These techniques are used to reduce analysis duration by increasing timestep sizes and improve stability of the coupling and are recommended over solver sub-cycling. In general, very small timesteps are recommended for first-order time integration schemes, as numerical damping may be large for explicitly coupled simulations with large timesteps.

Computational cost
The relative times for each computational component will vary depending on the mesh size of each participant and the boundary/ initial conditions of each model. For a small scale structure modelled in FOAMySees with realistic materials (concrete, steel) and small displacements, comprised of 300 elastic shell elements and 304 nonlinear frame elements with a fiber based section, resulting in 659 FEM nodes with 9,960 branch nodes coupled to 35,960 CFD patch surface nodes within an OpenFOAM model with 4,061,836 points and 3,917,848 cells, the mapping and structural analysis routines comprise about 16 percent of total computational   The FOAMySees API is a collection of files which together form the object OpenSeesPyInstance and the main coupling solution loop for fluid-structure. The coupled solution is initialized by calling the Python routine within Solid1Solver.py, which creates an OpenSeesPyInstance object. This object contains the OpenSeesPy model and metadata about the model, as well as arrays through which data is transferred to and from OpenSeesPy at each timestep. In addition to this, there are sub-routines which write and read checkpoints (OpenSeesPyInstance.stepForward), iterate for a solution (OpenSeesPyInstance.iterate), and rotate branch nodes about their associated FEM nodes for branch mesh updates (OpenSeesPyInstance.RotateTreeBranch).
FOAMySees is structured in such a way that the definition of component OpenFOAM and OpenSeesPy models for coupled FSI analysis does not vary greatly from the definition of individual structural models with OpenSeesPy or individual fluid models within OpenFOAM. Users may define their OpenSeesPy model as usual within the file Solid/buildOpenSeesModelInThisFile.py. This model is coupled with the fluid model contained with the Fluid/subfolder. Model settings such as time-stepping settings and interface coupling options are controlled via files contained within the Solid/subfolder, particularly Solid/timeSettings.py and Solid/ geometrySettings.py. Preliminary analysis such as gravity loading with dynamic relaxation or seismic non-linear time history analysis may be defined as usual within the file FOAMySees/ prelimAnalysis.py. This sub-routine is called before the fluidstructure interaction analysis loop, and the final state of the OpenSeesPy model at the end of the preliminary analysis is utilized as the initial state within the FSI model, including damage represented by material model variable time-histories, initial stresses and deflections.

Coupled analysis file structure
FOAMySees maintains a standard OpenFOAM file structure for a CFD model with dynamic mesh capabilities. A few additional files besides the those required for a standard multiphase OpenFOAM simulation necessary for initializing connections to coupling routines in preCICE are also included. These files include preciceDict, which is stored within the OpenFOAM case subfolder Fluid/system/, and precice-config.xml. The preCICE coupling adapter for OpenFOAM, which reads and writes forces and displacements, is called via the function object preciceAdapterFunctionObject which is defined in the file Fluid/system/controlDict. This function object is available for use within OpenFOAM only after installing preCICE and the associated OpenFOAM adapter (available at https://github.com/precice/openfoamadapter. Note: for incompressible multi-phase OpenFOAM solvers, a modified adapter found at https://github.com/

Frontiers in Built Environment
frontiersin.org moaxm/openfoam-adapter is recommended for proper calculation of interface forces, and was utilized in this study). Coupling settings are managed and defined in the file preciceconfig.xml, which includes the coupling timestep data and partners, residual convergence tolerances, coupling acceleration techniques (Aitken, Broyden, IQN-ILS, IQN-IMVJ), and data exchange mechanisms for coupling (sockets vs. mpi).
Initialization of coupling partners is completed by running a program batch script in bash on the Linux operating system. This script performs all necessary OpenFOAM pre-processing steps, decomposes the case amongst desired computational processors and nodes, and initializes the parallelized OpenFOAM case for coupling with the OpenSeesPy model. It then performs all pre-processing for the OpenSeesPy model, and initializes the OpenSeesPy model and FOAMySees interface for coupling with the existing OpenFOAM model. At this point, a preliminary OpenSeesPy analysis may be performed including any combination of gravity loading, preloading of structure with dead loads, and earthquake loading prior to initialization of the fluid-structure interaction simulation. The coupling between OpenSeesPy and OpenFOAM is controlled via the FOAMySees python class instance initialized at the beginning of the coupled analysis. This object also controls OpenSees configuration, model progression in time, and communication of displacement and force vectors.

Test cases 1.5.1 Hydrostatic Force FSI case
A quasi-static FSI case was developed using a beam with fixed-fixed end conditions. The domain was 10 m in length, 1 m in height, and 0.1 m in thickness. Water was initialized at a depth of 0.5 m, and the modulus of elasticity of the elastic beam was varied to verify expected forces are applied from OpenFOAM to the structural nodes within OpenSeesPy, and that correct displacements are returned from OpenSeesPy to OpenFOAM.

Frontiers in Built Environment frontiersin.org
The modulus of elasticity of the beam was deemed the least consequential variable to change with respect to other variables. Varying modulus of elasticity between cases allows for a consistent domain size, VOF phase-fraction, beam mass, beam element size, cell size, and beam cross section dimensions while providing a linear change in displacement with respect to the elastic modulus. A diagram of the case showing reduction of the fluid loading to a distributed load is shown in Figure 8A. The expected center span displacement for a fixed-fixed Euler-Bernoulli beam under a uniform distributed load is WL 4 /384EI. With an equivalent distributed load of 490.5 N/m (0.5 m × 0.1 m × 9.81 m/s 2 × 1000 kg/m 3 ), a length of 10 m, a moment of inertia of 100 m 4 , and with elastic modulus of 5E+09 Pa, the center span displacement is 2.55469E-08 m. Center span displacements for beams with softer moduluses of elasticity are shown in Table 3, along with the calculated center

Fixed-end frame under shifting hydrodynamic gravity loading
The same domain geometry as described in the Hydrostatic Force FSI Test Case section was utilized to examine the effect of large displacements upon the shifting of water sitting atop fixed-end frames. Frames for this test case possessed a moment of inertia of 8.33E-3 m 4 , a length of 10 m, a cross sectional area of 0.1 m 2 , and elastic modulus ranging from 1E + 08 Pa to 1E + 06 Pa. These properties were chosen to amplify the effects of fluid pooling upon the deflection imparted against frame elements. As the deflection increases, the loading distribution becomes more concentrated near the beam mid-span, as fluid redistributes itself naturally to find a free surface that results in equilibrium between the internal structural forces and external fluid pressure. An illustration of this behavior is shown in Figure 8B. Since the geometry, sections, materials, and loading chosen in this example result in large displacements, multiple phenomena contribute to the behavior of the frame elements, including frame rotational stiffening effects due to tension. OpenSees-only simulations were run with uniform distributed loads due to the fluid pressure in its initial configuration to assess the affects of axial load (i.e., tension stiffening) in this example case. The effect of axial force upon rotational stiffening is evident in Figure 9A-all calculated deflection curves, most notably the curve with "static" loading, tend to decrease in magnitude beyond a deflection to length ratio of 1:100 for this given case. For sections with different ratios of moment of inertia to area, this behavioral inflection point due to frame stiffening effects will also differ.
As the deflection to length ratio increases, so does the effect of the pooling of the water at the center of the frame upon the expected maximum deflection. A prediction loading function for determining quasi-static loading of beams under pooling water at various deflections was developed, as calculated by Eq. 6. These loads were applied dynamically to frame elements in OpenSeesPy identical to those which were used in the FSI simulations, in efforts to provide an independent benchmark against which FOAMySees calculated loads and results could be compared. The prediction loading function estimates the fluid elevation required above a particular deformed frame shape to achieve quasi-static equilibrium between the beam in its deformed shape and the "pooled" water sitting atop the beam. This calculation is completed after every timestep within a dynamic analysis to account for the current deformed shape.
where g gravity δ x ( ) shape function t domain thickness ρ f fluid density SWL still water level L domain length V f fluid volume SWL*L*t ( ) Results from this study are shown via curves termed "OpenSeesPy, Shifting Load", shown in Figures 9A, B. Figure 9B is simply a normalized Figure 9A-each curve was normalized by the "OpenSeesPy, Uniform Load" results calculated in OpenSeesPy to isolate the effects of load redistribution due to fluid pooling upon the frame near the center of its span. A fluid with a higher density and a lower still water level was utilized to both simulate a case within FOAMySees and compare with the quasi-static load shape calculation described previously in this section. This case is labeled as "OpenSeesPy, Shifting Load, H = 0.1 m, ρ f = 3000 kg/m 3 " in Figure 9B. Results are compared with a solution provided by Bonus, 2023, which utilized Material Point Method (MPM) for simulation of both fluid and structure solutions. Figures 9A, B both show strong agreement with the solutions provided by Bonus-slight variations between both solutions exist due to the To demonstrate the dynamic nature of the FOAMySees code, snapshots were selected from a case within this example, particularly the case with the largest displacements and therefore the largest displacement to length ratio. In Figure 10A, four points in time during the simulation are marked as points A, B, C, and D respectively. In sub-figures Figures 10B-E, free surface elevation snapshots are displayed, showing the redistribution of water along the length of the frame and the settlement of the frame into a loading configuration that resulted in much larger deflections than what would occur for a uniform distributed load case.

Laminar multiphase hydrodynamic impact
Winter, (2019) is an analytical FSI case comprised of a 2D multiphase fluid dam break within a tank with initial conditions as shown in Figure 11A which impinges upon a flexible flap affixed to the floor of the tank. The fluid properties correspond to water and air, so that ρ 1 = 1.0 g/cm 3 (1000 kg/m 3 ) and μ 1 = 0.01 g/s/cm while ρ 2 = 0.001 g/cm 3  (1 kg/m 3 ) and μ 2 = 0.0001 g/s/cm.The structure is elastic, 1.2 cm wide by 8 cm tall, with density ρ 0 = 2.5 g/cm 3 (2500 kg/m 3 ), Young's modulus E = 10 7 g/cm/s 2 (1 MPa) and Poisson's ratio = 0.3. Results show strong correlation between the coupled FEM-CFD solution (FOAMySees) and historical results for the same geometry and initial conditions modelled with finite elements entirely or using the particle finite element method. The case was modelled in STAR-CCM + using 8-node bricks with non-linear geometry model enabled to compare results with commercial software. The lateral displacement time history for the top left corner of the flap is shown in Figure 11B. Results from Marti et al. (2006); Idelsohn et al. (2008) Figure 11C. Snapshots from various points along the solution progression show good agreement with the original solution provided by Walhorn.

Post-earthquake laminar multiphase hydrodynamic impact
Identical domain geometry, fluid properties, and fluid initial conditions as described in the Laminar Multiphase Hydrodynamic Impact section were utilized to demonstrate the multi-hazard analysis capabilities of FOAMySees. A fiber section with "non-linearBeamColumn" formulation frame finite elements were used to represent the flap. Material properties of the flap were modified to allow for material yielding and hardening. The flap material chosen was the OpenSeesPy Uniaxial Hardening material, with Young's modulus E = 10 7 g/cm/s 2 (1 MPa), Yield Stress = 10 6 g/cm/s 2 (100 Kpa), and hardening parameters H iso = 0.15 and H kin = 0.25.
The analysis case was run twice; first without a preliminary earthquake analysis to establish a baseline displacement time history with a non-linear material with which to compare, and again with a ground motion excitation preceding the fluid-structure interaction analysis. A plot of the horizontal acceleration over time and flap tip response to the excitation are shown in Figures 12A, B respectively. Ground motion was stopped after 3 s and free vibration of the flap was allowed until 15 s. Numerical damping was utilized in the seismic analysis through the means of changing the Newmark integration coefficients to γ = 0.66 and β = 0.33. The final deformed shape of the structure at the time of termination of the preliminary seismic analysis resulting in a flap tip displacement of −0.0045 m was utilized as the initial configuration of the structure within the fluid-structure analysis. Both cases are plotted alongside the original FOAMySees solution for the Walhorn validation case in Figure 12C. The results show that in this case both material nonlinearity and material history variables play a role in the flap tip displacement response to impact from the breaking dam.

Three dimensional structure under turbulent breaking solitary wave impingement with real structural sections and materials
In order to evaluate the efficacy of the force and displacement mapping routines for complex three dimensional geometries, a 1:5 full scale structure studied experimentally by Sullivan (2021)

Structural model description
The structure was 1.016 m by 1.016 m in plan from column center to column center, comprised of steel frame elements, steel gusset plates, and concrete filled steel tubes. HSS2 × 2 × 1/8 elements were framed horizontally between columns with their center-line at heights of 1.3208 m, 1.8288 m, and 2.3368 m HSS1.5 × 1.5 × 1/4 elements were utilized for chevron bracing from 1.3208 m to 1.8288 m along all four sides of the structure. Panels spanned between HSS2 × 2 elements forming diaphragms within the structural plan footprint. For exact material properties of the experimental specimen and exact structural dimensions, see Sullivan (2021). Uni-axial material properties for structural steel ("Steel02" model, with a yield stress of 344.75 MPa, FIGURE 12 (A) ground motion utilized in preliminary seismic analysis (B) flap tip displacement resulting from base ground motion excitation prior to fluid-structure analysis, (C) comparison between FOAMySees Laminar hydrodynamic impact solutions-Demonstration of material hardening effects due to preliminary seismic loading upon flap tip displacement.

Frontiers in Built Environment
frontiersin.org initial elastic tangent of 200 GPa, strain-hardening ratio of 0.1, and isotropic hardening parameters of a1, a2, and a3 of 18.0, 0.925, and 0.15, respectively) and concrete ("Concrete02" model, with concrete compressive strength at 28 days of −49.64 MPa, concrete strain at maximum strength of −0.00326, concrete crushing strength of −9.93 MPa, concrete strain at crushing strength of −0.01631, ratio between unloading slope and initial slope of 0.1, tensile strength of 4.39 MPa, and tension softening stiffness of 2 GPa) respectively were chosen for modelling materials within OpenSeePy, and gusset plates were not modelled. EqualDOF commands were utilized for connection of structural elements of different formulations within OpenSeesPy. The structure was fixed at its base. Fiber sections with elastic uniaxial materials and 64 fibers each were utilized to represent the composite sections of the columns, which were standard 10.16 cm (4 in) steel pipes with 1.27 cm (1/2 in) thick walls filled with concrete. Elastic beam sections were utilized for the HSS components of the frame. Elastic MITC4 shell elements were utilized for the modelling of the panels spanning across each story, which were given a thickness of 12.7 mm (1/2 in) and material properties of elastic steel. The model was given Rayleigh damping in OpenSeesPy with a value of 7.5% from the frequency of the first structural mode (f 1 ) to five times that frequency (5f 1 ) with Rayleigh mass coefficients of α mass = 0.0 and Rayleigh stiffness coefficients of β tangent 0.0, β initial ζ*5f1−ζ*f1 π*(5f 2 1 −f 2 1 ) , and β committed = 0.0.

Fluid model
The Large Wave Flume is 104.0 m long, 3.7 m wide, and 4.6 m deep, with waves generated by a piston-type wavemaker on the upstream end of the flume. The model was initialized with a still water level (SWL) of 2.0 m, and the structural model and CFD patch surfaces were positioned approximately 40.77 m from the neutral position of the wave maker. A paddle-generated wave with a maximum crest height of 1.45 m and a celerity (the velocity with which a wave advances) of 5.82 m/s within a still water level of 2 m was sent down the flume toward the structure. As the Oregon State University O.H. Hinsdale Wave Research Laboratory Large Wave Flume CFD model with paddle driven waves utilized in this study has been used previously and its information is available in literature, details are omitted here for brevity. See Lewis et al. (2022) for detailed information about the flume geometry and bathymetry, initial conditions of the water within the flume, location of the specimen in the flume, paddle driven wave properties, validation of wave height and velocity with experimental results, CFD boundary conditions utilized, turbulence model properties, and detailed OpenFOAM modelling procedures. The geometry of the structure (elevated structure rather than a concrete shear core) and mesh motion (pointDisplacement) boundary conditions were the only components changed in OpenFOAM between the model utilized in the previously cited study and the model utilized presently. For detailed information about OpenFOAM modelling procedures utilized for paddle generated waves and recommendations for CFD modelling of complex structures, see Winter (2019).

Branch meshes
Two branch meshes of roughly 25,000 nodes each were utilized in this example. The first branch mesh was utilized for displacement transfer from FOAMySees to OpenFOAM, and consisted of a point cloud identical to that of the CFD surface patch nodal locations. The second branch mesh was utilized for force transfer from OpenFOAM to FOAMySees and consisted of a point cloud identical to that of the CFD patch surface face centers. The purpose of this was to ensure direct mapping of forces from their calculated locations (CFD patch surface face centers) to branch nodes with FOAMySees, and direct mapping of FOAMySees calculated branch displacements to the CFD surface patch nodes to demonstrate scalability of FOAMySees without the need for coupling mesh convergence studies for this example. FEM node to

Frontiers in Built Environment
frontiersin.org branch node relationships were generated automatically by means of a clustering algorithm, namely, the KDTree functions within the scipy Python package. The OpenFOAM patch surface node locations were loaded into FOAMySees as a point cloud, and were each clustered to the OpenSeesPy FEM node which was closest in absolute distance. This operation was also completed with the CFD patch surface face centers providing two branch groups for each FEM node, with each controlling either displacement or force transfer to preCICE for application to the OpenFOAM model. A select number of FEM nodes were chosen for coupling to CFD patch surfaces as the full geometry of the structure was not represented in OpenFOAM. The CFD patch surface utilized in this analysis with the OpenSeesPy model used in FOAMySees overlaid is shown in Figure 13A and the displacement branch groupings utilized are shown in Figure 13B.

Computational cost
The model was run for 4 s of simulation time at a timestep of 0.0005 s on 1 UW HYAK Klone HPC Node with 40 processors for a total computational time of 14 h. In comparison to a geometrically

Frontiers in Built Environment
frontiersin.org similar model utilized for benchmarking computational cost with respect to an equivalent CFD model (9,960 branch nodes coupled to 35,960 CFD patch surface nodes) the first second of the simulation with the present model (50,000 branch nodes coupled to 50,000 CFD patch nodes) took roughly 3.5 h, which is nearly twice the computational time of the model described in Section 1.3.7 with nearly identical node and cell counts within the CFD models compared (4,034,512 points and 3,894,298 cells versus 4,061,836 points and 3,917,848 cells). For more information, see Section 1.3.7.

Results
Mesh motion of OpenSeesPy, the FOAMySees displacement branch mesh, and OpenFOAM surface patches along with fluid free surface (isosurface at α = 0.5) overlaid for selected times of the simulation are shown in Figures 14A-C. In Figure 14D a time history of the displacement of the top left corner on the upstream face of the structure in OpenSeesPy and its associated location in OpenFOAM are plotted alongside each other to demonstrate proper transfer of displacement during the simulation. To demonstrate proper transfer of force between coupled models, the forces applied to the branch nodes within FOAMySees were output and summed for comparison with the force calculated within OpenFOAM via functionObjects. The streamwise force time histories for both OpenFOAM and FOAMySees are shown along with experimentally measured forces from Sullivan (2021) for the duration of wave impingement in Figure 14E, and errors for force and impulse transfer across interfaces is shown in Figure 14F.

Conclusion
The presented work offers a highly scalable, versatile, opensource methodology for numerical simulation of non-linear fluid-structure interaction. By allowing for geometric discrepancies between structural and fluid solvers and compensating for this through the use of branch node data communication meshes, high-resolution reduced order centerline element and shell element comprised finite element models may be utilized to simulate the response of structures to fluid loading in three dimensions. This methodology has been validated through several static analysis cases, one high energy dam break case involving a flexible rubber flap placed within the path of flow, and a medium-scale case representing a 1:5 scale structure comprised of steel and concrete subject to breaking solitary wave impingement tested previously at Oregon State University. Additionally, the high energy dam break case was run with seismic excitation preceding the fluidstructure interaction simulation to demonstrate damage state coherence between termination of the preliminary structural analysis subroutines and the initial state of the fluid-structure analysis. More testing of the software and coupling methodology is necessary to determine rate of convergence of solution, memory usage, limits of the branch node data communication mesh method, and timestep stability regimes of the presented coupled solution technique with respect to limiting factors within both the finite element and computational fluid dynamics methodologies. This includes investigating feasibility and advantages of utilizing iterative techniques for interface manifold acceleration convergence such as Aitken Under-Relaxation and Iterative-Quasi-Newton Inverse-Least-Squares (IQN-ILS). Furthermore, more work must be done to validate the solution technique offered by the proposed API with experimental studies conducted at large scales. A detailed investigation of how turbulence modelling affects results was not considered in the present study, thus more research must be completed to assess the accuracy of results provided by FOAMySees when utilizing OpenFOAM fluid solvers which incorporate turbulence modelling.

Data availability statement
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

Publisher's note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.