AN INTEGRATED OPENFOAM MEMBRANE FLUID-STRUCTURE INTERACTION SOLVER

-Abstract. The scope of this paper is to present the design and verification of an integrated OpenFOAM membrane fluid-structure interaction (FSI) solver for small deflections, which employs the finite volume method (FVM) for solving the flow field and the finite area method (FAM) for solution of the membrane deflection. A key feature is that both the fluid and the solid solver operate on a common mesh geometry and are included into a single executable. Although the scope of applicability is narrow due to limitations of the membrane solver at its current state, positive verification results prove the practicability of the design, which allows for lightweight implementation as well as simple data transfers and post-processing.


Introduction
Fluid-structure interaction is a generic term for settings in which the flow-imposed load on a structure leads to its deformation. The flow field itself may in turn be influenced by the structural deformation. If no notable so-called back-coupling on the flow field occurs, the setting is termed one-way coupling, while two-way coupling describes cases where the structural deformation has a non-negligible effect on the flow field. In many cases, the slenderness of surface structures implies a comparatively small stiffness and therefore high deformability, thus making them interesting for FSI studies. Numerical treatment of the said problems is usually demanding regarding computational requirements. Consequently, it has been gaining momentum in the past decades due to the increasing computing power available to the scientific community. This is reflected by the diverse range of publications that cover simulations of FSI problems involving surface structures [1][2][3][4][5][6]. In particular the work on FSI of membranes and plates by Torlak [7], who used a control-volume approach for discretization of the structural part, bears analogies to the work presented in this article.
Common terminology differentiates between FSI solvers of partitioned and monolithic type [8,9]. The former rely on distinct solver applications for fluid and structural part, that exchange boundary conditions at the coupling interface, while the latter treat both field problems combined in a single matrix. Usually, monolithic solvers stand out due to superior stability, but their general implementation causes limitations in applicability to specific field problems. The partitioned solution approach, on the other hand, allows for specialized solvers and meshes for fluid and structural part. In addition, a coupling tool is required for data transfer between the coupling participants, for handling communication and for interpolation between the different meshes. Time advancement in partitioned solvers can be approached in various ways. With small time steps, a single solution of fluid and structure field within each time step may be sufficient without compromising stability, which is termed weak coupling. More economically, strong or implicit coupling can be used with larger time steps while maintaining stability. In this case, within each time step, a nested coupling loop is executed with repeated solution of both coupling participants and exchange of boundary conditions, until compatibility between fluid and solid field is achieved. Due to the artificial added mass effect [10][11][12], special care must be taken when computing FSI with (nearly) incompressible fluids. Added mass refers to the fluid moving with the structure. Since this additional inertial contribution is not known beforehand during the solution process, it is over-or underestimated. The estimation error is termed artificial added mass and may lead to peaks in the pressure field and subsequently numerical instabilities.
The scientific community has used OpenFOAM for partitioned FSI investigations targeting a range of challenging problems. While the OpenFOAM software suite originally focused on flow simulations, also capable large-deflection structural FVM solvers [13] have been implemented. Hence, it is possible to use the framework both for the structural as well as for the fluid part, as is also intended in this article. In 2010, Jasak et al. studied a well-known FSI test case [14], initially proposed by Hron and Turek [15]. Among the newer publications treating both fluid and solid part with OpenFOAM, the work of Tuković et al. stands out, covering a wide range of test cases [16]. Furthermore, code-specific coupling tools or general coupling toolkits like preCICE [17] allow for partitioned FSI with virtually all kinds of solvers, under the condition that a corresponding interface exists.
Apart from the core FVM implementation, the OpenFOAM library features a finite area discretization method. Similar to the FVM, the FAM is a control-volume-based discretization approach, with the control volumes being discrete areas on a surface mesh. It is described in a detailed manner by Tuković and Jasak [18]. In OpenFOAM, both FVM and FAM are based on co-located computing nodes, placed at the volume center or area center, respectively. As a result, the framework can handle finite volume discretizations and finite area discretizations on a single mesh geometry, with the finite area mesh consisting of control volume surfaces from the finite volume mesh. This predestines the OpenFOAM library for the implementation of lightweight solvers for coupled surface and volume field problems. In the past, this capability has already been used for simulation of shallow granular flow [19], surface tension influence [18,20] and thin-film flows [21,22]. In addition, more recent OpenFOAM versions include FAM Kirchhoff plate solvers, which were first introduced by the solids4foam toolbox [23], allowing the computation of simple surface structures.
The scope of this publication is the presentation and implementation of a partitioned approach to FSI using the OpenFOAM framework. It is integrated into a single executable, using an FAM implementation of a membrane equation solver for the solid part and an FVM-based fluid part. Due to its design, it comprises some advantageous attributes of a monolithic solver. Subsequent to the description of the implementation in the following section, verification test cases and corresponding results will be presented.

Coupled Solver Implementation
Using the FVM for computation of a flow field, and the FAM for computation of a deflecting membrane, OpenFOAM can be used to perform fluid-structure interaction simulations using a single mesh geometry for both fluid and solid field. Other than in a monolithic solver, the finite area and finite volume discretizations lead to decoupled systems of equations that are solved independently and sequentially.
The key difference to the common partitioned solution approach is the use of a single mesh geometry and the combination of both solvers in one application, obviating the need for inter-process communication.
On the contrary, all information about the solution progress as well as field data is readily available to both solver subroutines. Fig. 2 depicts the adopted integrated FSI solver setup in contrast to the common partitioned approach with separate solver executables in Fig. 1  2.1. Flow Field Solution. Fluid flow is described by the specific transport quantities mass and momentum. The resulting equations are known as Navier-Stokes equations. On a moving grid, fluxes due to boundary movement have to be considered in the control volume balance. Equations (1) and (2) display the compressible continuity and Navier-Stokes equations for a control volume boundary moving at the velocity u b .
Pressure-velocity coupling is achieved using a combination of PISO [24] and SIMPLE [25] approaches, which is termed PIMPLE in OpenFOAM and iterates through two nested loops within each time step. An outer loop resolves the convective non-linearity, while an inner loop includes successive pressurecorrections. Pressure-velocity coupling according to the PIMPLE scheme is versatile and stable for a wide range of flows, even with comparatively large time steps. By reduction of the number of iteration of either of the loops to one, the basic SIMPLE or PISO behavior can be recovered.

Membrane Equation Solution
. The vital attribute of surface structures like membranes is the low thickness. Their extension in thickness dimension is by orders of magnitude smaller than in the other dimensions. Adopting the assumption of simple distributions (e.g. constant or linear) for stresses and strains in thickness direction permits the integration to resultant quantities, which simplifies the governing equations from three to two-dimensional form. This in turn implies a considerable reduction in the computational effort to solve for their deflection and avoids difficulties of large aspect ratio meshes. The most elaborate theories on surface structures apply for shells with arbitrary curvature in threedimensional space, as depicted in Fig. 3. Within this study, only very simple plain membranes with constant thickness and subject to only normal loading are considered. Membrane theory assumes that the thickness is especially low, such that bending moments -in contrary to e.g. plate or shell theoryare negligible. Instead, only forces due to pretension and curvature counteract the imposed normal load, expressed here as static pressure p, as shown in Fig. 4. The membrane equation for small deflections is given in integral form as equation (3), with the pretension σ 0 , density ρ and thickness d, while m At the boundaries, a fixed support is prescribed by applying a zero deflection Dirichlet condition. As the OpenFOAM FAM implementation provides implicit discretizations for both the Laplacian operator and the second time derivative, the membrane equation can be solved directly for the deflection w using the code: + fam::d2dt2(rho/sigma0, w) 7 ); ¦ ¥ The chosen discretization for the Laplacian was Gauss linear corrected. For the second time derivative in transient cases, only Euler was applicable, which is of first order accuracy. In stationary cases, the time discretization scheme was set to steadyState.

Common Mesh and Data
Transfers. Fig. 5 shows a sketch of a finite volume mesh with an overlaid finite area mesh. The finite areas constitute boundary faces of the finite volume mesh, and the finite area boundary consists of lines between mesh vertices. Centroids of the finite area computing cells are marked as black dots. They coincide with the boundary face centers of the finite volume mesh, which renders data transfers trivial, since no interpolation is necessary. Furthermore, conservation of transferred quantities is intrinsically guaranteed.
On the downside, independent refinement of the finite volume and finite area meshes is not possible. As the finite area mesh is two-dimensional only, it may be assumed that the number of degrees of freedom in the three-dimensional finite volume mesh by far exceeds the number of solid degrees of freedom. Consequently, the computational requirements for the solid part likely are comparatively low, even if the mesh density as dictated by the fluid part is higher than actually required for sufficient accuracy of the solid part. Therefore, the lack of independence in mesh density is deemed to be an only minor disadvantage for coupled solvers adopting the presented design.

2.4.
Mesh Adaption. The OpenFOAM library contains a range of dynamic meshing capabilities, including topological changes [14,26]. For the scope of this study, only small deflections of the membrane, represented by a boundary of the three-dimensional finite volume mesh, are considered. Therefore, mesh smoothing is sufficient in order to absorb the boundary deformation without compromising accuracy or stability due to bad quality cells. The membrane deflection is computed at the area-centered finite area nodes and at the boundary nodes placed on the outer boundary of the finite area domain. For mesh adaption, the first step is the interpolation of the structural deflection to the mesh vertices constituting the finite areas' corners. With the surface mesh vertices deflection known after interpolation, the inner volume mesh has to be adapted in order to smoothly absorb the boundary deflection.
A desirable property of mesh smoothing methods is to approximately conserve the boundary layer resolution, such that the boundary movement is mainly absorbed by inner mesh cells with higher wall distance and larger volume. Even at large deformations, diffusion-based smoothing approaches allow for this by specification of the diffusion coefficient depending on the distance to the deforming boundary [26,32]. The diffusion equation to be solved for the mesh movement x is given as equation (4).
By choice of a large diffusivity Γ at locations where the mesh is to be preserved, and a small value where the deformation is to be absorbed, the deformation can be favorably distributed. In this study, the boundary wall distance method with linear dependence on the inverse distance was adopted, yielding slightly better results in preliminary simulations than exponential or quadratic dependence.

Time
Stepping and Coupling Loop. Convergence of the inner coupling loop is achieved when a compatible solution of fluid and solid part has been found. This is the case, when the participating solvers and the data transfers converge, i.e. the deflection and its derivatives with respect to time as well as the surface load field are consistent between fluid and solid part. Furthermore, a minimum and maximum number of coupling iterations may be specified, such that advancement to the next time step may artificially be delayed or enforced. Although this should ideally not be necessary, it may be useful during case setup when adequate values for the convergence criteria have not yet been found. Weak coupling can be imposed by setting the maximum number of coupling iterations to one and switching off relaxation of transferred quantities. In order to boost convergence of the coupling loop, both adaptive Aitken relaxation [33] and deflection prediction have been implemented, which have proven to be reliable improvements in the literature [2,12,[34][35][36]. While other, usually more complex methods may yield even better results [8], relaxation using Aitken method often provides reasonable improvements and is straight-forward to implement into an existing coupling loop. Relaxation by a constant factor has also been incorporated for comparison purposes. According to literature, the optimal value of the relaxation factor depends on several characteristics of the treated problem [37], hence it cannot generally be expected to be as efficient as more advanced techniques like adaptive Aitken relaxation.
Starting a time step with the previous time step deflection corresponds to zero order deflection prediction (equation (5)) and gives a very basic estimate. With better predictions, substantial speedup of the inner coupling loop can be achieved [34,36]. Based on the assumption of linear or parabolic deflection behavior with respect to time, the formulas for first-order and second-order prediction accuracy as given in equations (6) and (7) for constant time step sizes can be obtained.
Prediction of higher order may be prone to spurious high-frequency disturbances [8]. Hence, the ideal order of prediction remains problem-dependent.
2.6. Convergence Criteria. In total, there are four boolean convergence flags: one for fluid and solid solution each, and a further one for both of the exchanged boundary conditions, deflection and surface load field. Convergence of the fluid and solid solvers is judged by the remainder after solution of the discretized governing equations. In addition, compatibility between both participants must be achieved, implying that data transfers must have converged. This is assessed by the change of the transferred quantity fields between subsequent coupling iterations, which must be small compared to the change between time steps. In order to simplify the quantification by condensing the field differences into a single number, mathematical norms of the difference fields are computed, such that the resulting condition for convergence may be written as equation (8). Indices t mark time steps, while indices n stand for the coupling iterations within a time step. The type of norm applied is denominated by the subscripts a and b.
Two norms adopted in the implementation are the maximum absolute norm (peak-to-peak value) and the root mean square value (RMS), corresponding to a weighted euclidean norm. The latter bears the advantage of giving a good representation of the average order of magnitude with emphasis on the larger values. On the other hand, the peak-to-peak value corresponds to the upper bound of the absolute differences, implying that although it may not be good at representing the mean order of magnitude, it effectively prohibits smoothing out of unphysical peaks. Convergence indication must work reliably in the sense that neither stalling nor premature convergence indication should occur. While it is hard to find an indicator that works absolutely reliably, for this study good results were obtained when adopting the RMS value for quantification of the change between time steps and the maximum absolute value for quantification of the change between coupling iterations. This ultimately leads to equation (9) for the relative error.
The deflection field and its changes can be expected to be rather smooth, due to damping of the motion by structural inertia. On the other hand, the load field, due to fine turbulent motions, is potentially more prone to show individual local differences and therefore larger peak-to-peak differences between iterations values, depending on the fluid velocity discretization. For the scope of this study, the convergence criterion for the deflection was set to a relative error according to equation (9) of 0.0002, while for the surface load field, a criterion as high as 0.2 had to be chosen in order to prevent significant delay in convergence without noticeable gain in accuracy. In both cases, for the computation of the relative errors, the computed transfer quantities without under-relaxation were used as reference, in order to avoid erroneous premature convergence indication caused by strong under-relaxation.

Test Cases Setup
Verifications were individually conducted for the membrane solver and for the coupled FSI solver. Welldefined test cases on simple geometries were set up in order to allow the comparison and to reduce further erroneous influences. For testing the membrane solution only, the coupled solver was configured to apply a predefined load instead of fluid pressure from the neighboring finite volume cells on the membrane.
3.1. Stationary Membrane Test Case. For a range of stationary cases, solutions to the membrane equation can be derived by analytically [38,39]. In order to check the stationary membrane solution, the maximum deflection of a circular and quadrilateral membrane under uniform loading, occurring at the center, is compared to analytical results. The formulas for membranes with diameter D or side length L, respectively, are given in equations 10 and 11.
While the quadrilateral membrane employs a regular Cartesian grid, the circular membrane was discretized either with a semi-regular quad-dominant or with a completely irregular triangular mesh, as depicted in Fig.7. The control area grid density indicates the fraction of the side length or diameter, respectively, which was set as characteristic area (side) length for grid generation. By analysis of different meshes, the required grid density for a given accuracy can be estimated.

Transient Membrane Test Case.
In order to verify the transient behavior, a reference solution was generated using ANSYS Mechanical and compared to the results obtained with OpenFOAM. A time and space-independent load was applied onto a quadrilateral membrane of 1 m side length without initial deflection. Due to second-order accurate time integration, the reference solution was obtained using a time step size of 4 ms, while for OpenFOAM the time step size was set to only 0.25 ms. The membrane was modeled with a mesh of 40 times 40 SHELL281 second-order finite elements using ANSYS and 80 times 80 finite areas in OpenFOAM, respectively.

Membrane Solver Oscillation Test
Case. The exact eigenfrequencies and normal modes of quadrilateral membranes of side length L can be derived analytically [40] as given in equation (12).
Hence, from a comparison of analytical to numerical solution with varying time step sizes, the influence of numerical damping on the time dependent solution may be assessed. The oscillation test case starts with a deflection according to either the first or the third normal mode. From that on, without a load boundary condition applied onto the membrane, the time-dependent deflection of the membrane is simulated. The analytically derived first and third eigenfrequency for the quadrilateral membrane under consideration are f eigen,1 = 5 √ 2 Hz f eigen,3 = 15 √ 2 Hz which corresponds to a cycle time of roughly 0.141 s and 0.047 s, respectively. Simulation runs resolving each cycle with 100, 200, and 400 time steps give insight into the required time stepping size for a targeted maximum numerical damping of the oscillation amplitude.

Coupled Solver Test Case.
Other than for the stationary test cases, no suitable coupled problem with a reference solution exists. Therefore, the FSI capabilities of the ANSYS software suite were relied on for creating a reference solution for the academic test case described in the following. Fig. 8 displays the cuboid domain geometry, which is aligned with the coordinate system. Very similar fluid domain meshes with 0.5×10 6 and 0.56×10 6 control volumes were used for OpenFOAM and ANSYS, respectively. The time step size was set to 0.05 ms in all simulations. A nozzle of quadratic cross-section of 0.04 m side length points downwards on a quadrilateral membrane of 0.2 m side length with a material density of 1000 kg m −3 . The outflow boundary is at one side of the cuboid. Fluid with a density of 1 kg m −3 and a kinematic viscosity of 10 −4 m 2 s −1 enters with a velocity of 20 m s −1 , leading to a nozzle outlet Reynolds number of 8000. Although turbulence has to be expected under these conditions, turbulence modeling was disabled, because differing details of turbulence model implementations might lead to spurious deviations between the results of the OpenFOAM and ANSYS simulation. Deviations from exact physical behavior were accepted in favor of good comparability between the simulation frameworks and reduced computational requirements.
The geometry of the test case with central impingement at the membrane leads to shear forces likely nearly canceling over the membrane in an integral sense. Hence, neglecting the shear forces in the modeling of the solid part is supposed to be a minor source of error in the investigated setting.
Except for the structural part of the ANSYS solution with second-order accuracy in time, all computations were conducted with first-order accuracy in time, since no higher-order discretization scheme was readily available for the FAM discretization of the second time derivative in the OpenFOAM membrane equation solver. For spatial discretization, on the other hand, schemes of second-order accuracy were chosen, in particular linear upwind for convection. The PIMPLE fluid solver configuration was set up to use two pressure correction iterations, and as many as necessary outer iterations in order to fulfill the flow convergence criteria. In line with suggestions from literature [2,3], different configurations for deflection transfer relaxation were tested, as specified in Table 2, while the surface load field relaxation was always disabled by setting the relaxation factor to 1.

Results and Discussion
All results were obtained on desktop machines using a single processor core. As no dedicated computing cluster was used, no wall clock time comparisons of the simulations are given, because these would not necessarily reflect the computational cost.

Membrane Solver Stationary Test
Cases. Fig. 9 shows the relative deviation between analytical and numerical maximum deflection for a range of control area grid densities on a quadrilateral and circular membrane. The results comply with the expectation that finer meshes lead to better results. For the quad-dominant circular membrane, the error changes sign at a grid density of 20, causing an apparent increase in accuracy for this data point. In general, regular grids as the one employed for discretization of the quadrilateral membrane seem to be advantageous when compared to irregular grids. However, a switch in discretization scheme of the Laplacian in equation (3) from Gauss linear corrected to Gauss linear fourth, with potentially increased accuracy of the line-normal gradient, did not have an effect Figure 8. Coupled test case brick geometry with inlet and outlet marked by arrows. The membrane surface structure lies centered beneath the nozzle on the results. The reasons are not clear to the authors, but it may indicate that the error stems mainly from the handling of boundary conditions and not from inside the finite area domain. In all three cases, a deviation of less than 2 % is reached with a control area grid density of as little as 20, making the implementation seemingly good enough for a range of scenarios with moderate accuracy requirements.

4.2.
Membrane Solver Transient Test Case. The results of the transient test case under constant loading obtained with OpenFOAM in comparison to the results obtained with ANSYS is depicted in Fig. 10 and 11. The former, displaying the central point deflection over time, shows that initially the deviation between both curves is quite small, but accumulates over time. Apparently, a slight deviation in oscillation frequency leads to a growing phase shift with increasing number of oscillation cycles. Apart from that, also the amplitude appears to be damped in the OpenFOAM solution, although in general the absolute value of the deflection is captured well. In Fig. 11, the snapshot of the cross-sectional deflection after 140 ms more clearly displays the differences across the width of the membrane. Due to different scaling of the y-axis, the results appear more pronounced in Fig. 11 when compared to Fig. 10.
As the results from the membrane oscillation test case in the following section suggest, the differences are likely due to the use of the highly dissipative first order Euler time discretisation which is used in OpenFOAM, in contrast to the second order scheme in ANSYS. Considering that the deviations of the OpenFOAM membrane solution to the ANSYS reference solution are initially small and accumulate slowly over time, the authors conclude that the results are acceptable. Table 1 summarizes the findings of the oscillation test cases conducted for the first and third normal mode. For two of the investigated cases, Fig. 12  quadratic membrane circular membrane, quad mesh circular membrane, tri mesh   In all cases, as expected for first-order accuracy in time, the damping due to numerical error drops proportionally by reduction of the time step size. For higher normal modes, damping at a given time step size is stronger, but apparently solely due to the relatively coarser resolution in time. The relative spatial resolution of the half-wave -thirty finite areas for the first and ten finite areas for the third normal mode -apparently has no effect. Hence, the significant damping can be attributed to the time discretization of first order. In order to allow for higher accuracy at larger time step sizes, a discretization with second-order accuracy for the second time derivative in the OpenFOAM FAM implementation is desirable. Otherwise, for resolution of fine-scaled oscillations, the time step size dictated by alleviation of numerical damping avoidably increases the computational cost. Coupled Solver Test Case. The coupled OpenFOAM FSI solver was run using several configurations, in order to assess the impact of different coupling configuration parameters on coupling cycles and on computation time, as shown in Table 2. Differences in the results among the OpenFOAM simulations were negligible, except for the one-way coupling case without mesh deformation. This indicates that the adopted convergence criteria for fluid and solid part as well as for convergence of data transfers were chosen sufficiently strict. The Courant number of the flow solution showed a peak value of around 1.2 with a domain average of 0.09 for fully developed flow. Maximum y + values at the membrane were close to 2.4, while the average was as low as 1.5. Since no turbulence modeling was used, this is deemed to not have any impact on the simulation result. In all two-way coupled OpenFOAM cases, approximately 84 % of the computing time were spent for the flow field solution and the remaining 16 % for mesh movement related computations, including flux correction. The computing time spent for the membrane deflection and data transfer between fluid and solid part was negligible. This confirms, that at least for this test case the mesh dependency is to be considered only a minor drawback of the coupling approach.

Membrane Solver Oscillation Test Case.
In line with expectations, the results in Table 2 indicate that higher order prediction effectively decreases the total computational cost due to effective reduction of the number of required coupling iterations per time step. Even though the time step size was quite low, second-order prediction proved notably better than first-order prediction. Deflection relaxation using the Aitken method, on the other hand, only provided an improvement of roughly 15 % when compared to a constant relaxation factor close to one. Albeit in literature conditions were described in which under-relaxation by a constant factor slightly outperforms the Aitken method [3], most of the times it showed notable advantages when compared to cases where a lower constant under-relaxation had to be applied for stability [2,3,35]. This suggests that for general applications the Aitken method is preferable, since the additional computational cost is negligible. Fig. 13 compares the one-and two-way coupling results for the center point deflection over time along with the reference solution computed using ANSYS. Interestingly, even though the relative deviation of the membrane peaks at a maximum of only roughly 1 %, already notable differences between two-and one-way coupling (with and without dynamic mesh) become apparent.  Figure 13. Simulation results for central point deflection over time with one-and twoway coupling. In spite of the relatively small absolute membrane deflection, two-way coupling appears to be mandatory for accurate results. fails to replicate the dynamic movement of the membrane central point, not only with respect to the amplitude, but also to fidelity of the movement. Even though the two-way coupling case does not exactly replicate the reference solution, there is in general good agreement between the ANSYS and two-way OpenFOAM results. The differences may at least partially be explained by the weaknesses of the membrane solver revealed in the previous sections, especially by the use of the dissipative first-order time scheme in OpenFOAM, in contrast to the second-order scheme applied in the structural solver of the reference solution. Furthermore, deviations might be attributed to differences in implementation details in the fluid parts.
In summary, the initially small and slowly increasing deviations indicate that the coupled solver in the current state is applicable in the considered scenario for the prediction of the membrane motion. Due to the rather special nature of the test case, it is hard to tell if there are further limitations, which remained undiscovered in this study. For example, these might arise from shear forces acting on the membrane if there is notable cross-flow. Nonetheless, it could be shown that the concept of an integrated FSI solver according to the presented design is viable, although the current state suffers from several weaknesses.

Conclusion and Outlook
An integrated partitioned membrane FSI solver using OpenFOAM has been presented. It is based on coupling of finite area and finite volume discretizations on a common mesh geometry for solid and fluid part, respectively. The main advantage of the adopted integrated design into a single executable is the seamless interaction between solid and fluid solution routines with fast and conservative data transfers. Furthermore, the whole result data is available in OpenFOAM data format, which allows for combined post-processing of fluid and solid part. On the downside, the coupled solver's main limitation is the narrow scope of applicability of the solid part to simplistic membranes with small deflections. In addition, the common mesh geometry prohibits independent mesh refinement of fluid and solid part.
Verification results for the FAM membrane solver indicate that a mesh length scale of around a twentieth of the membrane length scale suffices in order to adequately resolve coarse deflection patterns. A major limitation at the presented state is the lack of an FAM discretization for the second time derivative of higher than first-order accuracy. This may lead to notable damping of oscillations in the solid part even if the time step size is chosen two orders of magnitude smaller than the oscillation which is to be resolved.
For verification of the coupled solver, a reference solution using a commercial FSI framework was generated. While the test case is to be considered academic only, the agreement between the two-way coupled solver and reference solution are in good accordance, although the fidelity of the reference solution is not correctly reproduced. A solution obtained using one-way coupling (neglecting mesh deformation) showed significantly larger deviations to the reference solution than the solution obtained using two-way coupling, although the membrane deflection was two orders of magnitude lower than its lateral dimension. While for the considered test case, deflection prediction significantly reduced the required computing power, Aitken relaxation for deflection transfer did not yield notable improvements when compared to constant relaxation with a factor close to one. In summary, the verification results suggest that the implementation bears the potential to efficiently solve simplistic FSI problems, if a coarse resolution of the flow-induced motions is acceptable. More notably, the feasibility of the coupled solver design could be proven.
A possible topic for which the FSI solver in the current state may already be suitable is the simulation of synthetic jet generators with small membrane deflections. The FAM implementation of more elaborate surface structures, like plates or membranes with large deflection, is necessary in order to make the approach applicable for a broader range of test cases. A starting point would be the incorporation of the plate equation solvers available in recent OpenFOAM versions. In order to integrate the solver more seamlessly into the OpenFOAM hierarchy, the implementation of the solid part bound to a finite area region appears reasonable.
Furthermore, the characteristic feature of the coupled solver, that fluid and solid part are contained in a single executable bears the advantage of each solver having full insight into the other coupling participant's convergence progress. This may be exploited to adaptively tighten the convergence criteria according to the convergence progress during the inner coupling iterations of each time step, reducing the computational requirements.