On aerodynamic force computation in fluid–structure interaction problems — Comparison of different approaches
Introduction
The fluid–structure interaction (FSI) is a very interesting multi-physical problem with many technical applications as aeroelasticity, hydroelasticity, operating life of compressor blades, see [1], [2], as well as many biomechanical applications, e.g. hemodynamics or human phonation, [3], [4]. The understanding of the FSI phenomena and their further subsequent optimization depend on an accurate solution of these FSI problems. Since the FSI problems have very rarely an analytical solution, reliable and fast numerical simulations are of great interest. However, numerical simulation of FSI is a computational challenge since the governing equations for solid and fluid regions are different and the algorithm needs to solve the locations of solid–fluid interfaces simultaneously with the dynamics in both regions where the kinematic (e.g., non-slipping) and dynamic (e.g., stress matched along the normal to solid–fluid interface) boundary conditions are imposed at the interface.
The FSI problem is the non-linear problem consisting of two subproblems, the fluid flow problem and the elastic structure deformation problem, coupled by the boundary conditions at the common interface. The coupling between the two subproblems often leads to a highly nonlinear problem, which needs to be addressed by a robust and efficient solution technique, see [5]. The numerical treatment of FSI problems is complicated due to the need for mixed description — while the spatial (Eulerian) description is suitable for the fluid, the natural description for the solid part is the material (Lagrangian) description. Thus usually a mixed arbitrary Lagrangian–Eulerian (ALE) description is employed introducing an additional nonlinearity into the resulting equations. In general, the numerical schemes for solving the FSI problems may be classified into two approaches: the monolithic method, see e.g. [6], [7], [8], [9], [10] or the segregated method, see e.g. [11], [12], [13], [14]. The monolithic solvers are usually considered to be more robust, whereas computationally more demanding. Their advantage is that interface conditions can be automatically satisfied as these are included in the fully coupled (weak) formulation, see e.g. [9]. In a monolithic approach, the system of nonlinear algebraic equations arising from the discretizations is solved typically using a variant of Newton’s method, see e.g. [6], [8]. In [7] a block preconditioner is presented for the efficient solution of the linear systems that arise when employing Newton’s method.
Even though the monolithic – or fully coupled – approach is usually considered to be more robust, still the segregated approach, or partitioned scheme, is widely used for multi-physics problems due to its modularity: the specialized well-tuned solvers for the subproblems can be coupled only by the exchange of the data at the interface. In practice, usually the Dirichlet–Neumann (DN) partitioned scheme is used, see e.g. [11], where a simple and robust fixed-point FSI solver with dynamic relaxation was presented suitable for a wide range of applications. On the other hand the DN partitioned schemes are sensitive to the artificially added-mass effect which can cause the coupling iteration to diverge, see e.g. [15], [16]. In the last decade, the interface quasi-Newton method was developed in order to control the added-mass effect and simultaneously to speed up convergence, see e.g. [17]. Although such approach has several advantages, still in many applications the standard DN coupling is used, see e.g.[18], where the numerical investigation of free dynamic response of a rigid, symmetric airfoil was performed using standard commercial software. It was confirmed in [18], in particular, that the choice of the numerical parameters – such as domain size, grid size, and time step – can have a big influence on the accuracy of the numerical results.
In all the above mentioned approaches, the numerical procedure used for evaluation of the aerodynamic quantities seems to be somehow overlooked. In many papers the procedure of the evaluation of the aerodynamic forces (AF) is not specified even though this is one of the key components of the algorithm and there exist several possibilities of how the AF can be numerically approximated, [19]. The AF is mathematically defined as the surface integral of the AF density represented by the fluid stress tensor. From the mathematical point of view the evaluation of the surface integral is not well established for the standard choice of the functional spaces and for pressure and velocity, respectively. From the numerical evaluation point of view, if we omit the possible discontinuities of the pressure, the viscous part of the tensor determined by the symmetric part of the velocity gradient needs to be approximated e.g. by numerical differentiation formula and it is called the derived quantity, i.e. the quantity which can be obtained by postprocessing of the primary variables, see [19], [20]. However the numerical differentiation, or more precisely the straightforward extrapolation of the stress tensor components from the interior domain on its boundary, can significantly pollute the solution by numerical errors and it possibly reduces the theoretical accuracy of the method.
In literature there are known two approaches applicable for the AF evaluation within the finite element method (FEM) framework of how to overcome the above mentioned possible loss of accuracy. The first approach is motivated by the superconvergence phenomenon observed for FE approximations of a scalar linear problem, see e.g. [21], [22], [23]. Here, it was shown that derivatives of the approximate solution at appropriately chosen points of the elements (e.g. in triangle centroids) have a higher accuracy compared to the standard error estimates. This phenomenon is called superconvergence. One way how to exploit is to compute quantities of interest at these points [23] or to use the averaging of derived quantities directly computed from numerical solution over local patch of surrounding elements, see [24]. A technique followed here is to use the least-square method instead of pure averaging, see e.g. [25]. These procedures are common components of the FEM commercial softwares like Comsol or Abaqus, etc. [26], [27], and they improve theoretical accuracy for sufficiently fine meshes, see [23]. The second alternative approach of the AF evaluation is provided by their weak reformulation, which is sometimes called the consistent force formulation, see [28]. In this case the evaluation of the forces is rewritten with the aid of the Green theorem in combination with a suitable test function in the form using only volume integrals consistent with the weak form of the problem. Therefore, this method retains the accuracy of FE velocity approximation, see error analysis in [20]. The application of this method for a simplified FSI problem was published e.g. in [15], [29]. This approach seems to be also available in the Comsol toolbox, see [26], where however the detailed mathematical description is missing.
There are several FSI problems motivated by biomechanical applications e.g. by the flow-induced vibrations of vocal folds, see e.g. [30], [31]. Such a problem can be modeled as a FSI problem formed by the elastic subproblem described using the linear elasticity theory and by the fluid flow modeled by the incompressible Navier–Stokes equations in the ALE formulation allowing to incorporate the effects of the time dependent fluid domain, see [32]. The FEM is used for the numerical approximation of both FSI subproblems and the strongly coupled partitioned approach is implemented in order to increase accuracy and robustness of the numerical method especially for the simulation of the flutter onset, see [15], [31], [33]. Let us mention that in this case the classical Theodorsen method is not applicable due to complicated geometry of vocal folds, but the linear stability theory, which in the discrete version is based on the eigenvalues determination of associated system matrix, was applied e.g. in [30], [34] using several simplifying assumptions (as simplified vocal folds geometry, the potential flow model and a simplified AF). This simplified approach was compared with the FSI simulation based on the continuum flow model discretized by the FEM in [35], [36] with a good yet partial correspondence. The determination of aeroelastic stability boundary represented by the flutter velocity is an important task also in engineering areas like e.g. airfoils and turbomachines, [1], [2].
In this paper we pay the main attention to the problem of numerical approximation of the AF namely in the FSI problems, where the effect of the accuracy of approximation of AF can be expected to be substantially stronger due to the substantial dependence of the interface position on the acting AF and vice versa. Here, the three possibilities of AF calculations are described motivated by the above mentioned techniques, and these numerical procedures are tested on several benchmark problems similar to those in [37]. First, the classical DFG benchmark problem is considered, see [38]. However, for the FSI problems the evaluation of the AF can be strongly influenced by the substantial dependence of the interface position on the acting AF and vice versa, which brings motivation for the second test case of harmonically oscillating structure. Further, the AF evaluation techniques are applied for solution of a relevant FSI problem motivated by [36] and the influence of the applied techniques on the flutter velocity is tested. The main novelty of the present paper lies in significantly extended comparison of AF calculations studied for the benchmark fluid problems in [37] to the FSI problems motivated by the flow-induced vibrations of vocal folds, see e.g. [30], [31], [34].
The paper is structured as follows: The mathematical description of the FSI problem is given in the first section of this article. The second section contains brief description of the numerical scheme with a special attention paid to the three considered techniques of aerodynamic force calculation (within the FEM). The numerical results presented in the third section compare different AF calculations on four test cases. The first two cases consider static and dynamic benchmark problems of a cylinder in a cross-flow, while the next two FSI simulations investigate the structural energy dissipation and accumulation including determination of the critical inlet velocity for the considered problem setting.
Section snippets
Mathematical model
Let us consider a two-dimensional FSI problem consisting of an elastic body represented by domain and flow domain represented by domain and their interface , see Fig. 1, where the time dependences are highlighted by lower index . At , both the domains are given as part of initial conditions of the FSI problem and they are called reference domains, i.e. , . Since the Lagrange coordinates are used for description of the elastic body, the time index is further
Numerical model
For the spatial discretization of the subproblems (1), (5) the FEM is utilized. For the purpose of time discretization the time interval is divided into equidistant time steps with and the Newmark and backward differentiation formula (BDF) is used.
Numerical results
The described aerodynamic forces computation procedures are tested on four tests and the results are compared. First, these techniques are applied for computation of the drag and the lift coefficients for the flow around the cylinder and their convergence on refined meshes is investigated. Furthermore, a test with the oscillating cylinder is realized testing the methodology on moving meshes. This test is followed by two FSI cases motivated by simulation of human vocal folds. In the first FSI
Conclusion
This paper presented numerical techniques applicable for evaluation of aerodynamic forces in the context of a FSI problem. The mathematical model of a fluid–structure interaction (FSI) problem was presented and its numerical solution based on the finite element method was described, where the so-called minielement, i.e. P1-bubble/P1 element, was used for approximation of the velocity and pressure. Special attention was paid to different numerical realizations of aerodynamic force (AF)
CRediT authorship contribution statement
Jan Valášek: Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing – original draft, Writing – review & editing. Petr Sváček: Conceptualization, Funding acquisition, Methodology, Supervision, Writing – review & editing.
Acknowledgments
Authors acknowledge support from the ESIF, Czech Republic, EU Operational Programme Research, Development and Education, Czech Republic, and from the Center of Advanced Aerospace Technology, Czech Republic (CZ.02.1.01/0.0/0.0/16_019/0000826), Faculty of Mechanical Engineering, Czech Technical University in Prague as well as support provided by Premium Academiae of Šárka Nečasová provided by the Czech Academy of Sciences .
References (48)
- et al.
Special issue: Advances in computational methods for fluid–structure interaction and coupled problems - preface
Comput. Methods Appl. Mech. Engrg.
(2001) An efficient solver for the fully coupled solution of large-displacement fluid–structure interaction problems
Comput. Methods Appl. Mech. Engrg.
(2004)- et al.
An efficient preconditioner for monolithically-coupled large-displacement fluid–structure interaction problems with pseudo-solid mesh updates
J. Comput. Phys.
(2012) - et al.
Performance of a new partitioned procedure versus a monolithic procedure in fluid–structure interaction
Comput. Struct.
(2009) - et al.
Validation of recipes for the recovery of stresses and derivatives by a computer-based approach
Math. Comput. Modelling
(1994) On energy conservation for finite element approximation of flow-induced airfoil vibrations
Math. Comput. Simulation
(2010)- et al.
Lagrangian-eulerian finite element formulation for incompressible viscous flows
Comput. Methods Appl. Mech. Engrg.
(1981) - et al.
Stabilized finite element schemes with LBB-stable elements for incompressible flows
J. Comput. Appl. Math.
(2005) - et al.
A Modern Course in Aeroelasticity
(2004) - et al.
On the use of a flux-splitting scheme in the numerical flutter analysis of a low-pressure turbine stage
Acta Polytech.
(2021)