A moving-grid approach for fluid-structure interaction problems with hybrid lattice Boltzmann method
Giovanni Di Ilio, Daniele Chiappini, Stefano Ubertini, Gino Bella, Sauro Succi
AA moving-grid approach for fluid-structure interaction problems withhybrid lattice Boltzmann method
G. Di Ilio ∗1 , D. Chiappini , S. Ubertini , G. Bella , and S. Succi University of Rome ’Niccolò Cusano’, Via don Carlo Gnocchi 3, 00166, Rome, Italy University of Tuscia, Largo dell’Università snc, 01100, Viterbo, Italy University of Rome ’Tor Vergata’, Via del Politecnico 1, 00133, Rome, Italy Istituto Applicazioni Calcolo, CNR, Via dei Taurini 19, 00185, Rome, ItalyJune 25, 2020
Abstract
In this paper, we propose a hybrid lattice Boltzmann method (HLBM) for solving fluid-structure inter-action problems. The proposed numerical approach is applied to model the flow induced by a vibratingthin lamina submerged in a viscous quiescent fluid. The hydrodynamic force exerted by the fluid on thesolid body is described by means of a complex hydrodynamic function, whose real and imaginary parts aredetermined via parametric analysis. Numerical results are validated by comparison with those from othernumerical as well as experimental works available in literature. The proposed hybrid approach enhances thecapability of lattice Boltzmann methods to solve fluid dynamic problems involving moving geometries.
During the last decades an increasing attention has been addressed to model fluid-structure interaction (FSI)problems where the mutual actions played in between a viscous fluid and an immersed object determine theevolution of both fluid and solid motion. One of the most studied test-case is the oscillation of a lamina intoa quiescent fluid, which may have practical applications in several technical branches, such as atomic forcemicroscopy [1, 2], sensors and actuators based on micro-mechanical oscillators [3, 4, 5, 6], cooling devices forelectronics [7], underwater propulsion [8, 9, 10], micro-electro-mechanical systems [11, 12] and energy harvestingthrough smart materials [13, 14].Due to the large interest on this class of problems, several experiments as well as numerical studies have beencarried out in order to characterize the behavior of such a system and predict the forces exerted on the structureby the fluid. Interactions between the fluid and the structure are commonly treated in terms of added massand hydrodynamic damping. These components may act both in-phase and out-of-phase with respect to theacceleration of the oscillating body.The works available in literature span a wide range of Reynolds numbers and oscillation amplitudes. In par-ticular, such analyses concern very low Reynolds number (Stokes flow regime), where the inertial effects areneglected [15], as well as high frequency/high oscillations regimes, where non linearities and second order effectsplay a significant role in force exchange process [16].The effects due to the oscillatory motion of objects submerged in a fluid are typically described by means ofthe Keulegan-Carpenter (KC) number, which is a non-dimensional quantity defined as the ratio between theoscillation amplitude and the characteristic length of the body (usually the cross-section), multiplied by 2 π .Therefore, the low Reynolds number regime corresponds to the limit of KC →
0. This case has been firstlystudied by Stokes [17] and Tuck [15], who analyzed it numerically through linearized Navier-Stokes equationsfor small amplitude oscillations. By adopting a boundary integral formulation for the stream function, Tuckdetermined the numerical solution for the hydrodynamic load on rigid bodies of different geometries submergedin an unbounded fluid. Approaches similar to [15] were also employed in [18] to study the effects of the presenceof a solid surface in the proximity of the oscillating lamina, and in [19], where the case of a vibrating laminaunder a free-surface is studied. Under the same assumption of infinitely small amplitude oscillations, Sader [20]presented a detailed theoretical analysis of the frequency response of a cantilever beam immersed in a viscousunbounded fluid and excited by an arbitrary driving force. In his work, he provided a mathematical practical ∗ Electronic address: [email protected] ; Corresponding author a r X i v : . [ phy s i c s . c o m p - ph ] J un ormulation of that FSI problem, by introducing a complex hydrodynamic function, whose real and imaginaryparts correspond to the added mass and the viscous damping, respectively, acting on the oscillating body.The findings obtained in all of the aforementioned studies are valid for values of KC close to zero. However,at higher KC, vortex shedding and convective effects become predominant and start to dominate the fluid flowphenomena. Therefore, for these cases, the simplified approaches based on linearized Navier-Stokes equationsis no more valid. In order to cope with this issues, many experimental tests have been carried out and severalnumerical models have been proposed. Aureli and Porfiri [21] and Bidkar et al. [22] studied the nonlinear vibra-tions of cantilever beams of rectangular cross-section undergoing large amplitude oscillations. In a recent work,Aureli et al. [23] extended the validity of the classical hydrodynamic function developed for small amplitudeoscillation by providing a correction which takes into account nonlinear effects. Their results were obtainedvia computational analyses and validated by comparison with an ad hoc designed experimental study. Morerecently, Tafuni and Sahin have presented results obtained through Smoothed Particle Hydrodynamics for thecase of a lamina undergoing large amplitude oscillations in unbounded domain [24] and under a free-surface[25]. In particular, in [24] the authors propose a novel formulation for the hydrodynamic function. Some worksconcerning vibration of laminae in viscous fluids were also performed with lattice Boltzmann methods (LBMs)[26, 27, 28, 29, 30, 31, 32].Focusing on LBMs, previous works have been developed within the framework of standard schemes, based onthe Bhatnagar-Gross-Krook approximation (LBGK) on regular cartesian grids. In particular, Falcucci et al. [33]were the first studying this kind of problem with a LBM and they proposed an innovative refill procedure for thelattice sites in the proximity of the oscillating lamina. Although the presented method was proven to representsatisfactorily the fluid flow physics, the numerical prediction of the hydrodynamic forces slightly overestimatesliterature results. Further, a combined lattice Boltzmann and finite element method for FSI problems has beenproposed by De Rosis et al. [34]. In [35], the same authors presented a coupled lattice Boltzmann-finite elementapproach with immersed boundary method (IB-LBM). They later applied such a method to the case of a thinlamina in a quiescent viscous fluid [36]. Their results show a significant improvement, in terms of accuracy,in the prediction of the hydrodynamic load acting on the lamina, with respect to the previous lattice Boltz-mann simulations. In Shi and Sader [37], an ad-hoc LBM implementation was developed to analyze harmonicoscillatory Stokes flows in the frequency domain. The proposed method was validated by simulating the one-dimensional oscillatory Couette flow between two plates, the two-dimensional flow generated by an oscillatingcircular cylinder and the three-dimensional flow induced by an oscillating sphere. Moreover, Colosqui et al. [38]investigated on the high-frequency oscillations of electromechanical resonators operating in gaseous media, andthey obtained a good agreement with experimental observations.In order to enhance the capability of the lattice Boltzmann method to simulate FSI problems, we propose anextension of the hybrid lattice Boltzmann method [39, 40] to moving grids (MG-HLBM). In this method, theregion around the body is treated by means of a finite-volume LB scheme, while the outer domain is solved viaa traditional LBGK approach on structured cartesian grid.The benefits deriving from HLBM have been already shown in a recent work [39]. Here, we explore the capa-bilities of such a method to deal with moving unstructured grids into a fixed Cartesian domain. The resultsare compared with those available in literature, showing a good agreement and a significant improvement withrespect to previous works using original LBM implementation. Such an improvement is mainly due to tworeasons. First, the proposed moving-grid approach allows an efficient mesh refinement. The local refinement,which is realized by means of an unstructured grid and is applied in the proximity of the solid body, followsthe motion of the object. This leads to a computationally efficient and accurate method, since the outer fixeddomain can be represented by a course mesh, while the solid geometry is described by a fine mesh. This featureof the MG-HLBM represents an advantage with respect to the numerical methods based on a fixed mesh, whichare, in general, computationally expensive, since they require a dynamical re-meshing to capture accurately theflow physics around the moving object. Second, the presence of an unstructured halo in the surrounding of theoscillating lamina allows an accurate evaluation of the forces, since computational nodes are located on the realbody surface. This prevents from the implementation of error-prone interpolation strategies.The work is organized as follow. In Sec. 2 the characteristic parameters for the problem under study areprovided. In Sec. 3 we briefly recall the key points of the HLBM and the fundamental equations. In Sec. 4the refill procedure for the lattice nodes in the proximity of the oscillating lamina is presented. In Sec. 5 wedescribe the numerical setup used for simulations. In Sec. 6 the results of this work are presented and discussed.Finally, in Sec. 7, we conduct a qualitative analysis of the numerical performance of the MG-HLBM and wediscuss some basic implementation aspects, while in Sec. 8 we summarize the key points of the work in theconclusions. 2 Problem statement
We consider the unsteady flow generated by a thin, rigid lamina undergoing transverse harmonic oscillations ina viscous, incompressible, quiescent fluid. The lamina is assumed to be infinitely long, such that the analysisis conducted within a two-dimensional framework. The cross-section of the lamina is a rectangle with length L and thickness T , with an aspect ratio T /L equal to 1/100. At each time t the vertical position of the lamina isgiven by a sinusoidal function y ( t ) = Asin ( ωt ), where A denotes the oscillation amplitude and ω is the radianfrequency.The specific problem of interest is governed by two parameters, that are the non-dimensional amplitude ofoscillation ε = A/L and the frequency parameter β = ρωL / (2 πµ ), respectively [33, 21, 41]. In particular, thenon-dimensional amplitude of oscillation is directly related to the KC number as κ = 2 πε .The Reynolds number based on the length of the lamina L and the maximum value of the lamina velocity, V max = Aω , can be expressed as a function of ε and β , as follows: Re = 2 πβε .In this work, we perform numerical simulations in the range of oscillatory Reynolds numbers Re = 2 . − . ε = [0.02, 0.03, 0.04, 0.05, 0.075, 0.1] and β = [20, 50, 100, 150,200, 250, 300]. In this regime, effects of fluid inertia are not negligible, and vorticity generated by the movingboundary leads to nonlinear inertial and damping effects [22, 21]. In particular, the set of control parameters( β , ε ) chosen for each simulation satisfies the following approximate correlation, provided by Aureli et al. [23]: β < . ε − . . (1)Relation (1) represents a reference boundary within which the force exerted by the fluid on the lamina isexpected to be a purely harmonic time function. For a cantilever beam moving in a viscous fluid under harmonic base excitation, the hydrodynamic load perunit length due to the motion of the fluid around the beam can be expressed, in the frequency domain, as [20]:ˆ F ( ω ) = F e iφ = π ρω L Θ( β, ε ) A (2)where the superimposed hat denotes a phasor quantity, and F and φ represent the force amplitude and thephase shift between the force and the displacement of an harmonic response, respectively. In equation (2),Θ( β, ε ) indicates a complex hydrodynamic function. By following [21, 23], such a quantity can be expressed asthe sum of two contributions: Θ( β, ε ) = Γ( β ) + ∆( β, ε ) . (3)The term Γ( β ) corresponds to the complex hydrodynamic function proposed in [20] for infinitely small amplitudeoscillations, while ∆( β, ε ) is the correction term for finite amplitude oscillations. For the range of controlparameters β and ε considered in this work, the following semianalytical formulations have been proven to wellapproximate the hydrodynamic function [23]:Γ( β ) = 1 .
02 + 2 . β − / − i . β − / (4)∆( β, ε ) = − i . β / ε (5)where i is the imaginary unit.To extract the hydrodynamic function components, for each simulation, that is, for each set of control parame-ters, we assume that the force response is purely harmonic. Therefore, in the same fashion of [23, 33], a singleharmonic sine model of the form F ( t ) = F sin( ωt + φ ) is used for the least square fitting of the force timehistory. After discarding the first cycle of oscillation, we consider three cycles of oscillations to identify thevalues of F and φ , while the radian frequency ω is provided in the fitting model as an input. From these twoquantities we then compute the real and the imaginary parts of Θ( β, ε ). The moving grid approach proposed in this work is based on a lattice Boltzmann method applied to hybridstructured and unstructured grids recently proposed by the authors (HLBM) [39]. For clarity of exposition, inthis section the key points of such a method are recalled.The HLBM combines the standard single-time relaxation lattice Boltzmann scheme with an unstructured finite-volume lattice Boltzmann formulation [42, 43, 44, 45, 46].The standard scheme is applied on a uniformly spaced lattice domain while the finite-volume approach is solved3n a unstructured mesh of triangular elements. The two mesh components overlap to each other, entirelycovering the computational domain.As far as the standard lattice Boltzmann method is concerned, this is based on the following equation: f i ( x + c i ∆ t s , t + ∆ t s ) − f i ( x , t ) = − ∆ t s τ s [ f i ( x , t ) − f eqi ( x , t )] , (6)where τ s and ∆ t s are the relaxation time and the time step, respectively, related to the ’structured’ scheme.In equation 6, f i ( x , t ) is the distribution function, representing the probability of finding a fluid particle atposition x and time t that is moving along the i -th lattice direction with a discrete speed c i . The equilibriumdistribution functions f eqi ( x , t ) are given by the following expression: f eqi ( x , t ) = w i ρ ( x , t ) (cid:26) c i · u ( x , t ) c s + [ c i · u ( x , t )] c s − [ u ( x , t )] c s (cid:27) (7)where c s is the lattice speed of sound, the parameters w i are a set of weights normalized to unity, and ρ ( x , t )and u ( x , t ) are the fluid density and velocity, respectively, which are given by the first two moments of thedistribution function. It can be proven that, for an ideal incompressible fluid flow, equation 6 reproduces theNavier-Stokes equations, when pressure is p = ρc s and kinematic viscosity is ν = c s ( τ s − ∆ t s ).The finite-volume lattice Boltzmann method is a cell-vertex type scheme which is expressed by the followingequation, for each node P of the mesh: f i ( P, t + ∆ t u ) = f i ( P, t ) + ∆ t u K X k =0 S ik f i ( P k , t )+ − ∆ t u τ u K X k =0 C ik [ f i ( P k , t ) − f eqi ( P k , t )] , (8)where τ u and ∆ t u are the relaxation time and the time step, respectively, related to the ’unstructured’ scheme[42, 43]. In equation 8 k = 0 denotes the pivotal node P and the summations run over the nodes P k connectedto P . The quantities S ik and C ik represent the streaming and collisional matrices of the i -th population relatedto the k -th node, respectively. The equilibrium distribution functions are defined by equation (7). For theunstructured lattice Boltzmann method the theoretical kinematic viscosity is ν = c s τ u [43].The value of the lattice speed of sound, c s , is the same in both the unstructured and the standard approaches.In this work, the D Q c s = 1 / √ n times, being ∆ t s = n ∆ t u . The exchange of infor-mation between the two mesh components, in terms of distribution functions, takes place at defined interpolationnodes, by applying the following set of equations:˜ f s i = f eq, u i + 2 ( τ s − ∆ t s )2 τ s − ∆ t s f neq, u i (9)˜ f u i = f eq, s i + (cid:18) − ∆ t s τ s (cid:19) f neq, s i + ∆ t u K X k =0 S ik (cid:2) f eq,?ik + f neq,?ik (cid:3) − ∆ t u τ u K X k =0 C ik f neq,?ik (10)where superscripts s and u refer to the ’structured’ and the ’unstructured’ nodes, respectively, and f neqi is thenon-equilibrium distribution function. The quantities f eq,?ik and f neq,?ik in the summation terms of equation (10)are defined as follows: f eq,?ik = f eq, s ik , f neq,?ik = (cid:18) − ∆ t s τ s (cid:19) f neq, s ik (11)or f eq,?ik = f eq, u ik , f neq,?ik = f neq, u ik (12)depending on whether the k -th node is an interpolation node or not, respectively. The distribution functions f eq, u i and f neq, u i of the right-hand side of equation (9) and f eq, s i and f neq, s i of the right-hand side of equation(10) are evaluated by interpolation procedure.Equations (9) and (10) represent the post-collision step for the interpolation nodes of structured and unstruc-tured mesh, respectively. Such relations are obtained by imposing consistency of viscosity in the flow field andcontinuity of velocity, density and stresses across the interface.4 Refill procedure
In the numerical approach proposed in this work, the unstructured mesh surrounding the lamina is enabled tomove within the structured fixed domain. At each time interval the hybrid method is applied and the flow fieldis computed all over the overlapping grid system. To accomplish this aim, we propose a refill procedure forthe unstructured lattice sites, which is performed at the beginning of each time interval. Such a procedure isbased on the attribution of both macroscopic quantities and distribution functions to the relocated unstructurednodes.Figure 1 illustrates a generic overlapping grid system at different times. In order to describe the procedure, wefirst recall the functional definition for the nodal points of each grid. According to the nomenclature adoptedin [39], nodes are classified as discretization , interpolation and unused . As far as the unstructured mesh isconcerned, all the boundary nodes are considered as interpolation nodes when this boundary represents aninterface with the structured mesh (gray-filled nodes in Figure 1). All of the interior unstructured nodes areconsidered as discretization nodes. Regarding the structured mesh, a fictitious curve, which is defined within theunstructured grid boundaries, separates the discretization nodes from the unused ones. Among the structured unused nodes, those with at least one discretization node as neighbor, in any of the allowed lattice directions,are redefined as interpolation nodes (dark gray-filled nodes in Figure 1).Let’s then consider the unstructured mesh at a given time t . After one time step ∆ t s , the unstructured nodesare relocated in a new position. Since the motion of the lamina is purely translational, all nodes move withsame velocity along the vertical direction. Therefore, unstructured nodes at time t + ∆ t s are classified as type1 , type 2 or type 3 depending on their location with respect to the mesh configuration at time t . In particular, type 1 refers to those nodes which are out of the unstructured grid at level time t and whose correspondingposition is located within the physical computational domain, namely in the interior of a structured elementdefined by four computed nodes. Type 2 is attributed to each node lying out of the unstructured grid at leveltime t and whose corresponding position is located inside the body, where no computation is performed. Thenodes in the new configuration whose corresponding position at level time t lies inside the unstructured grid isconsidered as type 3 . We remark that unstructured nodes defining the solid-wall of the body are not consideredin such a classification. (3)(2)(1)configuration at time t configuration at time t +∆ t s Figure 1: Schematic representation of a symmetric portion of the mesh configuration at two different timelevels: t (left) and t + ∆ t s (right). Dashed line is the fictitious line which is drawn in order to define thestructured interpolation nodes. Gray-filled nodes and dark gray-filled nodes represent the interpolation pointsfor the unstructured and the structured grid, respectively. The functional classification for the unstructurednodes is represented: type 1 , (1); type 2 , (2); type 3 , (3).To ensure consistency of the fluid flow, each type of unstructured node is initialized, in the new position x + ∆ x , at the beginning of each structured time step, according to a specific rule.For nodes of type 1 , values of non-equilibrium distribution functions f neq, u i are computed by means of thefollowing equation: f neq, u i = (cid:18) − ∆ t s τ s (cid:19) f neq, s i . (13)where the non-equilibrium distribution functions f neq, s i are evaluated via a bilinear Lagrange interpolationperformed over the four donor nodes constituting the structured element which encloses the unstructured nodeunder consideration. The structured donor nodes are indeed interpolation or discretized nodes, thus their5ontributions are known. Moreover, the macroscopic variables are computed by applying the same interpolationscheme, thus allowing the reconstruction of the equilibrium distribution functions.As far as unstructured nodes of type 2 are concerned, at the beginning of each structured time step we setthe same macroscopic variables and distribution functions determined for the same nodes in the previous meshconfiguration. This approximation is valid as the maximum value for the velocity of the lamina is much smallerthan a lattice unit (in this case is of the order of 0.01 in lattice units).For nodes of type 3 , velocity components, density and distribution functions f u i are evaluated via barycentricinterpolation by considering the three vertexes of the triangular element in the previous mesh configurationenclosing the same node at new position. Once macroscopic variables are computed, it is then possible todetermine also the equilibrium distribution functions.Regarding unstructured nodes defining the solid-wall of the body, we set the velocity as a boundary conditionand we impose the same density and the same distribution functions f u i of the same node in previous meshconfiguration.In addition, a dynamic redefinition of the structured interpolation nodes is performed at each time iteration.In fact, the fictitious line designed to identify such nodes moves jointly with the unstructured grid. Therefore,for each structured interpolation node overcome by the moving fictitious line, a new one must be defined.In order to test the refill procedure, several unstructured meshes involving all of the three described typesof node have been considered. However, we emphasize that the refinement level of the unstructured meshesemployed for simulating the fluid flow problem under consideration is such that the minimum distance betweentwo nodes is greater than the distance covered by any node within one time interval. Therefore, in practice,nodes of type 2 are never defined in this particular case. The computational domain is constituted by a square region inside which the lamina is defined. No-slip wallboundary conditions are imposed at the lateral sides of the domain, while on top and bottom boundaries weset outflow conditions.In order to ensure an adequate representation of the flow field, we conduct a sensitivity analysis to select thecomputational domain size. Specifically, for a fixed pair of the parameters ( β , ε ), the following square domainsizes are tested: (1.5 × L , (3 × L , (4 × L , (5 × L and (10 × L . The length of the lamina is kept toa constant value. Figure 2 shows the time history of the resultant (non-dimensional) force exerted by the fluidon the lamina. 6 4000 8000 12000 16000-4-2024 Time-step F o r ce [l. u .] (1 . × . L (3 × L (4 × L (5 × L (10 × L Figure 2: Hydrodynamic force time history for selected domain sizes: a) (1.5 × L , black dashed line; b)(3 × L , blue dotted line; c) (4 × L , violet dash dotted line; d) (5 × L , red solid line; e) (10 × L , orangedensely dotted line.Although the numerical discrepancies observed by inspecting Figure 2 are small, the results, in terms ofhydrodynamic function components vary significantly among the cases. In Table 1 we report the percentageerror related to the case of domain size equal to (5 × L . The convergence study shows that the differencesDomain size Err % Re [Θ( β, ε )] Err % − Im [Θ( β, ε )](1.5 × L × L × L × L × L β = 100 and ε = 0 .
1. The relative error represents the deviation, in percent, from the valuesdetermined with the domain size equal to (5 × L .between the cases (5 × L and (10 × L can be considered negligible. Therefore, we select a final domainwith size (5 × L .In order to further limit the computational cost, the structured domain is subdivided in four refinement levels.The length of the lamina is fixed at 800 lattice units for all the simulations, while the relaxation time at thefinest structured level is varied between 0.7 and 0.9. For this set of parameters, the maximum value of thelamina speed results to be about 0.025 lattice units, which is much below the compressibility limit. In Figure 3 we report the numerical findings for the real and the imaginary components of the hydrodynamicfunction, in comparison with the semyanalitical expression (equation 3), with results obtained by Falcucci etal. [33] via standard LB approach, and with those from De Rosis and Lévêque [36] , which use a combinedIB-LBM. The results are here reported as a function of the frequency parameter β .7 50 100 150 200 250 30000 . . . β − I m [ Θ ] , R e [ Θ ] (a) . . . β − I m [ Θ ] , R e [ Θ ] (b) . . . β − I m [ Θ ] , R e [ Θ ] (c) . . . β − I m [ Θ ] , R e [ Θ ] (d) . . . β − I m [ Θ ] , R e [ Θ ] (e) . . . . β − I m [ Θ ] , R e [ Θ ] (f) Figure 3: Real part Re [Θ( β, ε )] (gray) and imaginary part − Im [Θ( β, ε )] (black) of the complex hydrodynamicfunction. The results determined via MG-HLBM (circles) are compared with the semianalytical formula 3provided by Aureli et al. [23] (lines), results from De Rosis and Lévêque [36] (triangles, panels (d),(e),(f)), andthe numerical findings from Falcucci et al. [33] (squares, panels (a),(d),(e),(f)). (a): ε = 0 .
02; (b): ε = 0 . ε = 0 .
04; (d): ε = 0 .
05; (e): ε = 0 . ε = 0 . Re [Θ( β, ε )] is nearly constant with the variations of parameter ε , whilethe imaginary part − Im [Θ( β, ε )] increases whit ε for each fixed value of β , as suggested by the reference8iterature [23].Our results are in line with those from the semianalytical expression and with results from [36]. We observe thatthe general dependence of viscous damping and added mass on the control parameters is properly estimatedby the present MG-HLBM. On the other hand, Figure 3 show significant differences between the MG-HLBMresults and those obtained via standard LBM, which appear to be overestimated. In spite of an increasing ofthe level of complexity introduced by the presence of the unstructured grid, the MG-HLBM results to be moreaccurate than the standard LBM in the prediction of the hydrodynamic load for the case under consideration.A crucial aspect is represented by the direct computation of the forces at the real boundary of the solid body,which is not allowed in the standard LBM. Also, the presence of an unstructured grid allows an high and flexiblelevel of refinement in proximity of the solid body, thus leading to a reduction of the computational cost withrespect to the case of uniform grid, as demonstrated in [39, 40].For a further qualitative comparison between MG-HLBM and literature results, we report a complete overviewof experimental and numerical findings in aggregated form in Figures 4 and 5, for a broad range of controlparameters. The real and the imaginary parts of the hydrodynamic function are reported as functions of thenon-dimensional amplitude of oscillation ε .0.0001 0.001 0.01 0.1 1 1001234567 ε R e [ Θ ] ε R e [ Θ ] Figure 4: Aggregated results for the real part Re [Θ( β, ε )] of the complex hydrodynamic function (color online).Black circles: present study; gray pentagons: [23]; magenta stars: [33]; red squared: [21]; blue asterisks: [22];orange triangles: [36]; violet diamonds: [47]; green crosshairs: [24]; brown crosses: [48, 49]9.0001 0.001 0.01 0.1 1 10012345678 ε − I m [ Θ ] ε − I m [ Θ ] Figure 5: Aggregated results for the imaginary part − Im [Θ( β, ε )] of the complex hydrodynamic function (coloronline). Black circles: present study; gray pentagons: [23]; magenta stars: [33]; red squared: [21]; blue asterisks:[22]; orange triangles: [36]; violet diamonds: [47]; green crosshairs: [24]; brown crosses: [48, 49].We observe a good agreement with the considered literature, with the results obtained by the proposedmethod overlapping those available from previous works. In fact, the values for both components of the hydro-dynamic function cast within similar ranges of those obtained by other numerical and experimental methods.In particular, our results suggest that variations of added mass effects are negligible for different values of theamplitude of oscillations. This is in contrast with findings obtained via standard LBM [33]. Moreover, consis-tently with the results reported in the scientific literature, the dependence of the imaginary part on β tends toreduce as ε is increased.Finally, in Figure 6 we present the evolution of velocity and pressure fields, respectively, for a representativecase of the fluid dynamic problem under analysis ( β = 100, ε = 0 . a) (b) (c)(d) (e) (f) Figure 6: Detail of instantaneous velocity (a-c) and pressure (d-f) fields for β = 100 and ε = 0 . β and ε within the ranges analyzed. In this section, we analyze the numerical performance of the moving grid approach for the hybrid latticeBoltzmann method and, finally, we discuss some implementation aspects.To asses the numerical efficiency of the MG-HLBM, we evaluate the computational cost required by the proposedmethod in comparison with the standard LBM. To this aim, we consider the simplified case of an hybrid meshcomposed by the overlapping of a uniformly spaced structured region, which covers the entire computationaldomain, with the same unstructured mesh previously used in this study. In other words, we do not considerlocal refinements for the structured mesh. In such a scenario and for the case of static unstructured mesh, wecan apply the following equation, introduced in [39], to quantify the differences, in terms of computational cost,11etween hybrid and standard LB approaches: G = T h T eq = 1 + χs ng . (14)In equation 14, G is the speedup (gain) indicator, which is defined as the ratio of the wall-clock time per algorithmiteration associated to the hybrid method ( T h ) over the wall-clock time per algorithm iteration with referenceto an equivalent LBM implementation on uniform spaced mesh ( T eq ), namely a standard LBM with similarresolution around the solid body. A value of G less than 1 indicates that, by definition, the computational cost ofthe hybrid method is lower than that of standard LBM, when a similar resolution around the body is considered.The speedup indicator G depends on four parameters which are, respectively: the amount of unstructured nodeswith respect to the number of structured ones ( χ ), the ratio between the number of structured nodes whichwould be required by an equivalent LBM and the number of structured nodes actually used within the hybridmesh ( g ), the ratio of the standard LBM processing speed to the specific processing speed associated to theunstructured method ( s ), and the number of sub-iterations required by the unstructured method ( n ). Forthe present case, the percentage number of unstructured nodes χ is equal to 0.244%. The value of χ refers tothe specific mesh under consideration. However, such a value is representative of a typical unstructured meshrequired by the hybrid method and it is in line with those reported also in [39, 40]. In fact, we stress that thehybrid method requires only a small portion of the computational domain to be covered by the unstructuredmesh, namely the region around the solid body. This represents one of the key elements of the proposedMG-HLBM. Further, n is set to 50. We emphasize here that the hybrid mesh adopted in this work has beendesigned on the basis of either a sensitivity analysis, as presented in Sec. 5, and a stability analysis. As faras stability aspects are concerned, we recall that, due to the explicit time integration, both the standard LBand the unstructured method must satisfy the following condition: ∆ t < τ . While such a condition is alwaysverified for the standard scheme, this is not necessarily the case for the unstructured one. This represents aclear limitation for the present hybrid method since it constrains, eventually, the value of n , namely the numberof sub-iterations required by the unstructured method in order to synchronize the solution. The parameter s was evaluated to be roughly equal to 4. Such an average value was obtained by actually measuring the codecomputational performance for several unstructured meshes. In addition, the number of nodes required by theequivalent standard LBM to have a similar resolution around the lamina is such that g ’ .
78. With this setof parameters, we find G ’ .
84. This result indicates that, at each time-step, the hybrid method providesa gain in efficiency with respect to the standard LBM, when a similar resolution around the solid body isconsidered. Notwithstanding a relatively low speedup, we emphasize that the code developed in this work stillrequires further optimization. Therefore, the value of G obtained for the specific case under analysis should beconsidered as an illustrative example.To complete the assessment of the MG-HLBM numerical performance we must include in the analysis theadditional computational cost required by the mesh-motion and refill routine. In particular, such a routineis performed, at each time-step, by means of three sequential steps, that are: 1) motion of the unstructuredmesh, 2) re-classification of interpolation nodes and donor elements, on structured and unstructured mesh,respectively, 3) refill procedure. Despite this routine being laborious for the general case of arbitrary/inducedmotion of the unstructured mesh, its computational realization becomes particularly efficient when the motionof the unstructured mesh is imposed, such as in the present case, since it is known a-priori the location of theunstructured nodes at each time-step. Another important aspect is related to the number of unstructured nodes.As observed also in [39, 40], the number of unstructured nodes required by the hybrid method is significantlylower than the number of structured nodes used within the whole computational domain. For instance, in thepresent case, the number of unstructured nodes is less than 1% of the total number of structured nodes. Thisimplies that, the interpolations needed to be performed within the refill procedure involve only a relatively lowamount of computer operations. In order to quantify the influence, in terms of computational performance, ofthe mesh-motion and refill procedure, we measured the computational time associated with such a routine. Asa result, we found out that its processing speed is about two times higher than the processing speed associatedto the unstructured method. Therefore, we can take into account of this additional computational cost byconsidering a modified number of unstructured sub-iterations n in equation 14, such that: n = n + 0.5. Theoverall computed value of the speedup indicator results then to be fairly the same of the one computed forthe static case ( G ’ . Conclusion
In this work, we extend the HLBM toward a moving grids approach and we explore its capability through theparametric study of a thin lamina undergoing transverse harmonic oscillations in a viscous quiescent fluid. Thenumerical findings are in line with those presented in literature and obtained by different numerical strategiesas well as experimental approaches. The performed analysis demonstrates that the proposed method is ableto properly predict the main features of the fluid flow induced by the motion of the lamina. In particular, weobserve a significant improvement, in terms of accuracy, with respect to previous works performed via standardlattice Boltzmann method on cartesian grid. The proposed method represents a viable alternative also to otherapproaches such as the IB-LBM. In this regard, the MG-HLBM presents some key advantage. First of all,the hybrid strategy allows the possibility of applying an efficient local refinement. Second, the presence of anunstructured body-fitted mesh allows the direct evaluation of the forces at the real boundary of the body, thusleading to accurate results while preserving the computational efficiency of standard LBM. These features makethe MG-HLBM particularly appealing for multi-scale problems. On the contrary, the IB-LBM is based on asubstantially different concept, therefore the choice between the use of this or the present method should bebased, in general, on the peculiarity of the problem under analysis. Also, the MG-HLBM and, in general, thehybrid method, can be regarded as a complex boundary condition for solid walls, which can be extended todeformable objects. This can be the subject of future works. These aspects are instrumental when dealing withfluid-structure interaction problems, involving complex geometries.
This work was supported by the Italian Ministry of Education, University and Research under PRIN grantNo. 20154EHYW9 "Combined numerical and experimental methodology for fluid structure interaction in freesurface flows under impulsive loading".The numerical simulations were performed on
Zeus
HPC facility, at the University of Naples "Parthenope";
Zeus
HPC has been realized through the Italian Government Grant PAC01_00119
MITO - Informazioni Multime-diali per Oggetti Territoriali , with Prof. Elio Jannelli as the Scientific Responsible.One of the authors (S. Succi) acknowledges funding from the European Research Council under the Euro-pean Union’s Horizon 2020 Framework Programme (No. FP/2014-2020) ERC Grant Agreement No. 739964(COPMAT). 13 eferences [1] A. Maali, C. Hurth, R. Boisgard, C. Jai, T. Cohen-Bouhacina, J. Aimè, Hydrodynamics of oscillatingatomic force microscopy cantilevers in viscous fluids, Journal of Applied Physics 97 (2005) 074907.[2] S. Kirstein, M. Mertesdorf, M. Schoenhoff, The influence of a viscous fluid on the vibration dynamics ofscanning near-field optical microscopy fiber probes and atomic force microscopy cantilevers, Journal ofApplied Physics 84 (1998) 1782–1790.[3] M. Kimber, S. Garimella, A. Raman, Local heat transfer coefficients induced by piezoelectrically actuatedvibrating cantilevers, Transactions of the ASME Journal of Heat Transfer 129 (2007) 1168–1176.[4] M. Kimber, R. Lonergan, S. Garimella, Experimental study of aerodynamic damping in arrays of vibratingcantilevers, Journal of Fluids and Structures 5 (2009) 1334–1347.[5] C. Castille, I. Dufour, C. Lucat, Longitudinal vibration mode of piezoelectric thick-film cantilever-basedsensors in liquid media, Applied Physics Letters 96 (2010) 154102. doi:10.1063/1.3387753 .[6] H. Hosaka, K. Itao, S. Kuroda, Damping characteristics of beam-shaped micro-oscillators, Sensors andActuators A: Physical 49 (1995) 87–95.[7] M. Kimber, S. Garimella, Measurement and prediction of the cooling characteristics of a generalized vi-brating piezoelectric fan, International Journal of Heat and Mass Transfer 52 (2009) 4470–4478.[8] M. Aureli, V. Kopman, M. Porfiri, Free-Locomotion of Underwater Vehicles Actuated by Ionic PolymerMetal Composites, IEEE/ASME Transactions on Mechatronics 15 (2010) 603–614.[9] K. Abdelnour, E. Mancia, S. Peterson, M. Porfiri, Hydrodynamics of underwater propulsors based on ionicpolymer metal composites: a numerical study, Smart Materials and Structures 18 (2009) 085006.[10] C. Pooley, J. Yeomans, Lattice Boltzmann simulation techniques for simulating microscopic swimmers,Computer Physics Communications 179 (2008) 159–164.[11] R. Batra, M. Porfiri, D. Spinello, Electromechanical model of electrically actuated narrow microbeams,Journal of Microelectromechanical Systems 15 (2006) 1175–1189.[12] R. Batra, M. Porfiri, D. Spinello, Review of modeling electrostatically actuated microelectromechanicalsystems, Smart Materials and Structures 16 (2007) R23.[13] Y. Cha, L. Shen, M. Porfiri, Energy harvesting from underwater torsional vibrations of a patterned ionicpolymer metal composite, Smart Materials and Structures 22 (2013) 055027.[14] M. Aureli, C. Prince, M. Porfiri, S. Peterson, Energy harvesting from base excitation of ionic polymer metalcomposites in fluid environments, Smart Materials and Structures 19 (2010) 015003.[15] E. Tuck, Calculation of unsteady flows due to small motions of cylinders in a viscous fluid, Journal ofEngineering Mathematics 3 (1969) 29–44.[16] A. Facci, M. Porfiri, Nonlinear hydrodynamic damping of sharp-edged cantilevers in viscous fluids under-going multi-harmonic base excitation, Journal of Applied Physics 112 (2012) 124908.[17] G. Stokes, On the effect of the internal friction of fluids on the motion of pendulums, Transactions of theCambridge Philosophical Society 9 (1851) 1–141.[18] C. Green, J. Sader, Small amplitude oscillations of a thin beam immersed in a viscous fluid near a solidsurface, Physics of Fluids 17 (2005) 073102.[19] G. Di Ilio, I. Sahin, A. Tafuni, Unsteady Stokes Flow for a Vibrating Cantilever Under a Free-Surface,Proceedings of the ASME 2014 International Mechanical Engineering Congress and Exposition 9 (2014)IMECE2014–36929. doi:10.1115/IMECE2014-36929 .[20] J. Sader, Frequency response of cantilever beams immersed in viscous fluids with applications to the atomicforce microscope, Journal of Applied Physics 84 (1998) 64–76.[21] M. Aureli, M. Porfiri, Low frequency and large amplitude oscillations of cantilevers in viscous fluids, Appl.Phys. Lett. 96 (2010) 164102. doi:10.1063/1.3405720 .1422] R. Bidkar, M. Kimber, A. Raman, A. Bajaj, S. Garimella, Nonlinear aerodynamic damping of sharp-edged flexible beams oscillating at low Keulegan–Carpenter numbers, J. Fluid Mech. 634 (2009) 269–289. doi:10.1017/S0022112009007228 .[23] M. Aureli, M. Basaran, M. Porfiri, Nonlinear finite amplitude vibrations of sharp-edged beams in viscousfluids, Journal of sounds and vibrations 331 (2012) 1624–1654. doi:10.1016/j.jsv.2011.12.007 .[24] A. Tafuni, I. Sahin, Non-linear hydrodynamics of thin laminae undergoing large harmonic oscillations in aviscous fluid, Journal of Fluids and Structures 52 (2015) 101–117. doi:10.1016/j.jfluidstructs.2014.10.004 .[25] A. Tafuni, I. Sahin, Hydrodynamic loads on vibrating cantilevers under a free surface in viscous fluidswith SPH, Proceedings of the ASME 2013 International Mechanical Engineering Congress and Exposition7 (2013) IMECE2013–63792. doi:10.1115/IMECE2013-63792 .[26] S. Ansumali, I. Karlin, H. Öttinger, Thermodynamic theory of incompressible hydrodynamics, Phys. Rev.Lett. 94 (2005) 080602.[27] S. Succi, The Lattice Boltzmann Equation for Fluid Dynamics and Beyond, Clarendon Press, 2001.[28] S. Succi, The Lattice Boltzmann Equation for Complex States of Flowing Matter, Oxford University Press,2018.[29] Y. Qian, D. D’Humières, P. Lallemand, Lattice BGK models for Navier-Stokes equation, Europhys. Lett.17 (1992) 479–484.[30] R. Benzi, S. Succi, M. Vergassola, The Lattice Boltzmann Equation: Theory and Applications, Phys. Rep.222 (1992) 145–197. doi:10.1016/0370-1573(92)90090-M .[31] F. Higuera, S. Succi, R. Benzi, Lattice gas dynamics with enhanced collisions, Europhys. Lett. 9 (1989)345–349.[32] G. Di Ilio, B. Dorschner, G. Bella, S. Succi, I. Karlin, Simulation of turbulent flows with the entropicmultirelaxation time lattice Boltzmann method on body-fitted meshes, J. Fluid Mech. 849 (2018) 35–56. doi:10.1017/jfm.2018.413 .[33] G. Falcucci, M. Aureli, S. Ubertini, M. Porfiri, Transverse harmonic oscillations of laminae in viscous fluids:a lattice Boltzmann study, Phil. Trans. R. Soc. A 369 (2011) 2456–2466. doi:10.1098/rsta.2011.0062 .[34] A. D. Rosis, G. Falcucci, S. Ubertini, F. Ubertini, S. Succi, Lattice Boltzmann analysis of fluid-structureinteraction with moving boundaries, Communications in Computational Physics 13 (2013) 823–834.[35] A. D. Rosis, G. Falcucci, S. Ubertini, F. Ubertini, Aeroelastic study of flexible flapping wings by a cou-pled lattice Boltzmann-finite element approach with immersed boundary method, Journal of Fluids andStructures 49 (2014) 516–533.[36] A. D. Rosis, E. Lévêque, Harmonic oscillations of a thin lamina in a quiescent viscous fluid: A numericalinvestigation within the framework of the lattice Boltzmann method, Computers and Structures 157 (2015)209–217.[37] Y. Shi, J. Sader, Lattice Boltzmann method for oscillatory Stokes flow with applications to micro- andnanodevices, Phys. Rev. E 81 (2010) 036706. doi:10.1103/PhysRevE.81.036706 .[38] C. Colosqui, D. Karabacak, K. Ekinci, V. Yakhot, Lattice Boltzmann simulation of electromechanicalresonators in gaseous media, J. Fluid Mech. 652 (2010) 241–257. doi:10.1017/S0022112010000042 .[39] G. Di Ilio, D. Chiappini, S. Ubertini, G. Bella, S. Succi, Hybrid lattice Boltzmann method on overlappinggrids, Phys. Rev. E 95 (2017) 013309. doi:10.1103/PhysRevE.95.013309 .[40] G. Di Ilio, D. Chiappini, S. Ubertini, G. Bella, S. Succi, Fluid flow around NACA 0012 airfoil at low-Reynolds numbers with hybrid lattice Boltzmann method, Computers and Fluids 166 (2018) 200–208. doi:10.1016/j.compfluid.2018.02.014 .[41] C. Wang, On high-frequency oscillatory viscous flows, J. Fluid Mech. 32 (1968) 55–68. doi:10.1017/S0022112068000583 .[42] S. Ubertini, G. Bella, S. Succi, Lattice Boltzmann method on unstructured grids: Further developments,Phys. Rev. E 68 (2003) 016701. doi:10.1103/PhysRevE.68.016701 .1543] S. Ubertini, S. Succi, G. Bella, Lattice Boltzmann schemes without coordinates, Phil. Trans. R. Soc. Lond.A 362 (2004) 1763–1771. doi:10.1098/rsta.2004.1413 .[44] S. Ubertini, S. Succi, A Generalised Lattice Boltzmann Equation on Unstructured Grids, Commun. Com-put. Phys. 3 (2007) 342–356.[45] A. Zarghami, C. Biscarini, S. Succi, S. Ubertini, Hydrodynamics in Porous Media: A Finite Volume LatticeBoltzmann Study, Journal of Scientific Computing 59 (2014) 80–103.[46] A. Zarghami, S. D. Francesco, C. Biscarini, Porous substrate effects on thermal flows through a REV-scalefinite volume lattice Boltzmann model, International Journal of Modern Physics C 25 (2014) 1350086.[47] M. Jalalisendi, R. Panciroli, Y. Cha, M. Porfiri, A particle image velocimetry study of the flow physicsgenerated by a thin lamina oscillating in a viscous fluid, Journal of Applied Physics 115 (2014) 054901. doi:10.1063/1.4863721 .[48] J. Graham, The forces on sharp-edged cylinders in oscillatory flow at low Keulegan–Carpenter numbers,Journal of Fluid Mechanics 97 (1980) 331–346. doi:10.1017/S0022112080002595doi:10.1017/S0022112080002595