A weakly compressible hybridizable discontinuous Galerkin formulation for fluid-structure interaction problems
Andrea La Spina, Martin Kronbichler, Matteo Giacomini, Wolfgang A. Wall, Antonio Huerta
AA weakly compressiblehybridizable discontinuous Galerkin formulationfor fluid-structure interaction problems
Andrea La Spina , , ∗ , Martin Kronbichler , Matteo Giacomini ,Wolfgang A. Wall and Antonio Huerta Lehrstuhl f¨ur Numerische Mechanik (LNM), Technische Universit¨at M¨unchen (TUM),Garching b. M¨unchen, Germany. Laboratori de Calcul Numeric (LaC`aN), ETS de Ingenieros de Caminos, Canales y Puertos,Universitat Polit`ecnica de Catalunya, Barcelona, Spain. Centre Internacional de M`etodes Num`erics a l’Enginyeria (CIMNE), Barcelona, Spain.
Abstract
A scheme for the solution of fluid-structure interaction (FSI) problems withweakly compressible flows is proposed in this work. A novel hybridizable discontinu-ous Galerkin (HDG) method is derived for the discretization of the fluid equations,while the standard continuous Galerkin (CG) approach is adopted for the structuralproblem. The chosen HDG solver combines robustness of discontinuous Galerkin(DG) approaches in advection-dominated flows with higher order accuracy and effi-cient implementations. Two coupling strategies are examined in this contribution,namely a partitioned Dirichlet–Neumann scheme in the context of hybrid HDG-CGdiscretizations and a monolithic approach based on Nitsche’s method, exploiting thedefinition of the numerical flux and the trace of the solution to impose the couplingconditions. Numerical experiments show optimal convergence of the HDG and CGprimal and mixed variables and superconvergence of the postprocessed fluid velocity.The robustness and the efficiency of the proposed weakly compressible formulation,in comparison to a fully incompressible one, are also highlighted on a selection of twoand three dimensional FSI benchmark problems.
Keywords: fluid-structure interaction; finite elements; hybridizable discontin-uous Galerkin; weakly compressible flows; Navier–Stokes equations; Nitsche’smethod. ∗ Corresponding author: Andrea La Spina. E-mail: [email protected] 1 of 49 a r X i v : . [ phy s i c s . f l u - dyn ] S e p Introduction
The simulation of the interaction of fluid flows with flexible structures is of great interest inmany engineering fields and it has been extensively investigated over the last decades. Thenumerical techniques developed for the solution of this challenging multiphysics problemcan be catalogued based on many different aspects, for instance, with respect to the spatialdiscretization, the kinematical description, and the coupling of the fluid and the structuresubproblems.Among many other techniques developed so far, the finite element method is one ofthe most successful spatial discretization approaches for the solution of the partial differ-ential equations (PDEs) underlying many physical phenomena, including fluid-structureinteraction. The standard continuous Galerkin method provides computationally efficientdiscretizations with a very limited number of degrees of freedom (DOFs) for the solutionof elasticity problems. For an overview on CG methods for solid mechanics, the inter-ested reader is referred to [54]. On the other hand, the interest in discontinuous Galerkinmethods has increased over the last decades in the computational fluid dynamics commu-nity [12, 29, 43] because of their distinctive properties, such as the inherited stabilization ofthe convection terms in conservation laws, the ability to construct high order discretiza-tions on unstructured meshes and the flexibility in performing p -adaptivity in additionto the classical h -adaptivity. More recently, hybridizable discontinuous Galerkin methodshave gained a lot of attention owing to their reduced computational costs with respect toclassical matrix-based DG approaches, thanks to the reduced number of global DOFs in theassociated linear systems, especially for high-degree polynomial approximations. Moreover,the possibility to obtain a superconvergent solution through an efficient element-by-elementpostprocessing allows to obtain an improved approximation of the solution and to driveefficient degree adaptivity procedures [19,20,46]. In the context of flow problems, the HDGmethod has been successfully applied for the discretization of fully compressible flows [42]as well as incompressible flows [17, 19, 40]. The strong enforcement of the symmetry of thestress tensor via Voigt notation to retrieve the optimal convergence of the mixed variableand to ensure the superconvergence of the postprocessed solution without additional en-richment of the discrete spaces has been proposed in [17, 45] and it is exploited also in theproposed formulation for weakly compressible flows. The solution of FSI problems withincompressible flows by means of the HDG method for both the fluid and the structure hasbeen formulated in [47, 48]. However, these formulations are computationally much moreexpensive than the one proposed here and they moreover fail to provide an optimal con-vergent structural strain field (and therefore a superconvergent displacement field), losingtherefore one of the key advantages of the HDG method.A successful coupling of DG and CG methods for the solution of fluid-structure inter-action problems has been proposed in [15], where the high-order accuracy in time givenby the implicit-explicit Runge–Kutta method has been demonstrated on a non-trivial testproblem. As opposed to [15], the formulation proposed here aims to couple the HDG2ethod for the discretization of the fluid equations and the CG method for the solution ofthe structural problem and a high-order accuracy is demonstrated for the hybrid spatialdiscretization on a problem with manufactured solution. A first attempt to couple HDGand CG discretizations has been proposed in [41], while an improved minimally-intrusiveHDG-CG coupling has been formulated in [33] for the solution of multi-material structuralproblems, involving compressible and nearly incompressible solids.An important aspect in the simulation of multiphysics problems is the choice of anappropriate kinematical description. Pure fluid problems are commonly solved with anEulerian description, i.e., the computational mesh is fixed and the fluid moves with re-spect to the grid. This approach facilitates the treatment of large distortion in the fluidmotion and it is in particular indispensable in case of turbulent flows. On the contrary,the Lagrangian description is usually used in structural mechanics and with this approachthe nodes of the mesh follow the associated material particles during the motion, facili-tating the tracking of free surfaces and interfaces between different materials. ArbitraryLagrangian–Eulerian (ALE) algorithms [27] aim to combine the advantages of the classicalkinematical descriptions by introducing a computational mesh which can move with a ve-locity independent of the velocity of the material particles. This technique is particularlyuseful for flow problems in the presence of mobile and deforming boundaries, as it happensin fluid-structure interaction [34, 44, 53].In the formulation derived here, the Arbitrary Lagrangian–Eulerian method is used forthe description of the fluid flow, while the total Lagrangian method is adopted for thedescription of the structural motion.Regarding the coupling of the single fields, partitioned schemes solve one subproblemper time and exchange the interface information between the fluid and the structure. Theexchange of the interface state is performed just once per time step in the so-called loosely-coupled staggered approaches [11], while the solution of the flow and structural problemsare repeated within one time step until a convergence criterion is satisfied in the strongly-coupled staggered approaches [31]. The partitioned schemes may suffer several stabilityand convergence issues, but on the other hand they allow the use of well established andoptimized single-field solvers. In particular, these schemes are affected by a detrimentalphenomenon, defined in literature “artificial added mass effect” and analyzed in [3, 14]in the context of incompressible flows. In a recent work [32], it is analytically demon-strated how the introduction of a weak compressibility in the fluid formulation alleviatesthe constraints of the instability condition of the artificial added mass effect, thanks tothe reduction of the maximal eigenvalue of the so-called added mass operator. Moreover,in comparison to a fully incompressible solver, a significant reduction of the coupling it-erations and the computational time is observed. It is worth highlighting that alternativeapproaches to relax the incompressibility constraint based on the artificial compressibility,as the ones proposed in [1, 2] in the context of DG methods, provide a strongly consis-tent approximation of the incompressible Navier–Stokes equations. Embedding such fluidsolvers in a partitioned FSI code would not offer any beneficial contribution against the3rtificial added mass effect. A general framework for constructing high-order partitionedsolvers based on implicit-explicit Runge–Kutta methods for the solution of multiphysicsproblems (including fluid-structure interaction) has been introduced in [25], where fourconsistent predictors are proposed, leading to different partitioned solvers that preservethe theoretical order of accuracy of the temporal integration scheme. As opposed to thepartitioned strategies, in a monolithic framework [23,34], a unique solver is in charge of thesolution of the complete system of nonlinear equations. This approach exhibits a higher ro-bustness and it is usually faster compared to partitioned strategies in challenging cases, butit requires an ad hoc implementation. In addition, the adoption of efficient preconditionersis often needed for the solution of computationally demanding problems.In the present contribution, a strongly-coupled staggered scheme based on the Dirichlet–Neumann partitioning is revisited for the hybrid HDG-CG coupling and a novel FSI mono-lithic approach based on Nitsche’s method is also introduced. It is worth noting that thepartitioned algorithm adopted here considers a strong Gauss–Seidel-type predictor, whichis acknowledged in [25] to be most stable among the techniques therein analyzed. A weakcompressibility is then considered in the fluid formulation in order to provide an improvedrobustness and efficiency for the coupled solver.The present article is organized as follows. First, the novel HDG formulation forweakly compressible flows is derived in section 2 for pure fluid problems with an Ar-bitrary Lagrangian–Eulerian description. In section 3 the standard CG formulation fornonlinear elastodynamics is briefly recalled and the two HDG-CG coupling strategies forfluid-structure interaction are presented. Section 4 is devoted to the numerical validationof the pure fluid formulation first and the coupling strategies for fluid-structure interactionafterwards. Finally, in section 5, the results of this work are summarized. In this section, the governing equations of unsteady weakly compressible flows are firstpresented with regards to fixed domains and then formulated according to the ArbitraryLagrangian–Eulerian description to deal with moving domains. A brief overview on Voigtnotation is provided to handle symmetric tensors and the HDG formulation of the localand global problems is derived together with a postprocessing procedure to construct animproved approximation of the solution.
Let Ω x ∈ R n sd be a fixed open bounded domain with boundary ∂ Ω x = Γ Dx ∪ Γ Nx withΓ Dx ∩ Γ Nx = ∅ and n sd being the number of spatial dimensions and let T end > ∂ρ∂t + ∇ · ( ρ υ ) = 0 in Ω x × (0 , T end ) ,∂ρ υ ∂t + ∇ · ( ρ υ ⊗ υ ) − ∇ · σ = ρ b in Ω x × (0 , T end ) ,p ( ρ ) = 0 in Ω x × (0 , T end ) , (1)where ρ represents the fluid density, p the pressure, υ the velocity field, σ the Cauchystress tensor and b an applied body force. For a Newtonian fluid, it is assumed that thestress tensor and the augmented strain rate tensor are linearly related, therefore σ = − p I n sd + 2 µ ∇ S υ + λ ( ∇ · υ ) I n sd . (2)The operator ∇ S := (cid:0) ∇ + ∇ T (cid:1) returns the symmetric part of the gradient while I n sd denotes the n sd × n sd identity matrix. In equation (2), µ is the dynamic viscosity and λ theso-called second coefficient of viscosity. Stokes’ hypothesis states the following relationshipbetween the two material variables: λ = − µ. (3)The conservation of energy is taken into account by an equation of state p ( ρ ) = 0 defining arelationship between the fluid pressure and the density. For weakly compressible Newtonianfluids, a linear relationship is used [51]: ρ = ρ + ε ( p − p ) , (4)where ε is a (small) constant isothermal compressibility coefficient, while ρ denotes themass density evaluated at the reference pressure p . Equation (4) has been considered fornumerical simulations of weakly compressible flows in long tubes, such as waxy crude oil [52]and polymer extrusion [50]. Though equation (4) is very attractive for its simplicity, otherpressure-density relations have been used in literature, like the Murnaghan–Tait model [32].In order to derive the Arbitrary Lagrangian–Eulerian description of the flow, the ALEconvective velocity, defined as the velocity of the fluid relative to the moving backgroundmesh (whose velocity is indicated here with a ), is introduced: c = υ − a . (5)Using the ALE time derivative (i.e., the time derivative with respect to the referenceconfiguration), the governing equations of the fluid problem under analysis on a deforming5omain Ω can be written as: ∂ρ∂t + ρ ∇ · a + ∇ · ( ρ c ) = 0 in Ω × (0 , T end ) ,∂ρ υ ∂t + ρ υ ∇ · a + ∇ · ( ρ υ ⊗ c ) − ∇ · σ = ρ b in Ω × (0 , T end ) ,p ( ρ ) = 0 in Ω × (0 , T end ) ,ρ = ρ in Ω × (0) ,ρ υ = ρ υ in Ω × (0) ,ρ = ρ D on Γ D × (0 , T end ) ,ρ υ = ρ υ D on Γ D × (0 , T end ) , σ n = t N on Γ N × (0 , T end ) . (6)The pair ( ρ , ρ υ ) defines the initial conditions for the density and the momentum fields,while the quantities (cid:0) ρ D , ρ υ D (cid:1) and t N denote the Dirichlet and the Neumann boundarydata applied on Γ D and Γ N , respectively. Finally, n denotes the outward-pointing unitnormal vector to the corresponding boundary.Four types of boundary conditions are considered, namely inflow, outflow, no-slip andfree-slip conditions. At the inflow, the momentum profile is imposed via a Dirichlet bound-ary condition. At the outflow, the pressure (and therefore the density according to (4))is set to the given data, while the other quantities are extrapolated. For the no-slip con-dition, each velocity component (and therefore each momentum component) is forced tobe zero via a Dirichlet boundary condition, while the density is extrapolated. For thefree-slip condition instead, only the normal component of the velocity is set to zero, whilethe tangential component remains unconstrained. In this section, a hybridizable discontinuous Galerkin method is proposed for the solutionof weakly compressible flow problems satisfying equations (6). The so-called broken com-putational domain is defined by partitioning the fluid domain Ω in n el disjoint subdomainsΩ e : Ω = n el (cid:91) e =1 Ω e , Ω e ∩ Ω f = ∅ for e (cid:54) = f. (7)The internal element boundaries define the internal interfaceΓ := (cid:34) n el (cid:91) e =1 ∂ Ω e (cid:35) \ ∂ Ω (8)6nd the union of the internal interface with the boundary faces belonging to Γ N constitutesthe mesh skeleton, on which the hybrid variables are defined.The L scalar products for vector-valued functions in the elements and in their bound-aries are denoted in the following as:( w , υ ) Ω e := (cid:90) Ω e w · υ d Ω , (cid:104) w , υ (cid:105) ∂ Ω e := (cid:88) Γ i ⊂ ∂ Ω e (cid:90) Γ i w · υ d Γ . (9)Given the discontinuous character of the HDG variables, the jump operator (cid:74) · (cid:75) sums valuesfrom two adjacent elements Ω e and Ω f [36]: (cid:74) (cid:12) (cid:75) = (cid:12) e + (cid:12) f . (10)A key ingredient to preserve the convergence properties of the HDG method and to allowthe use of the same polynomial degree for the approximation of the primal and the mixedvariables without loss of accuracy is the adoption of the well-known Voigt notation [17, 45].The Voigt notation allows to strongly enforce the symmetry of the stress tensor in (2)by rearranging its diagonal and off-diagonal terms (according to the ordering in [13]) andstoring only the m sd independent components, with m sd = n sd ( n sd + 1) /
2. The dynamicviscosity and the second coefficient of viscosity can be embedded in the matrix D := (cid:20) µ I n sd + λ J n sd n sd × q sd q sd × n sd µ I q sd (cid:21) , (11)with q sd = m sd − n sd and J n sd denoting the n sd × n sd matrix of all ones. Hence, introducingthe operator ∇ S := (cid:20) ∂/∂x ∂/∂y ∂/∂y ∂/∂x (cid:21) T in 2D , ∂/∂x ∂/∂y ∂/∂z ∂/∂y ∂/∂x ∂/∂z ∂/∂z ∂/∂x ∂/∂y T in 3D , (12)and the vector E := (cid:20) n sd × q sd × (cid:21) , (13)the stress tensor can be expressed in Voigt notation as: σ V = − E p ( ρ ) + D ∇ S υ . (14)The normal component of the stress can be computed pre-multiplying σ V by n T , with n := (cid:20) n x n y n y n x (cid:21) T in 2D , n x n y n z n y n x n z n z n x n y T in 3D . (15)7n HDG methods, the solution of the overall problem is split into two phases [5,6,37–40].In the first stage, a local element-by-element Dirichlet problem is introduced to compute( L , ρ, ρ υ ) as a function of the unknown hybrid variables ( ˆ ρ, (cid:99) ρ υ ). With the notation justintroduced, the local problems can be written as: L + D ∇ S υ = in Ω e × (0 , T end ) ,∂ρ∂t + ρ ∇ · a + ∇ · ( ρ c ) = 0 in Ω e × (0 , T end ) ,∂ρ υ ∂t + ρ υ ∇ · a + ∇ · ( ρ υ ⊗ c )+ ∇ T S (cid:0) D L + E p ( ρ ) (cid:1) = ρ b in Ω e × (0 , T end ) ,ρ = ρ in Ω e × (0) ,ρ υ = ρ υ in Ω e × (0) ,ρ = ρ D on ∂ Ω e ∩ Γ D × (0 , T end ) ,ρ υ = ρ υ D on ∂ Ω e ∩ Γ D × (0 , T end ) ,ρ = ˆ ρ on ∂ Ω e \ Γ D × (0 , T end ) ,ρ υ = (cid:99) ρ υ on ∂ Ω e \ Γ D × (0 , T end ) , (16)for e = 1 , . . . , n el . The variable L denotes the aforementioned mixed variable, that allowsto reduce the second-order problem (6) to a system of first-order equations. In the secondstage, the traces of the density ˆ ρ and the momentum (cid:99) ρ υ are computed through the solutionof the global problem: (cid:74) ρ n (cid:75) = on Γ × (0 , T end ) , (cid:74) ρ υ ⊗ n (cid:75) = on Γ × (0 , T end ) , (cid:74) ρ c ⊗ n (cid:75) = on Γ × (0 , T end ) , (cid:74) ( ρ υ ⊗ c ) n (cid:75) = on Γ × (0 , T end ) , (cid:113) n T (cid:0) D L + E p ( ρ ) (cid:1) (cid:121) = on Γ × (0 , T end ) , − n T (cid:0) D L + E p ( ρ ) (cid:1) = t N on Γ N × (0 , T end ) . (17)These transmission conditions enforce the continuity of the primal variables ρ and ρ υ andthe normal fluxes across the interface Γ. The first two equations in (17) are automaticallysatisfied due to the unique definition of the hybrid variables on each face of the meshskeleton. Moreover, it is worth noting that, given the continuous nature of the velocity a of the moving background mesh, as opposed to the other HDG variables, no additionalconditions have to be enforced in the global problem.The following discrete functional spaces are introduced to derive the weak form of theproblem: W h (Ω) : = { w ∈ L (Ω) : w | Ω e ∈ P k (Ω e ) ∀ Ω e ⊂ Ω } , (18a) (cid:99) W h ( S ) : = { ˆ w ∈ L ( S ) : ˆ w | Γ i ∈ P k (cid:0) Γ i (cid:1) ∀ Γ i ⊂ S ⊆ Γ ∪ ∂ Ω } , (18b)8here P k (Ω e ) and P k (Γ i ) denote the spaces of polynomials of complete degree at most k in Ω e and on Γ i , respectively. The trace of the numerical normal fluxes, arising fromthe integration by parts of the terms under the divergence operator in (16), are defined asfollows: ρ c · n (cid:86) := (cid:40)(cid:0) ρ υ D − ρ D a (cid:1) · n + τ ρ (cid:0) ρ − ρ D (cid:1) on ∂ Ω e ∩ Γ D , ( (cid:99) ρ υ − ˆ ρ a ) · n + τ ρ ( ρ − ˆ ρ ) on ∂ Ω e \ Γ D , (19a)( ρ υ ⊗ c ) n (cid:86) := (cid:104) ρ υ D ⊗ (cid:0) ( ρ υ D /ρ D ) − a (cid:1)(cid:105) n + τ cρυ (cid:0) ρ υ − ρ υ D (cid:1) on ∂ Ω e ∩ Γ D , (cid:104) (cid:99) ρ υ ⊗ (cid:0) ( (cid:99) ρ υ / ˆ ρ ) − a (cid:1)(cid:105) n + τ cρυ ( ρ υ − (cid:99) ρ υ ) on ∂ Ω e \ Γ D , (19b) n T (cid:0) D L + E p ( ρ ) (cid:1) (cid:86) := (cid:40) n T (cid:0) D L + E p (cid:0) ρ D (cid:1)(cid:1) + τ dρυ (cid:0) ρ υ − ρ υ D (cid:1) on ∂ Ω e ∩ Γ D , n T (cid:0) D L + E p ( ˆ ρ ) (cid:1) + τ dρυ ( ρ υ − (cid:99) ρ υ ) on ∂ Ω e \ Γ D . (19c)These definitions are rather standard in the context of HDG methods and are closely relatedto other formulations published in literature. More precisely, the definition of the tracein the continuity equation is analogous to the one adopted for the solution of convection-diffusion problems in [37], whereas the definition of the traces in the momentum equationfollows the one proposed for the solution of flow problems in [17, 18]. However, the specificform of the definitions adopted in the context of the density-momentum formulation ofthe governing equations on deforming domains is original. The stabilization parameters τ ρ , τ cρυ and τ dρυ account for the compressibility, the convection and the diffusion effects,respectively, and they play a crucial role on the stability and the convergence of the HDGmethod [5, 6, 49]. Dimensional analysis provides a practical choice for the stabilizationparameters: τ ρ = C ρ ε | υ | , τ cρυ = C cρυ | υ | , τ dρυ = C dρυ µρ l , (20)with | υ | and l being a representative flow velocity and length scale, respectively, and C ρ , C cρυ and C dρυ denoting suitable positive scaling factors. It is empirically observedthat choosing the scaling factors in the range (1 ,
10) provides a good balance for thequality of the approximation of the primal, the mixed and the postprocessed variables,regardless of the polynomial degree, the type of element and the dimensionality of theproblem. These considerations are in agreement with the established results in the HDGliterature [5, 17, 28, 45]. Without loss of generality, a unique parameter τ ρυ = τ cρυ + τ dρυ ,taking into account both the convection and the diffusion effects, will be considered in thefollowing.With this definition of the numerical fluxes and expliciting all the unknowns, the dis-crete weak form of the local problems (16) reads: given ( ρ , ρ υ ) in Ω e × (0), (cid:0) ρ D , ρ υ D (cid:1) onΓ D and (cid:16) ˆ ρ h , (cid:99) ρ υ h (cid:17) on Γ ∪ Γ N , find (cid:0) L h , ρ h , ρ υ h (cid:1) ∈ (cid:2) W h (Ω e ) (cid:3) m sd × W h (Ω e ) × (cid:2) W h (Ω e ) (cid:3) n sd e = 1 , . . . , n el such that − (cid:0) L , L h (cid:1) Ω e + (cid:16) ∇ T S D L , ρ υ h ρ h (cid:17) Ω e = (cid:68) n T D L , ρ υ D ρ D (cid:69) ∂ Ω e ∩ Γ D + (cid:68) n T D L , (cid:99) ρ υ h ˆ ρ h (cid:69) ∂ Ω e \ Γ D , (21a) (cid:16) w, ∂ρ h ∂t (cid:17) Ω e + (cid:0) w, ρ h ∇ · a (cid:1) Ω e − (cid:0) ∇ w, ρ υ h − ρ h a (cid:1) Ω e + (cid:10) w, τ ρ ρ h (cid:11) ∂ Ω e = − (cid:10) w, (cid:0) ρ υ D − ρ D a (cid:1) · n − τ ρ ρ D (cid:11) ∂ Ω e ∩ Γ D − (cid:68) w, (cid:16) (cid:99) ρ υ h − ˆ ρ h a (cid:17) · n − τ ρ ˆ ρ h (cid:69) ∂ Ω e \ Γ D , (21b) (cid:16) w , ∂ρ υ h ∂t (cid:17) Ω e + (cid:0) w , ρ υ h ∇ · a (cid:1) Ω e − (cid:18) ∇ w , ρ υ h ⊗ (cid:18) ρ υ h ρ h − a (cid:19)(cid:19) Ω e + (cid:0) w , ∇ T S (cid:2) D L h + E p ( ρ h ) (cid:3)(cid:1) Ω e + (cid:10) w , τ ρυ ρ υ h (cid:11) ∂ Ω e − (cid:0) w , ρ h b (cid:1) Ω e = (cid:68) w , (cid:2) ρ υ D ⊗ (cid:0) ρ υ D ρ D − a (cid:1)(cid:3) n − τ ρυ ρ υ D (cid:69) ∂ Ω e ∩ Γ D (cid:68) w , (cid:2) (cid:99) ρ υ h ⊗ (cid:0) (cid:99) ρ υ h ˆ ρ h − a (cid:1)(cid:3) n − τ ρυ (cid:99) ρ υ h (cid:69) ∂ Ω e \ Γ D , (21c)for all ( L , w, w ) ∈ [ W h (Ω e )] m sd × W h (Ω e ) × [ W h (Ω e )] n sd .The discrete weak form of the global problem (17) instead reads: find ( ˆ ρ h , (cid:99) ρ υ h ) ∈ (cid:99) W h (Γ ∪ Γ N ) × [ (cid:99) W h (Γ ∪ Γ N )] n sd such that n el (cid:88) e =1 (cid:10) ˆ w, τ ρ ( ρ h − ˆ ρ h ) (cid:11) ∂ Ω e \ Γ D = 0 , (22a) − n el (cid:88) e =1 (cid:10) ˆ w , n T (cid:0) D L h + E p ( ˆ ρ h ) (cid:1) + τ ρυ ( ρ υ h − (cid:99) ρ υ h ) (cid:11) ∂ Ω e \ Γ D = n el (cid:88) e =1 (cid:10) ˆ w , t N (cid:11) ∂ Ω e ∩ Γ N , (22b)for all ( ˆ w, ˆ w ) ∈ (cid:99) W h (Γ ∪ Γ N ) × [ (cid:99) W h (Γ ∪ Γ N )] n sd . A key feature of the HDG method is the possibility to exploit the optimal convergenceof the mixed variable in order to construct a better approximation of the solution, con-verging in a superoptimal fashion. In HDG methods, the postprocessed variable is usually10he same physical quantity as the primal variable given by the solution of the local prob-lems, for instance the temperature in thermal problems, the displacement in elasticityproblems [7, 33, 45, 49] and the velocity in incompressible flow problems [17, 39, 40]. Thelocal postprocessing proposed in this contribution allows to construct a superconvergentvelocity field υ (cid:63) , although the primal variables are represented by the density ρ and themomentum ρ υ . To the best of the authors’ knowledge, such feature is not present in anyHDG formulation presented so far. The additional functional spaces are introduced: W (cid:63)h (Ω) : = { w (cid:63) ∈ L (Ω) : w (cid:63) | Ω e ∈ P k +1 (Ω e ) ∀ Ω e ⊂ Ω } , (23a) U h (Ω) : = { u ∈ L (Ω) : u | Ω e ∈ P (Ω e ) ∀ Ω e ⊂ Ω } . (23b)The discrete weak form of the local postprocessing then reads: given ( L h , ρ h , ρ υ h ) inΩ e , ( ρ D , ρ υ D ) on Γ D and ( ˆ ρ h , (cid:99) ρ υ h ) on Γ ∪ Γ N , find υ (cid:63)h ∈ [ W (cid:63)h (Ω e )] n sd for e = 1 , . . . , n el such that − (cid:16) ∇ S w (cid:63) , D ∇ S υ (cid:63)h (cid:17) Ω e = (cid:0) ∇ S w (cid:63) , L h (cid:1) Ω e , (24a) (cid:0) u T , υ (cid:63)h (cid:1) Ω e = (cid:16) u T , ρ υ h ρ h (cid:17) Ω e , (24b) (cid:0) u R , ∇ W υ (cid:63)h (cid:1) Ω e = (cid:68) u R , T ρ υ D ρ D (cid:69) ∂ Ω e ∩ Γ D + (cid:68) u R , T (cid:99) ρ υ h ˆ ρ h (cid:69) ∂ Ω e \ Γ D , (24c)for all ( w (cid:63) , u T , u R ) ∈ [ W (cid:63)h (Ω e )] n sd × [ U h (Ω e )] n sd × [ U h (Ω e )] q sd . The vorticity operator in (24)is defined in Voigt notations as ∇ W := (cid:104) − ∂/∂y ∂/∂x (cid:105) in 2D , − ∂/∂z ∂/∂y∂/∂z − ∂/∂x − ∂/∂y ∂/∂x in 3D , (25)while the matrix T accounts for the tangential direction to the boundary and it is definedas T := (cid:104) − n y n x (cid:105) in 2D , − n z n y n z − n x − n y n x in 3D . (26)The postprocessing presented here is formally similar to the one proposed in [17], but witha different definition of the matrix D (including here the second coefficient of viscosity λ )and with the presence of the ratio between the momentum and the density instead of thevelocity itself. The first equation in (24) directly follows from the definition of the mixed11ariable in (16) and is a least-squares fit to the accurate variable L h , while the last twoequations remove the underdetermination of the problem by constraining the rigid motions.It is worth recalling that the space U h , containing the functions of all ones in the elements,in the vector case selectively goes through the various components, hence providing n sd constraints for the translations and q sd constraints for the rotations in (24). In this section, the governing equations of nonlinear elastodynamics are first presentedtogether with the associated standard CG formulation. Then, the conditions to couple thefluid and the structural fields are briefly presented and two different coupling strategies,namely the Dirichlet–Neumann coupling and the Nitsche-based coupling, are proposedand their distinctive properties discussed. No specific indices have been used in section2 to refer to the fluid quantities in order ease the comprehension of the proposed HDGformulation for weakly compressible flows. In the following, however, the fluid and thestructural quantities are distinguished by means of the subscripts ( · ) F and ( · ) S . The strong form of the time-dependent nonlinear elastic problem can be written withrespect to the undeformed structural domain Ω S as: ρ S d u S dt − ∇ · P S = ρ S b S in Ω S × (0 , T end ) , u S = u S in Ω S × (0) ,d u S dt = ˙ u S in Ω S × (0) , u S = u D S on Γ D S × (0 , T end ) , P S n S = t N S on Γ N S × (0 , T end ) , (27)where u S represents the unknown displacement field, ρ S the structural density, P S the firstPiola–Kirchhoff stress tensor and b S an external body force per unit undeformed volume.The pair ( u S , ˙ u S ) defines the initial conditions for the displacement and the velocity, whilethe quantities u D S and t N S denote the Dirichlet and the Neumann boundary data applied onΓ D S and Γ N S , respectively. For hyperelastic materials, the first Piola–Kirchhoff stress tensoris defined as P S = ∂ψ S ∂ F S , (28)with F S being the deformation gradient and ψ S the strain energy density function. Theformer is derived from the displacement field as F S = ∇ u S + I n sd , (29)12hile the latter is defined as ψ S = (cid:40) µ S (cid:2) tr( F T S F S ) − n sd − | F S | ) (cid:3) + λ S [ln ( | F S | )] for Neo-Hooke ,µ S E S : E S + λ S [tr ( E S )] with E S = ( F T S F S − I n sd ) for St. Venant–Kirchhoff , (30)for the two popular material models used here. The Lam´e parameters µ S and λ S can beevaluated as functions of the Young modulus E S and the Poisson ratio ν S of the materialthrough the following relations: µ S = E S ν S ) and λ S = ν S E S (1 + ν S ) (1 − ν S ) . (31)The following discrete functional spaces are introduced: V h (Ω) : = { v ∈ H (Ω) : v | Ω e ∈ P k (Ω e ) ∀ Ω e ⊂ Ω , v | Γ D = u D S } , (32a) V h (Ω) : = { v ∈ H (Ω) : v | Ω e ∈ P k (Ω e ) ∀ Ω e ⊂ Ω , v | Γ D = 0 } . (32b)As usual in CG methods, the weak form of the problem is obtained by multiplying thegoverning equation with the test functions and integrating by parts the term with secondorder derivatives. The discrete weak form of the structural problem then reads: given( u S , ˙ u S ) in Ω S × (0), find u h S ∈ [ V h (Ω S )] n sd such that (cid:16) v , ρ S d u h S dt (cid:17) Ω S + (cid:0) ∇ v , P h S (cid:1) Ω S = ( v , ρ S b S ) Ω S + (cid:10) v , t N S (cid:11) Γ N S , (33)for all v ∈ [ V h (Ω S )] n sd . Of course, P h S features a nonlinear dependence on the displacementthrough (30) and a Newton–Raphson procedure is utilized for the solution of (33). In order to couple the fluid problem, whose weak form has been derived section in 2, andthe structural problem, whose weak form has been presented in section 3.1, kinematicand dynamic continuity conditions have to be enforced at the fluid-structure interfaceΓ I = Ω F ∩ Ω S . First, the no-slip condition υ F − υ S = (34)prohibits a fluid flow across the interface and a relative tangential movement of fluid andstructure at the interface. Here, the fluid velocity is evaluated as the ratio of the momentumand the density, while the structural velocity is computed as the time derivative of thestructural unknown displacement. Second, the traction equilibrium t F + t S = (35)13igure 1: Degrees of freedom of the coupled HDG-CG discretization using polynomialapproximation of degree k = 2 in the HDG fluid subdomain Ω F (in blue) and in the CGstructural subdomain Ω S (in red).states the equilibrium of the fluid and the structural forces at the interface. The way inwhich the coupling conditions (34) and (35) are imposed differs depending on the couplingstrategy adopted.It is worth mentioning that the deformation of the fluid computational mesh is evaluatedas a function of the structural displacement at the interface, by means of a unique ALEmapping d = ϕ ( u S ) . (36)The grid motion strategy is an artificial problem and does not affect the physics of thecoupled problem. Its sole purpose is to generate a proper mesh for the solution of the fluidproblem. The velocity of the computational mesh, which is independent of the velocity ofthe material particles, is computed as the time derivative of the mesh displacement.In this contribution, special attention is devoted to the spatial discretization of theweakly compressible flow problem by means of the HDG method and to the couplingof the fluid field with the structural one, discretized by means of the CG method. Asketch of the heterogeneous HDG-CG discretization is exemplarily shown in Figure 1. Thetemporal discretization does not represent the focus of this work and the implicit backwarddifferentiation formulas (BDF) are adopted here for the sake of simplicity. Better timeintegration schemes have been developed in literature, such as the generalized- α methodfor the fluid [26] and the structure [4], and the possibility of independently choosing themin order to meet the needs of the individual fields has been exploited in [34].For the sake of readability, the superscript ( · ) h associated with the numerical approx-imation of the unknowns will be henceforth omitted. Moreover, it is assumed that thefluid quantities refer to the deformed domain Ω F and the structural quantities refer to theundeformed domain Ω S . 14 .3 A partitioned Dirichlet–Neumann algorithm for the HDG-CG coupling The first coupling strategy is a partitioned Dirichlet–Neumann coupling, since it builds aDirichlet-to-Neumann map by taking the structural solution at the interface as a Dirichletboundary condition for the fluid problem and imposing the fluid normal flux as a Neumannboundary condition for the structural problem. This represents a popular partitionedscheme for the solution of FSI problems and it has been used and analyzed for instancein [8,31]. In [32] it has moreover been shown that the introduction of a weak compressibilityin the fluid field alleviates the constraints of the instability condition of the artificial addedmass effect and, in comparison to a fully incompressible solver, it reduces the number ofcoupling iterations required and it leads to an increase of the dynamic relaxation parameter.In this contribution, the partitioned Dirichlet–Neumann scheme is revisited for thecoupling of weakly compressible flows and elastic structures in which the fluid and thestructural problems are discretized by means of the HDG and the CG method, respectively.The solver coupling at each time steps can be schematized in the following steps:1. Predict the structural displacement u S on Γ I , assuming for instance a constant dis-placement, velocity or acceleration field.2. Update the fluid mesh configuration by means of the ALE mapping (36).3. Solve the fluid problem on the newly deformed domain by imposing the velocitycompatibility (34) as a Dirichlet-type boundary condition on the local problems as ρ υ F ρ F − d u S dt = on Γ I . (37)The weak form of the fluid local problems then reads: given ( ρ F , ρ υ F ) in Ω e F × (0),( ρ D F , ρ υ D F ) on Γ D F , ( ˆ ρ F , (cid:99) ρ υ F ) on Γ F ∪ Γ N F and ( ˆ ρ F , υ S ) on Γ I , find ( L F , ρ F , ρ υ F ) ∈ [ W (Ω e F )] m sd × W (Ω e F ) × [ W (Ω e F )] n sd for e = 1 , . . . , n el F such that − ( L , L F ) Ω e F + (cid:16) ∇ T S D F L , ρ υ F ρ F (cid:17) Ω e F = (cid:68) n T F D F L , ρ υ D F ρ D F (cid:69) ∂ Ω e F ∩ Γ D F + (cid:68) n T F D F L , (cid:99) ρ υ F ˆ ρ F (cid:69) ∂ Ω e F \ Γ D F \ Γ I + (cid:68) n T F D F L , d u S dt (cid:69) ∂ Ω e F ∩ Γ I , (38a) (cid:16) w, ∂ρ F ∂t (cid:17) Ω e F +( w, ρ F ∇ · a F ) Ω e F − ( ∇ w, ρ υ F − ρ F a F ) Ω e F + (cid:10) w, τ ρ ρ F (cid:11) ∂ Ω e F = − (cid:10) w, ( ρ υ D F − ρ D F a F ) · n F − τ ρ ρ D F (cid:11) ∂ Ω e F ∩ Γ D F − (cid:10) w, ( (cid:99) ρ υ F − ˆ ρ F a F ) · n F − τ ρ ˆ ρ F (cid:11) ∂ Ω e F \ Γ D F \ Γ I − (cid:10) w, ˆ ρ F (cid:0) d u S dt − a F (cid:1) · n F − τ ρ ˆ ρ F (cid:11) ∂ Ω e F ∩ Γ I , (38b)15 w , ∂ρ υ F ∂t (cid:17) Ω e F +( w , ρ υ F ∇ · a F ) Ω e F − (cid:0) ∇ w , ρ υ F ⊗ (cid:0) ρ υ F ρ F − a F (cid:1)(cid:1) Ω e F + (cid:16) w , ∇ T S (cid:0) D F L F + E p F ( ρ F ) (cid:1)(cid:17) Ω e F + (cid:10) w , τ ρυ ρ υ F (cid:11) ∂ Ω e F − ( w , ρ F b F ) Ω e F = − (cid:68) w , (cid:2) ρ υ D F ⊗ (cid:0) ρ υ D F ρ D F − a F (cid:1)(cid:3) n F − τ ρυ ρ υ D F (cid:69) ∂ Ω e F ∩ Γ D F − (cid:68) w , (cid:2) (cid:99) ρ υ F ⊗ (cid:0) (cid:99) ρ υ F ˆ ρ F − a F (cid:1)(cid:3) n F − τ ρυ (cid:99) ρ υ F (cid:69) ∂ Ω e F \ Γ D F \ Γ I − (cid:68) w , (cid:2) ˆ ρ F d u S dt ⊗ (cid:0) d u S dt − a F (cid:1)(cid:3) n F − τ ρυ ˆ ρ F d u S dt (cid:69) ∂ Ω e F ∩ Γ I , (38c)for all ( L , w, w ) ∈ [ W (Ω e F )] m sd × W (Ω e F ) × [ W (Ω e F )] n sd . From a practical point ofview, equations (38) are obtained from (21) by replacing (cid:99) ρ υ F with ˆ ρ F d u S dt at thefluid-structure interface Γ I .The weak form of the fluid global problem instead reads: find ( ˆ ρ F , (cid:99) ρ υ F ) ∈ (cid:99) W (Γ F ∪ Γ N F ∪ Γ I ) × [ (cid:99) W (Γ F ∪ Γ N F )] n sd such that n el F (cid:88) e =1 (cid:10) ˆ w, τ ρ ( ρ F − ˆ ρ F ) (cid:11) ∂ Ω e F \ Γ D F = 0 , (39a) − n el F (cid:88) e =1 (cid:10) ˆ w , n T F (cid:0) D F L F + E p F ( ˆ ρ F ) (cid:1) + τ ρυ ( ρ υ F − (cid:99) ρ υ F ) (cid:11) ∂ Ω e F \ Γ D F \ Γ I = n el F (cid:88) e =1 (cid:104) ˆ w , t N F (cid:105) ∂ Ω e F ∩ Γ N F , (39b)for all ( ˆ w, ˆ w ) ∈ (cid:99) W (Γ F ∪ Γ N F ∪ Γ I ) × [ (cid:99) W (Γ F ∪ Γ N F )] n sd .4. Solve the structural problem by imposing the coupling condition (35) as a Neumann-type boundary condition, consistently expressed in terms of the first Piola–Kirchhoffstress tensor as ( P F − P S ) n S = on Γ I , (40)where the fluid Cauchy stress is transformed by means of the pull-back operation P F = −| F F | V − ( D F L F + E p F ) F − T F . (41)The deformation gradient F F is evaluated as in (29) but with respect to the fluid meshdisplacement. Given a m sd × − returns theassociated n sd × n sd symmetric tensor:V − := (cid:104) σ xx σ yy σ xy (cid:105) T −→ (cid:34) σ xx σ xy σ xy σ yy (cid:35) in 2D , (cid:104) σ xx σ yy σ zz σ xy σ xz σ yz (cid:105) T −→ σ xx σ xy σ xz σ xy σ yy σ yz σ xz σ yz σ zz in 3D . (42)16nalogously to the viscous stress in (19), the following trace of the numerical normalflux is introduced: P F n F (cid:86) := −| F F | V − ( D F L F + E p F ( ˆ ρ F )) F − T F n F − τ ρυ (cid:0) ρ υ F − ˆ ρ F d u S dt (cid:1) on Γ I . (43)The weak form of the structural problem then reads: given ( u S , ˙ u S ) in Ω S × (0) and t F on Γ I , find u S ∈ [ V (Ω S )] n sd such that (cid:16) v , ρ S d u S dt (cid:17) Ω S + (cid:0) ∇ v , P S (cid:1) Ω S + (cid:10) v , | F F | V − (cid:0) D F L F + E p F ( ˆ ρ F ) (cid:1) F − T F n S − τ ρυ (cid:0) ρ υ F − ˆ ρ F d u S dt (cid:1)(cid:11) Γ I = (cid:0) v , ρ S b S (cid:1) Ω S + (cid:10) v , t N S (cid:11) Γ N S , (44)for all v ∈ [ V (Ω S )] n sd .5. Check for convergence: continue with next time step if the algorithm is converged,otherwise return to step 2. The convergence is considered satisfied if r i +1 S = (cid:107) ˜ u i +1 S − u i S (cid:107) < η on Γ I , (45)with u i S denoting the structural interface displacement at the i -th coupling iterationand ˜ u i +1 S the newly computed one by solving (44). The parameter η represents insteada user-defined convergence tolerance.To accelerate the convergence of the fixed-point scheme, a relaxation of the structuralinterface displacement is performed u i +1 S = ω i ˜ u i +1 S + (1 − ω i ) u i S , (46)where the relaxation parameter ω is evaluated at each coupling iteration by meansof the Aitken ∆ method [31]: ω i +1 = − ω i ( r i +1 S ) T ( r i +2 S − r i +1 S ) (cid:107) r i +2 S − r i +1 S (cid:107) . (47) Remark 1.
The strategy presented in this section to couple the fluid and the structurecould also be implemented in a monolithic fashion, by simultaneously solving the problems (38) - (39) - (44) in a large linear system. An analogous method to couple HDG and CGdiscretization has been presented in [41] in the context of conjugate heat transfer problems.However, such a method results in a coupling of local and global degrees of freedom of theHDG problem with the ones of the CG discretization, making the implementation of thisstrategy in existing HDG and CG libraries rather intrusive. .4 A monolithic algorithm for the HDG-CG coupling based onNitsche’s method The second coupling strategy imposes the structural numerical normal flux as a Neumannboundary condition for the fluid problem and takes the fluid hybrid variables at the in-terface as a Dirichlet-type boundary condition for the structural problem. The resultinghybrid HDG-CG coupling does not affect the structure of the core CG and HDG matrices,thus leading to a minimally-intrusive implementation of this technique in existing finiteelement codes. This approach to couple HDG and CG discretization has been recentlypresented in [33] for the solution of elasticity problems involving compressible and nearlyincompressible solids.The key feature of the Nitsche-based coupling is that it allows to impose the couplingconditions solely in the global problem and, as a consequence, the local problems remainthe same as in the pure HDG case (equations (21)).With the notation introduced in section 3.3 for the coupled problem, the HDG localproblems read: given ( ρ F , ρ υ F ) in Ω e F × (0), ( ρ D F , ρ υ D F ) on Γ D F and ( ˆ ρ F , (cid:99) ρ υ F ) on Γ F ∪ Γ N F ∪ Γ I ,find ( L F , ρ F , ρ υ F ) ∈ [ W (Ω e F )] m sd × W (Ω e F ) × [ W (Ω e F )] n sd for e = 1 , . . . , n el F such that − (cid:0) L , L F (cid:1) Ω e F + (cid:0) ∇ T S D F L , ρ υ F ρ F (cid:1) Ω e F = (cid:10) n T F D F L , ρ υ D F ρ D F (cid:11) ∂ Ω e F ∩ Γ D F + (cid:10) n T F D F L , (cid:99) ρ υ F ˆ ρ F (cid:11) ∂ Ω e F \ Γ D F , (48a) (cid:16) w, ∂ρ F ∂t (cid:17) Ω e F + (cid:0) w, ρ F ∇ · a F (cid:1) Ω e F − (cid:0) ∇ w, ρ υ F − ρ F a F (cid:1) Ω e F + (cid:10) w, τ ρ ρ F (cid:11) ∂ Ω e F = − (cid:10) w, ( ρ υ D F − ρ D F a F ) · n F − τ ρ ρ D F (cid:11) ∂ Ω e F ∩ Γ D F − (cid:10) w, ( (cid:99) ρ υ F − ˆ ρ F a F ) · n F − τ ρ ˆ ρ F (cid:11) ∂ Ω e F \ Γ D F , (48b) (cid:16) w , ∂ρ υ F ∂t (cid:17) Ω e F + (cid:0) w , ρ υ F ∇ · a F (cid:1) Ω e F − (cid:0) ∇ w , ρ υ F ⊗ (cid:0) ρ υ F ρ F − a F (cid:1)(cid:1) Ω e F + (cid:0) w , ∇ T S (cid:0) D F L F + E p F ( ρ F ) (cid:1)(cid:1) Ω e F + (cid:10) w , τ ρυ ρ υ F (cid:11) ∂ Ω e F − (cid:0) w , ρ F b F (cid:1) Ω e F = − (cid:10) w , (cid:2) ρ υ D F ⊗ (cid:0) ρ υ D F ρ D F − a F (cid:1)(cid:3) n F − τ ρυ ρ υ D F (cid:11) ∂ Ω e F ∩ Γ D F − (cid:10) w , (cid:2) (cid:99) ρ υ F ⊗ (cid:0) (cid:99) ρ υ F ˆ ρ F − a F (cid:1)(cid:3) n F − τ ρυ (cid:99) ρ υ F (cid:11) ∂ Ω e F \ Γ D F , (48c)for all ( L , w, w ) ∈ [ W (Ω e F )] m sd × W (Ω e F ) × [ W (Ω e F )] n sd . It is worth recalling that equations(48) are identical to (21) and they are rewritten here for completeness with the specificnotation adopted for fluid-structure interaction problems.18he coupling condition (35) is imposed as a Neumann-type boundary condition on theglobal fluid problem and it is consistently expressed in terms of the Cauchy stress tensor.The dynamic equilibrium can therefore be written as( σ F − σ S ) n F = on Γ I , (49)where the first Piola–Kirchhoff stress is transformed by means of the push-forward opera-tion σ S = | F S | − P S F S T . (50)In order to impose the coupling conditions on the interface, the following definition ofthe trace of the CG numerical normal flux is adopted: P S n S (cid:86) := P S n S − γh (cid:0) d u S dt − (cid:99) ρ υ F ˆ ρ F (cid:1) on Γ I , (51)where h denotes a characteristic element size on Γ I and γ is a sufficiently large positivevalue, called in the following “Nitsche’s parameter”, that is commonly used to enforcecoercivity of the discrete bilinear form in CG discretizations with Nitsche’s imposition ofessential boundary conditions. The influence of the Nitsche parameter on the accuracy ofthe hybrid HDG-CG coupling has been investigated in a previous work [33] on a computa-tionally cheap scalar problem. The analysis revealed that for low values of γ an insufficientstabilization produces unreliable results, whereas values of γ above a certain lower boundensure the stability of the scheme, in agreement with established results in literature [22].However, although an estimation of such a lower bound can be obtained by solving anauxiliary generalized eigenvalue problem as suggested in [21], γ is problem-dependent andaffected both by the equation under analysis and the material parameters. In terms of thescaling of the penalty parameter with the polynomial degree, a k scaling can be expectedfrom the results established by the symmetric interior penalty community [9], a topic notfurther investigated in the present work.It is worth recalling that in the standard CG approach the numerical normal fluxes arenaturally equilibrated on the internal faces of the triangulation. The velocity compatibility(34) is weakly imposed as a Dirichlet-type boundary condition on the structural problem,exploiting the definition (51).The weak form of the global problem then reads: given ( u S , ˙ u S ) in Ω S × (0), find( ˆ ρ F , (cid:99) ρ υ F , u S ) ∈ (cid:99) W (Γ F ∪ Γ N F ∪ Γ I ) × [ (cid:99) W (Γ F ∪ Γ N F ∪ Γ I )] n sd × [ V (Ω S )] n sd such that n el F (cid:88) e =1 (cid:10) ˆ w, τ ρ ( ρ F − ˆ ρ F ) (cid:11) ∂ Ω e F \ Γ D F = 0 , (52a)19 n el F (cid:88) e =1 (cid:26)(cid:10) ˆ w , n T F (cid:0) D F L F + E p F ( ˆ ρ F ) (cid:1) + τ ρυ ( ρ υ F − (cid:99) ρ υ F ) (cid:11) ∂ Ω e F \ Γ D F + (cid:68) ˆ w , | F S | − P S F S T n F + γh (cid:16) d u S dt − (cid:99) ρ υ F ˆ ρ F (cid:17)(cid:69) ∂ Ω e F ∩ Γ I (cid:27) = n el F (cid:88) e =1 (cid:10) ˆ w , t N F (cid:11) ∂ Ω e F ∩ Γ N F , (52b) (cid:16) v , ρ S d u S dt (cid:17) Ω S + (cid:0) ∇ v , P S (cid:1) Ω S − (cid:68) v , P S n S − γh (cid:16) d u S dt − (cid:99) ρ υ F ˆ ρ F (cid:17)(cid:69) Γ I − (cid:68) ∂ P S ∂ ∇ u S ∇ v n S , d u S dt − (cid:99) ρ υ F ˆ ρ F (cid:69) Γ I = (cid:0) v , ρ S b S (cid:1) Ω S + (cid:10) v , t N S (cid:11) Γ N S , (52c)for all ( ˆ w, ˆ w , v ) ∈ (cid:99) W (Γ F ∪ Γ N F ∪ Γ I ) × [ (cid:99) W (Γ F ∪ Γ N F ∪ Γ I )] n sd × [ V (Ω S )] n sd .After standard finite element discretization and assembly, the following linear systemis obtained in terms of the increments of the global unknowns (cid:98) U F = (cid:2) ˆ ρ F (cid:99) ρ υ F (cid:3) T and u S : (cid:20) K FF K FS K SF K SS (cid:21) (cid:20) δ (cid:98) U F δ u S (cid:21) = (cid:20) f F f S (cid:21) , (53)with the left hand side matrices computed as (cid:2) K FF (cid:3) = n el F (cid:88) e =1 (cid:110)(cid:2) K (cid:98) U (cid:98) U (cid:3) e − (cid:2) K (cid:98) UL K (cid:98) UU (cid:3) e (cid:20) K LL K LU K UL K UU (cid:21) − e (cid:20) K L (cid:98) U K U (cid:98) U (cid:21) e (cid:111) , (54a) (cid:2) K FS (cid:3) = n el I (cid:88) e =1 (cid:110)(cid:2) K (cid:98) Uu (cid:3) e (cid:111) , (54b) (cid:2) K SF (cid:3) = n el I (cid:88) e =1 (cid:110)(cid:2) K u (cid:98) U (cid:3) e (cid:111) , (54c) (cid:2) K SS (cid:3) = n el S (cid:88) e =1 (cid:110)(cid:2) K uu (cid:3) e (cid:111) , (54d)and the right hand side vectors computed as (cid:2) f F (cid:3) = n el F (cid:88) e =1 (cid:110)(cid:2) f (cid:98) U (cid:3) e − (cid:2) K (cid:98) UL K (cid:98) UU (cid:3) e (cid:20) K LL K LU K UL K UU (cid:21) − e (cid:20) f L f U (cid:21) e (cid:111) , (55a) (cid:2) f S (cid:3) = n el S (cid:88) e =1 (cid:110)(cid:2) f u (cid:3) e (cid:111) . (55b)20he terms n el F and n el S denote the number of elements in the fluid and structural discretiza-tion, respectively, whereas n el I refers to the number of elements adjacent to the interface(belonging either to the fluid or the structural subdomain). In the formulas (54)–(55), thesummation over elements is understood as the usual assembly process, adding the localmatrices and vectors into the associated positions of the global matrices and vectors. Itcan be observed from (54) how K FF and K SS feature the usual structure of the matricesof the HDG and CG global problem, respectively, and they differ from the single-field ma-trices only for the inclusion of a small number of terms in K (cid:98) U (cid:98) U and K uu arising from thedefinition (51). The blocks K FS and K SF are responsible for the coupling and they simplystem from the linearization and discretization of − (cid:10) ˆ w , | F S | − P S F S T n F + γh − ( d u S /dt ) (cid:11) and (cid:10) ( ∂ P S /∂ ∇ u S ) ∇ vn S , (cid:99) ρ υ F ˆ ρ − F (cid:11) − (cid:10) v , γh − (cid:99) ρ υ F ˆ ρ − F (cid:11) along the interface, respectively. The vec-tors f F and f S represent the residuals of the fluid and the structural global problems. In thespirit of the Nitsche-based coupling of HDG and CG discretizations proposed in [33], thefluid HDG local problems remain unchanged and they require, at each Newton iteration,the solution of the linear systems (cid:20) K LL K LU K UL K UU (cid:21) e (cid:20) δ L F δ U F (cid:21) e = (cid:20) f L f U (cid:21) e − (cid:20) K L (cid:98) U K U (cid:98) U (cid:21) e (cid:104) δ (cid:98) U F (cid:105) e , (56)for e = 1 , . . . , n el F .Since here the coupling takes place only at a global level, unlike the monolithic versionof the Dirichlet–Neumann coupling mentioned in section 3.3, the communication of theHDG local matrices is not required outside the fluid block of the matrix and the righthand side in (53). This segregation between the HDG local DOFs and and the CG DOFsleads to the minimally-intrusive computer implementation already mentioned. Moreover,the HDG local problems and the HDG local postprocessing remain the same as in the pureHDG case and no special treatment of the interface elements is required. Last but notleast, since the Nitsche-based coupling solely relies on the hybrid variables to impose thecoupling conditions, the treatment of non-matching grids and/or non-uniform polynomialdegrees is easily handled without the need of introducing special projection operators. Remark 2.
The strategy presented in this section to couple the fluid and the structurecould in theory be implemented in a partitioned fashion, by alternating the solution of purefluid and structural problems and exchanging the interface information among the fields.This method can be referred to as “partitioned Neumann–Dirichlet coupling” and it hasbeen proposed for instance in [30] as a possible remedy for the so-called incompressibilitydilemma. However, as stated by the same authors of [30], such a method fails to solve realworld problems, since the response of stiff structures to varying interface displacements willbe too sensitive for any numerical approach to find the equilibrium. Numerical studies
In this section several numerical studies are presented to assess the performance of theproposed HDG-CG formulation for weakly compressible fluid-structure interaction. Theconvergence properties of the HDG formulation for weakly compressible flows introducedin section 2 are analyzed first with respect to a simple steady state Poiseuille flow equippedwith analytical solution and then with respect to an unsteady flow on a moving mesh whichexercises all the terms present in the fluid PDEs. The third example verifies the optimalconvergence of the HDG-CG coupling schemes on a problem with manufactured solutionand the following examples solve weakly compressible fluid-structure interaction problemson two and three dimensions.
The first numerical example considers a steady state isothermal Poiseuille flow of a weaklycompressible Newtonian fluid in a straight channel. The goal of this study is to show ona simple and physically meaningful example the convergence properties of the proposedHDG formulation, as well as its robustness with respect to the compressibility level. Inthe work [24], the authors derive an analytical solution by representing the primary flowvariables as asymptotic expansions of the compressibility coefficient, which is assumed to bea small parameter, and perturbing them with respect to the same coefficient. The solutionis then found up to the first order in ε . The study [24] considers also a pressure-dependentviscosity, but this feature has been neglected because unimportant in the context of thepresent contribution. Since this example concerns a pure flow problem, the subscript ( · ) F is omitted for brevity.The analytical solution in terms of velocity and pressure can be written as υ x ( x, y ) = 32 U (cid:104) − (cid:16) yR (cid:17) (cid:105) − µLU ρ R (cid:16) − xL (cid:17)(cid:104) − (cid:16) yR (cid:17) (cid:105) ε,υ y ( x, y ) = 0 ,p ( x, y ) = p + 3 µLUR (cid:16) − xL (cid:17) − µ U ρ R (cid:110) (cid:16) LR (cid:17) (cid:16) − xL (cid:17) − (cid:104) − (cid:16) yR (cid:17) (cid:105)(cid:111) ε, (57)where L and R represent the length and the half-height of the channel, respectively, while U denotes the mean velocity at the channel exit. The solution of the density is derivedfrom the equation of state (4) using the expression of the pressure in (57) and the solutionof the momentum is then obtained by multiplying the density just derived with the velocityfield in (57). No body forces in the momentum equation appear in this expansion. Dueto the asymptotic expansion of the solution in terms of ε up to the first order, a residual O ( ε ) is added to the right hand side of the continuity equation in (1) in the spirit of22 a) υ x ( x,
0) (b) p ( x, Figure 2: Axial velocity (left) and pressure (right) along the horizontal axis of the weaklycompressible Poiseuille flow with different values of the dimensionless compressibility co-efficient.manufactured solutions: R ( x, y ) = 27 µ U ρ R (cid:0) R − y (cid:1) (cid:8) R ( L − x ) ε − µU ρ − (cid:2) L − x ) − (cid:0) R − y (cid:1)(cid:3) ε (cid:9) . (58)From the expressions (57) it can be easily observed that the solution of a classic incompress-ible Poiseuille flow in a straight channel is fully recovered when ε →
0. In this situation theaxial velocity assumes a simple parabolic profile, while the pressure varies linearly alongthe channel.The fluid domain is the rectangle Ω = [0 , L ] × [ − R, R ], with L = 10 and R = 1. Theviscosity µ is considered equal to 1 as well as the mean velocity at the channel exit U and the reference density ρ , evaluated at the reference pressure p = 0. A dimensionlesscompressibility number can be defined as ε ∗ = 3 µLUρ R ε (59)and three different orders of magnitude are considered in the following studies, i.e., ε ∗ =[0 . , . , . U = 1 at the channel exit). Clearly, thesolution has a physical meaning only when the volumetric flow rate at the entrance of thechannel is positive, therefore ε ∗ = 1 constitutes an upper limit of validity of the analyticalsolution in terms of the compressibility coefficient. Although this last case lacks physical23 a) Mesh with m = 1 (b) Mesh with m = 2 (c) Mesh with m = 3 Figure 3: First three levels of refinement of the mesh used for the convergence studies ofthe weakly compressible Poiseuille flow.meaningfulness, it is interesting to investigate the performance of the numerical methodproposed also in this extreme case. With regards to the pressure, it varies quadraticallyalong the horizontal axis and it always reaches the reference value ( p = 0) at the channelexit corners. Moreover, the average pressure drop required to drive the flow decreases withthe compressibility.In the following experiments, Dirichlet boundary conditions, corresponding to the re-striction of the analytical solution to the domain boundary, are imposed on Γ D = ∂ Ω.Uniform meshes of triangular elements are considered for the whole domain by splitting aregular 2 m × m Cartesian grid (with m ranging from 1 to 5) into a total of 2 m +1 triangles,giving element sizes of h = L/ m . The first three levels of refinement of the mesh usedfor the convergence studies are shown in Figure 3. The degree of approximation k usedranges from 1 to 2 for ε ∗ = 0 .
01, from 1 to 3 for ε ∗ = 0 . ε ∗ = 1.The stabilization parameters are computed according to (20), considering | υ | = U and l = R as representative velocity and length, respectively. The scaling factors are chosenas C ρ = 3 .
33 and C dρυ = 1 (no convective effects are included in this flow configuration),returning the stabilization parameters τ ρ = 3 . /ε and τ ρυ = 1 for the density and themomentum, respectively. In Figure 4 the solution of the density and the momentum fieldobtained with the proposed HDG formulation using m = 5 and k = 2 is shown. Withregards to the density, its maximum variation from the reference value is about 1% for ε ∗ = 0 .
01, 10% for ε ∗ = 0 . ε ∗ = 1, for which the maximum value of thedensity reaches 1 .
50. In Figure 5 the improvement of the approximation of the velocityfield given by the local postprocessing described in section 2.3 is exemplarily shown for theintermediate compressibility coefficient ( ε ∗ = 0 .
1) using m = 2 and k = 1.The convergence of the error measured in the L norm as a function of the characteristicelement size h is represented in Figure 6, for the different compressibility coefficients con-sidered. It is worth noting that the error decreases by orders of magnitude when decreasing24 a) ρ with ε ∗ = 0 .
01 (b) ρ with ε ∗ = 0 . ρ with ε ∗ = 1(d) ρυ x with ε ∗ = 0 .
01 (e) ρυ x with ε ∗ = 0 . ρυ x with ε ∗ = 1 Figure 4: Approximation of the density and the momentum field of the weakly compressiblePoiseuille flow with different values of the dimensionless compressibility coefficient. (a) υ x with ε ∗ = 0 . υ (cid:63)x with ε ∗ = 0 . Figure 5: Approximation of the velocity and the postprocessed velocity field of the weaklycompressible Poiseuille flow with ε ∗ = 0 .
1. 25he compressibility coefficient because in the limit case of ε = 0 the solution belongs to thespace of quadratic polynomials P (Ω). On the other hand, the analytical solution is fullyrecovered up to machine precision with k ≥ ε , since the pressure (and therefore the density according to (4)) belongs to P (Ω) andthe velocity belongs to P (Ω), thus the momentum belongs to P (Ω). The convergencestudies are therefore performed with variable ranges of polynomial degrees for the differentcompressibility coefficients, in order to clearly visualize the convergence rates in the asymp-totic regime by avoiding the errors to reach the machine precision for excessively coarsemeshes. Optimal convergence rates (with order k + 1) are obtained for the mixed variable L thanks to the adoption of Voigt notation, strongly enforcing the symmetry of the stresstensor, and for the primary variables ρ and ρ υ , regardless of the compressibility coefficient ε . The optimal convergence of the mixed variable and the procedure presented in section2.3 to resolve the underdetermination of the rigid body motions allow the construction ofa superconvergent velocity field υ (cid:63) (converging therefore with order k + 2). A reductionof 0 . k = 1 and ε ∗ = 1, which represents however a limit case, notmeaningful from a physical point of view. The second numerical experiment considers a two dimensional fluid problem with manu-factured solution. The solution aims to exercise all the terms present in the fluid partialdifferential equations, including the time-dependent terms, and to tackle all the nonlinear-ities, i.e., the weak compressibility and the convection. Again, the subscript ( · ) F is omittedhere for brevity. The analytical solution of the problem reads υ x ( x, y, t ) = sin( πx ) sin( πy ) sin( πt ) − π ρ (cid:2) πx ) sin( πy ) cos( πt )+ (cid:0) sin(2 πx ) + +2 πx cos(2 πy ) (cid:1) sin ( πt ) (cid:3) ε, (60a) υ y ( x, y, t ) = cos( πx ) cos( πy ) sin( πt )+ π ρ (cid:2) πx ) cos( πy ) cos( πt ) − (cid:0) sin(2 πy )+2 πy cos(2 πx ) (cid:1) sin ( πt ) (cid:3) ε, (60b) p ( x, y, t ) = π cos( πx ) sin( πy ) sin( πt ) , (60c)26 .01 2.71 (a) L with ε ∗ = 0 . (b) L with ε ∗ = 0 . (c) L with ε ∗ = 1 (d) ρ with ε ∗ = 0 . (e) ρ with ε ∗ = 0 . (f) ρ with ε ∗ = 1 (g) ρ υ with ε ∗ = 0 . (h) ρ υ with ε ∗ = 0 . (i) ρ υ with ε ∗ = 1 (j) υ (cid:63) with ε ∗ = 0 . (k) υ (cid:63) with ε ∗ = 0 . (l) υ (cid:63) with ε ∗ = 1 Figure 6: Spatial convergence of the L -error of the mixed, primal and postprocessed vari-ables for the weakly compressible Poiseuille flow with different values of the dimensionlesscompressibility coefficient. 27nd it has been obtained by adding to a divergence-free velocity field specific O ( ε ) terms,such that the residual of the continuity equation in (6) R ( x, y, t ) = (cid:26) π ρ sin( πx ) sin( πy ) sin( πt ) (cid:2) πx ) sin( πy ) cos( πt )+ (cid:0) sin(2 πx )+2 πx cos(2 πy ) (cid:1) sin ( πt ) (cid:3) + π ρ cos( πx ) cos( πy ) sin( πt ) (cid:2) πx ) cos( πy ) cos( πt ) − (cid:0) sin(2 πy )+2 πy cos(2 πx ) (cid:1) sin ( πt ) (cid:3) + π ρ (cid:0) p − π cos( πx ) sin( πy ) sin( πt ) (cid:1)(cid:2) cos( πx ) sin( πy ) cos( πt )+ (cid:0) cos ( πx )+ cos ( πy ) − (cid:1) sin ( πt ) (cid:3)(cid:27) ε (61)is of second order in ε . A body force to cancel out any imbalance is then added to theright hand side of the momentum equation in (6). If ε = 0, the solution satisfies the fullyincompressible Navier–Stokes equations with no need of adding the term (61).In order to validate the proposed HDG method formulated in the ALE form in section2, two cases are considered:1. the problem is solved on a fixed mesh,2. the problem is solved on a moving mesh, whose displacement is described by a pre-defined function d ( x, y, t ), with d x ( x, y, t ) = 14 sin(2 πx ) (cid:2) − cos(2 πy ) (cid:3)(cid:2) − cos(2 πt ) (cid:3) ¯ d,d y ( x, y, t ) = 14 (cid:2) − cos(2 πx ) (cid:3) sin(2 πy ) (cid:2) − cos(2 πt ) (cid:3) ¯ d, (62)with ¯ d = 0 . a ( x, y, t ) is evaluated at the elementallevel through standard finite differentiation techniques. The fluid domain is the unit squareΩ = [0 , × [0 ,
1] and the same triangular pattern of the previous example is consideredfor the meshes used to perform the convergence studies. Figure 7 shows the third levelof refinement of the fixed (left) and moving (right) mesh at t = 0 .
5. The degrees ofapproximation considered for the convergence studies are k = [1 , , µ isconsidered equal to 0 .
1, while the reference density ρ is taken equal to 1 and evaluatedat the reference pressure p = 0, with a compressibility coefficient ε equal to 0 .
1. Thestabilization parameters are set as τ ρ = 10 /ε and τ ρυ = 1. The final time of the simulationis t = 0 . t = [2 − , − , − ] for28 a) Fixed mesh (b) Moving mesh at t = 0 . Figure 7: Third level of refinement of the fixed (left) and moving (right) mesh at t = 0 . k = [1 , , D = ∂ Ω are computed from the analytical solution (60).In Figure 8 the solution of the density and the momentum field obtained with theproposed HDG formulation using m = 3 and k = 1 (top) and k = 3 (bottom) on thefixed (left) and moving (right) mesh at the final time is shown. With this choice of thecompressibility coefficient, the maximum variation of the density from the reference valueis about ± L norm as a function of the characteristic element size h is presented in Figure 9, for thecase 1 (left) and the case 2 (right). Although the errors in case 2 (evaluated on a distortedmesh) are systematically larger than the corresponding ones in case 1 by about 0 .
5, 1and 2 orders of magnitude for the finest meshes considered for k equal to 1, 3 and 5,respectively, the optimal convergence rates are nicely preserved. The capability of themethod to preserve the optimal convergence on arbitrarily moving meshes is a crucialfeature in order to accurately solve multiphysics problems like fluid-structure interaction.29 a) ρ with k = 1 on fixed mesh (b) ρ with k = 1 on moving mesh(c) | ρ υ | with k = 1 on fixed mesh (d) | ρ υ | with k = 1 on moving mesh(e) ρ with k = 3 on fixed mesh (f) ρ with k = 3 on moving mesh(g) | ρ υ | with k = 3 on fixed mesh (h) | ρ υ | with k = 3 on moving mesh Figure 8: Approximation of the density and the momentum field of the fluid problem withmanufactured solution at t = 0 . .01 3.91 5.91 (a) L on fixed mesh (b) L on moving mesh (c) ρ on fixed mesh (d) ρ on moving mesh (e) ρ υ on fixed mesh (f) ρ υ on moving mesh Figure 9: Spatial convergence of the L -error of the mixed and primal variables evaluatedat t = 0 . .3 FSI problem with manufactured solution The third numerical example aims to verify the convergence properties of the two couplingstrategies presented in section 3 for the solution of fluid-structure interaction problems. Thegeneration of manufactured solutions for the verification of FSI formulations can be a quitecomplex task and a rigorous procedure to generate non-trivial solutions has been developedin [10]. The problem under analysis considers a two dimensional unsteady incompressibleflow interacting with a flexible structure and this flow configuration is simulated by choosinga sufficiently small compressibility coefficient. The compressibility level however does notplay an important role for the coupling techniques we aim to test.The fluid domain is the square Ω F = [0 , × [0 ,
1] while the structural domain is definedas Ω S = [0 , × [1 , . I = { ( x, y ) ∈ R | y = 1 } . The parameters considered to generate thesolution according to the procedure in [10] are: δ ( t ) = δ sin [2 π ( t + t )] ,f ( x, t ) = 1 + δ ( t ) [1 − cos (2 πx )] sin (2 πx ) ,K ( x, t ) = x (cid:2) − h ( t ) /δ ) (cid:3) ,j = 1 ,a ( x, y, t ) = [1 + 10 δ ( t ) cos (2 πx ) /
3] (1 − y ) ,b ( x, y, t ) = 1 , (63)with δ = 0 .
05 and t = 0 .
25. The function f ( x, t ) describes the deformed fluid-structureinterface, K ( x, t ) is a user supplied data characterizing the fluid flow, while the chosenparameter j = 1 provides a non zero fluid velocity profile along the bottom boundary.Moreover, a ( x, y, t ) and b ( x, y, t ) define the structural solution. Two auxiliary functionsare then derived to ease the notation: M ( x, t ) = (cid:90) x K ( z, t ) dz, L ( x, t ) = (cid:90) x K ( z, t ) zdz. (64)The structural displacement is defined as follows: u S x ( x, y, t ) = a ( x, y, t ) ,u S y ( x, y, t ) = b ( x, y, t ) ( f ( x, t ) − . (65)An appropriate body force canceling out the imbalance in the structural PDEs is added inthe right hand side of the equation in (27). The fluid velocity is then constructed as: υ F x ( x, y, t ) = ( j + 1) y j (cid:0) M [ f ( x, t ) , t ] − M ( y, t ) (cid:1) − jy j − (cid:0) L [ f ( x, t ) , t ] − L ( y, t ) (cid:1) ,υ F y ( x, y, t ) = y j (cid:0) f ( x, t ) − y (cid:1) K (cid:2) f ( x, t ) , t (cid:3) ∂f∂x ( x, t ) + ∂f∂t ( x, t ) . (66)The parameters chosen induce a sufficiently complex flow which is shown in Figure 10on the initial (left) and final (right) deformed configuration. The velocity field presents32 a) Initial configuration (b) Final configurationh Figure 10: Velocity field of the FSI problem with manufactured solution on the initial (left)and final (right) deformed configuration.two pronounced vortices circulating in opposite directions and it equals the structuralvelocity on the fluid-structure interface. The structural forces on the deformed interfaceare evaluated through the push-forward operation S ( x, t ) = (cid:2) | F S | − P S F T S n (cid:3)(cid:12)(cid:12)(cid:12) y = f ( x,t ) , (67)where n is the unit normal vector to the interface whose Cartesian components can beevaluated from the definition of the deformed interface f ( x, t ) as n x ( x, t ) = − ∂f∂x ( x, t ) (cid:104) (cid:16) ∂f∂x ( x, t ) (cid:17) (cid:105) − / ,n y ( x, t ) = (cid:104) (cid:16) ∂f∂x ( x, t ) (cid:17) (cid:105) − / . (68)The fluid pressure is defined as p F ( x, t ) = − AS y n x − BS y n y + BS x n x + CS x n y − Bn x − Cn x n y + An x n y + Bn y , (69)where A ( x, t ) = 2 ∂υ F x ∂x (cid:12)(cid:12)(cid:12) y = f ( x,t ) , B ( x, t ) = (cid:16) ∂υ F x ∂y + ∂υ F y ∂x (cid:17)(cid:12)(cid:12)(cid:12) y = f ( x,t ) , C ( x, t ) = 2 ∂υ F y ∂y (cid:12)(cid:12)(cid:12) y = f ( x,t ) . (70)33 a) Undeformed mesh (b) Initial deformed mesh (c) Final deformed mesh Figure 11: Fourth level of refinement of the undeformed (left), initial deformed (center)and final deformed (right) mesh used for the convergence studies of the FSI problem withmanufactured solution. In blue is represented the fluid domain, in red the structuraldomain and in black the fluid-structure interface.In order to ensure the continuity of the fluid and solid forces without introducing additionalterms at the interface (that would require substantial modification for most codes), anappropriate spatially varying viscosity is taken into account: µ F ( x, t ) = S x n y − S y n x − Bn x − Cn x n y + An x n y + Bn y . (71)The last ingredient is the introduction of a fluid body force satisfying the momentumequation in (6). The complexity of algebraic manipulation needed to compute this termis quite challenging and the symbolic math toolbox of MATLAB has been used for thispurpose.In order to test the coupling strategies formulated in section 3, two cases are considered:1. the problem is solved with the partitioned Dirichlet–Neumann coupling presented insection 3.3,2. the problem is solved with the monolithic Nitsche-based coupling presented in section3.4.The convergence studies are performed through uniform mesh refinement of the fluid andsolid domains with triangular elements and considering as degree of approximation k =[1 , , ρ F = 1, p F = 0 and ε F = 10 − . The Young modulus andthe Poisson ratio of the elastic structure are E S = 1 and ν S = 0 .
49. The stabilization34 a) Initial configuration (b) Final configurationh
Figure 12: Approximation of the velocity field of the FSI problem with manufacturedsolution on the initial (left) and final (right) deformed configuration.parameters are τ ρ = 1 /ε F and τ ρυ = 10 and the Nitsche parameter is set as γ = 10 forthe case 2. The time span considered is (0 , . t = 2 − . The initial conditions and the boundary conditions imposedon Γ D = ∂ (Ω F ∪ Ω S ) are computed from the analytical solution developed here.In Figure 12 the solution of the fluid velocity field obtained with the monolithic Nitsche-based coupling using m = 4 and k = 3 on the initial (left) and final (right) deformedconfiguration is shown.The convergence of the error of the HDG fluid solution and the CG structural solutionmeasured in the L norm as a function of the characteristic element size h is representedin Figure 13, for the case 1 (left) and the case 2 (right). Optimal convergence rates areobserved for all the variables in both subdomains. By comparing the plots on the left andthe ones on the right, it can be observed that the two coupling strategies provide almostidentical results. The fourth numerical experiment considers a convergent fluid channel containing a flexiblewall structure attached to its bottom. This example is inspired by [35] and its mainfeature is the strong coupling between the fluid and structural fields, given by their similardensities. The geometry and the boundary conditions of the problem are depicted in Figure14. 35 .81 2.81 3.91 (a) L F with partitionedDirichlet–Neumann cou-pling (b) L F with monolithicNitsche-based coupling (c) ρ F with partitionedDirichlet–Neumann cou-pling (d) ρ F with monolithicNitsche-based coupling (e) ρ υ F with partitionedDirichlet–Neumann cou-pling (f) ρ υ F with monolithicNitsche-based coupling (g) u S with partitionedDirichlet–Neumann cou-pling (h) u S with monolithicNitsche-based coupling Figure 13: Spatial convergence of the L -error of the mixed and primal variables evaluatedat t = 0 .
125 for the FSI problem with manufactured solution solved with the partitionedDirichlet–Neumann coupling (left) and the monolithic Nitsche-based coupling (right). ρ υ F = ρ F = ρ F ρυ F y = 0 ρ υ F = ρ υ F i n .
50 0 .
75 0 .
50 0 . . . . µ F = 0 . kg/ ( m · s ) and the reference density ρ F = 956 kg/m evaluated at the reference pressure p F = 0 N/m . Three different ordersof magnitude are considered for the compressibility coefficient, i.e., ε F = [0 . , . , . kg/ ( N · m ). The structure is modelled as a St. Venant–Kirchhoff material with density ρ S = 1500 kg/m , Young’s modulus E S = 2 . × N/m and Poisson’s ratio ν S = 0 . ρυ F x (0 , y, t ) = y (1 − y ) 1 − cos ( πt/ ρ F ¯ υ F if t ≤ s, y (1 − y ) ρ F ¯ υ F if t > s,ρυ F y (0 , y, t ) = 0 , (72)with ¯ υ F = 0 . m/s , resulting in a Reynolds number of about 100 after the flow hasbeen completely formed. No-slip and free-slip boundary conditions are considered on thebottom and on the top sides of the channel, respectively. At the channel exit, the fluiddensity is imposed equal to its reference value.The fluid domain is discretized with the HDG method and contains 772 triangularelements, while the structural domain is discretized with the CG method and contains 200triangular elements. A boundary layer mesh is constructed near the physical walls. Thetime interval studied is 25 s and the temporal integration is performed by means of theBDF2 method. Three different orders of magnitude of the time step size are considered,i.e., ∆ t = [0 . , . , . s . The stabilization parameters are τ ρ = 10 /ε F and τ ρυ = 1.Two cases are considered:1. the problem is solved with the partitioned Dirichlet–Neumann coupling presented insection 3.3,2. the problem is solved with the monolithic Nitsche-based coupling presented in section3.4.For the first case, the convergence tolerance for the coupling iterations is set to a verysmall value, i.e., η = 10 − . For the second case, the Nitsche parameter is set to γ = 100.A preliminary study, conducted with the monolithic scheme for the intermediate com-pressibility coefficient ( ε F = 0 . kg/ ( N · m )) and the largest time step size (∆ t = 0 . s ), has highlighted the need for adopting a higher-order polynomial degree of approxima-tion. Figure 15 shows how the resulting horizontal displacement of the center top point ofthe flexible wall is severely underestimated when computed with linear elements ( k = 1),whereas k = 2 and k = 3 produce a reliable response with a final structural displacementdiffering by less than 1%. Hence, a degree of approximation k = 2 is considered in thefollowing for both subdomains.The structural displacement obtained with both the partitioned and the monolithicscheme is shown in Figure 16 for the different compressibility coefficients. For the sake of37igure 15: x -displacement component of the center top point of the flexible wall withdifferent polynomial degrees of approximation computed with the monolithic Nitsche-basedcoupling with ε F = 0 . kg/ ( N · m ) and ∆ t = 0 . s . (a) u S x with partitioned Dirichlet–Neumann cou-pling (b) u S x with monolithic Nitsche-based coupling Figure 16: x -displacement component of the center top point of the flexible wall withdifferent values of the compressibility coefficient computed with the partitioned Dirichlet–Neumann coupling (left) and the monolithic Nitsche-based coupling (right) with ∆ t = 0 . s . 38 a) ρ F at t = 5 s with ε F = 0 . kg/ ( N · m ) (b) | ρ υ F | at t = 5 s with ε F = 0 . kg/ ( N · m )(c) ρ F at t = 10 s with ε F = 0 . kg/ ( N · m ) (d) | ρ υ F | at t = 10 s with ε F = 0 . kg/ ( N · m )(e) ρ F at t = 25 s with ε F = 0 . kg/ ( N · m ) (f) | ρ υ F | at t = 25 s with ε F = 0 . kg/ ( N · m ) Figure 17: Approximation of the density and the momentum field of the channel withflexible wall at different time instants with ε F = 0 . kg/ ( N · m ).validation and comparison, the displacement is also computed on the same mesh with afully incompressible solver, based on an equal order stabilized finite element formulation,with the same settings used for the partitioned Dirichlet–Neumann coupling. The resultsare added in the left panel of Figure 16. By comparing the left and the right panels, it canbe observed that the two coupling strategies provide almost identical results. Moreover,for the compressibility levels considered, the physical results are sufficiently close to thoseobtained by considering the flow fully incompressible.The approximation of the density and the momentum field obtained with the monolithicNitsche-based coupling with ε F = 0 . kg/ ( N · m ) is shown in Figure 17 at different timeinstants.The average value of the dynamic relaxation parameter ω avg (evaluated over all the39 t = 0 . s ∆ t = 0 . s ∆ t = 0 . sω avg i avg ω avg i avg ω avg i avg Incompressible 0.17 27.5 – – – – ε F = 0 . kg/ ( N · m ) 0.19 22.5 0.20 34.0 0.45 16.1 ε F = 0 . kg/ ( N · m ) 0.20 22.5 0.29 22.6 0.64 8.2 ε F = 0 . kg/ ( N · m ) 0.23 18.4 0.47 11.4 0.86 4.6Table 1: Average value of the dynamic relaxation parameter (evaluated over all the couplingiterations) and average number of coupling iterations (evaluated over all the time steps)needed by the Dirichlet–Neumann algorithm for the channel with flexible wall.coupling iterations) and the average number of coupling iterations i avg (evaluated over allthe time steps) needed by the Dirichlet–Neumann algorithm are summarized in Table 1.When using the fully incompressible solver, a large number of coupling iterations isrequired to satisfy the convergence criterion (45) for the largest time step size considered.Moreover, for smaller time step sizes, the coupling fails to converge. As expected fromthe analysis in [32], the weakly compressible formulation for the fluid field instead leadsto a smaller number of coupling iterations and to a larger value of the dynamic relaxationparameter, thanks to the reduction of the maximal eigenvalue of the so-called added massoperator. As also shown in [32], these beneficial effects are proportional to the compress-ibility coefficient and more pronounced for small time step sizes. On the one hand, themonolithic scheme, not requiring sub-iterations of the single-field problems, is observed tooutperform the partitioned one for small ε F and large ∆ t . On the other hand, the parti-tioned scheme becomes very competitive for larger ε F and smaller ∆ t . However, since thetwo algorithms have been implemented on distinct platforms, a quantitative comparisonof the computational cost cannot be fairly assessed in the context of the present work andit constitutes a topic of future research. The fifth numerical example considers an elastic wall embedded in a channel flow. Thisexample was set up by [16] and the goal of this study is to show the advantages of the weaklycompressible fluid formulation for the solution of FSI problems in a three dimensionalsetting. The geometry and the boundary conditions of the problem are depicted in Figure18. The dynamic viscosity of the fluid is µ F = 0 .
01 and the reference density ρ F = 1evaluated at the reference pressure p F = 0. Three different orders of magnitude areconsidered for the compressibility coefficient, i.e., ε F = [0 . , . , . ρ S = 1, Young’s modulus E S = 500 and40 υ F = ρ F = ρ F ρ υ F = ρ υ F = ρ υ F i n . . . .
05 0 . . . . . ρ υ F = ρ F = ρ F ρ υ F = ρ υ F = ρ υ F i n . . . .
05 0 . ν S = 0. The following parabolic momentum profile is imposed at the inlet: ρυ F x (0 , y, z, t ) = (cid:0) − y (cid:1) (cid:0) − z (cid:1) − cos ( πt/ ρ F ¯ υ F if t ≤ , (cid:0) − y (cid:1) (cid:0) − z (cid:1) ρ F ¯ υ F if t > ,ρυ F y (0 , y, z, t ) = 0 ,ρυ F z (0 , y, z, t ) = 0 , (73)with ¯ υ F = 0 .
1, resulting in a Reynolds number of about 5. No-slip boundary conditions areapplied on the top and bottom walls as well as on the lateral walls. At the channel exit, thefluid density is imposed equal to its reference value. The HDG fluid discretization contains25664 hexahedral elements, while the CG structural discretization contains 1536 hexahedralelements. The degree of approximation considered is k = 1 for both subdomains. The timeinterval studied is 10 s and the temporal integration is performed by means of the BDF2method. Three different orders of magnitude of the time step size are considered, i.e.,∆ t = [0 . , . , . τ ρ = 1 /ε F and τ ρυ = 1.The problem is solved with the partitioned Dirichlet–Neumann coupling presented insection 3.3, with convergence tolerance η = 10 − .The x component of the displacement of the center top point of the wall is shown inFigure 19 for the different compressibility coefficients considered. The final displacementdiffers from the one obtained with a fully incompressible solver with the same settings forthe coupling algorithm by about 6%.The approximation of the fluid velocity and pressure field and the structural displace-ment obtained with ε F = 0 .
01 and ∆ t = 0 . x -displacement component of the center top point of the flexible wall withdifferent values of the compressibility coefficient computed with the partitioned Dirichlet–Neumann coupling with ∆ t = 0 .
1. ∆ t = 0 . t = 0 .
01 ∆ t = 0 . ω avg i avg ω avg i avg ω avg i avg Incompressible 0.33 10.9 – – – – ε F = 0 .
001 0.37 10.1 0.43 7.7 0.54 2.2 ε F = 0 .
01 0.39 8.7 0.57 4.2 0.74 0.8 ε F = 0 . ω avg (evaluated over all thecoupling iterations) and the average number of coupling iterations i avg (evaluated over allthe time steps) needed by the partitioned coupling algorithm are summarized in Table2. The response obtained with a higher polynomial degree ( k = 2) in the fluid fieldwith ∆ t = 0 . ω avg = [0 . , . , .
44] and i avg = [11 . , . , .
5] for ε F = [0 . , . , . t and it fails to reach convergenceof the partitioned algorithm for smaller time step sizes. On the other hand, the weaklycompressible solver is always able to converge, exhibiting a decreasing number of coupling42 a) | υ F | at t = 10 with ε F = 0 . p F at t = 10 with ε F = 0 . | u S | at t = 10 with ε F = 0 . Figure 20: Approximation of the fluid velocity and pressure field and the structural dis-placement of the 3D channel with flexible wall at t = 10 with ε F = 0 . t . For ε F = 0 . t = 0 . ω = 1) and almost no coupling iterations are needed to meet the convergencecriterion. A hybrid HDG-CG formulation for the solution of fluid-structure interaction problems hasbeen proposed. Special emphasis has been devoted to the derivation and validation of anovel hybridizable discontinuous Galerkin approach for weakly compressible flows with theArbitrary Lagrangian–Eulerian description. A partitioned Dirichlet–Neumann scheme hasbeen revisited for the hybrid discretization and a minimally-intrusive monolithic Nitsche-based coupling has been proposed for the coupling of the HDG and the CG subproblems,by exploiting the definition of the numerical flux and the trace of the solution to impose thecoupling conditions. The numerical examples demonstrate the convergence of the HDGand CG primal and mixed variables with order k + 1 and the superconvergence of thefluid velocity with order k + 2, through an inexpensive element-by-element postprocessing.Moreover, the advantages of introducing a weak compressibility in the fluid have beenconfirmed on two and three dimensional FSI problems. In particular, the constraints thatthe incompressibility poses on the fluid-structure coupling are alleviated by the proposedformulation, resulting in a more robust and efficient FSI solver. The investigation ofthe computational efficiency of the proposed schemes in a high-performance framework issubject of future work. Acknowledgments
This work was supported by the European Education, Audiovisual and Culture ExecutiveAgency (EACEA) under the Erasmus Mundus Joint Doctorate “Simulation in Engineeringand Entrepreneurship Development” (SEED), FPA 2013-0043.
References [1] F. Bassi, A. Crivellini, D. Di Pietro, and S. Rebay. An artificial compressibility fluxfor the discontinuous galerkin solution of the incompressible navier–stokes equations.
Journal of Computational Physics , 218(2):794–815, 2006.[2] F. Bassi, A. Crivellini, D. A. Di Pietro, and S. Rebay. An implicit high-order discon-tinuous galerkin method for steady and unsteady incompressible flows.
Computers &Fluids , 36(10):1529–1546, 2007. Special Issue Dedicated to Professor Michele Napoli-tano on the Occasion of his 60th Birthday.443] P. Causin, J. Gerbeau, and F. Nobile. Added-mass effect in the design of partitionedalgorithms for fluid–structure problems.
Computer Methods in Applied Mechanics andEngineering , 194(42):4506 – 4527, 2005.[4] J. Chung and G. M. Hulbert. A time integration algorithm for structural dynamicswith improved numerical dissipation: The generalized- method.
Journal of AppliedMechanics , 60(2):371–375, 1993.[5] B. Cockburn, B. Dong, and J. Guzm´an. A superconvergent ldg-hybridizablegalerkin method for second-order elliptic problems.
Mathematics of Computation ,77(264):1887–1916, 2008.[6] B. Cockburn, J. Gopalakrishnan, and R. Lazarov. Unified hybridization of discon-tinuous galerkin, mixed, and continuous galerkin methods for second order ellipticproblems.
SIAM Journal on Numerical Analysis , 47(2):1319–1365, 2009.[7] B. Cockburn and K. Shi. Superconvergent HDG methods for linear elasticity withweakly symmetric stresses.
IMA Journal of Numerical Analysis , 33(3):747–770, 102012.[8] J. Degroote, K.-J. Bathe, and J. Vierendeels. Performance of a new partitioned pro-cedure versus a monolithic procedure in fluid-structure interaction.
Comput. Struct. ,87(11-12):793–801, June 2009.[9] Y. Epshteyn and B. Rivi`ere. Estimation of penalty parameters for symmetric inte-rior penalty galerkin methods.
Journal of Computational and Applied Mathematics ,206(2):843–872, 2007.[10] S. ´Etienne, A. Garon, and D. Pelletier. Some manufactured solutions for verificationof fluid-structure interaction codes.
Computers & Structures , 106-107:56 – 67, 2012.[11] C. Farhat, K. G. van der Zee, and P. Geuzaine. Provably second-order time-accurateloosely-coupled solution algorithms for transient nonlinear computational aeroelastic-ity.
Computer Methods in Applied Mechanics and Engineering , 195(17):1973 – 2001,2006. Fluid-Structure Interaction.[12] N. Fehn, W. A. Wall, and M. Kronbichler. A matrix-free high-order discontinuousgalerkin compressible navier-stokes solver: A performance comparison of compress-ible and incompressible formulations for turbulent incompressible flows.
InternationalJournal for Numerical Methods in Fluids , 89(3):71–102, 2019.[13] J. Fish and T. Belytschko.
A First Course in Finite Elements . John Wiley & Sons,2007.[14] C. F¨orster, W. A. Wall, and E. Ramm. Artificial added mass instabilities in sequentialstaggered coupling of nonlinear structures and incompressible viscous flows.
ComputerMethods in Applied Mechanics and Engineering , 196(7):1278 – 1293, 2007.4515] B. Froehle and P.-O. Persson. A high-order discontinuous galerkin method for fluid–structure interaction with efficient implicit–explicit time stepping.
Journal of Com-putational Physics , 272:455 – 470, 2014.[16] A. Gerstenberger and W. A. Wall. An embedded dirichlet formulation for 3d continua.
International Journal for Numerical Methods in Engineering , 82(5):537–563, 2010.[17] M. Giacomini, A. Karkoulias, R. Sevilla, and A. Huerta. A superconvergent hdgmethod for stokes flow with strongly enforced symmetry of the stress tensor.
Journalof Scientific Computing , 77(3):1679–1702, Dec 2018.[18] M. Giacomini, R. Sevilla, and A. Huerta.
Tutorial on Hybridizable Discontinu-ous Galerkin (HDG) Formulation for Incompressible Flow Problems , pages 163–201.Springer International Publishing, Cham, 2020.[19] G. Giorgiani, S. Fern´andez-M´endez, and A. Huerta. Hybridizable discontinuousgalerkin with degree adaptivity for the incompressible navier–stokes equations.
Com-puters & Fluids , 98:196–208, 2014.[20] G. Giorgio, F.-M. Sonia, and H. Antonio. Hybridizable discontinuous galerkin p-adaptivity for wave propagation problems.
International Journal for Numerical Meth-ods in Fluids , 72(12):1244–1262, 2018/07/05 2013.[21] M. Griebel and M. A. Schweitzer.
A Particle-Partition of Unity Method Part V:Boundary Conditions , pages 519–542. Springer Berlin Heidelberg, Berlin, Heidelberg,2003.[22] P. Hansbo. Nitsche’s method for interface problems in computational mechanics.
GAMM-Mitteilungen , 28(2):183–206, 2020/07/24 2005.[23] M. Heil. An efficient solver for the fully coupled solution of large-displacement fluid–structure interaction problems.
Computer Methods in Applied Mechanics and Engi-neering , 193(1):1 – 23, 2004.[24] K. D. Housiadas and G. C. Georgiou. New analytical solutions for weakly compressiblenewtonian poiseuille flows with pressure-dependent viscosity.
International Journal ofEngineering Science , 107:13–27, 2016.[25] D. Huang, P.-O. Persson, and M. Zahr. High-order, linearly stable, partitioned solversfor general multiphysics problems based on implicit–explicit runge–kutta schemes.
Computer Methods in Applied Mechanics and Engineering , 346:674–706, 2019.[26] K. E. Jansen, C. H. Whiting, and G. M. Hulbert. A generalized- method for integratingthe filtered navier–stokes equations with a stabilized finite element method.
ComputerMethods in Applied Mechanics and Engineering , 190(3):305 – 319, 2000.4627] D. Jean, H. Antonio, P. J.-Ph., and R.-F. A.
Arbitrary Lagrangian–Eulerian Methods ,chapter 14. American Cancer Society, 2004.[28] R. M. Kirby, S. J. Sherwin, and B. Cockburn. To cg or to hdg: A comparative study.
Journal of Scientific Computing , 51(1):183–212, 2012.[29] B. Krank, N. Fehn, W. A. Wall, and M. Kronbichler. A high-order semi-explicitdiscontinuous galerkin solver for 3d incompressible flow with application to dns andles of turbulent channel flow.
Journal of Computational Physics , 348:634 – 659, 2017.[30] U. K¨uttler, C. F¨orster, and W. A. Wall. A solution for the incompressibility dilemmain partitioned fluid–structure interaction with pure dirichlet fluid domains.
Computa-tional Mechanics , 38(4):417–429, Sep 2006.[31] U. K¨uttler and W. A. Wall. Fixed-point fluid–structure interaction solvers with dy-namic relaxation.
Computational Mechanics , 43(1):61–72, Dec 2008.[32] A. La Spina, C. F¨orster, M. Kronbichler, and W. A. Wall. On the role of (weak) com-pressibility for fluid-structure interaction solvers.
International Journal for NumericalMethods in Fluids , 92(2):129–147, 2020.[33] A. La Spina, M. Giacomini, and A. Huerta. Hybrid coupling of cg and hdg discretiza-tions based on nitsche’s method.
Computational Mechanics , 65(2):311–330, Feb 2020.[34] M. Mayr, T. Kl¨oppel, W. A. Wall, and M. W. Gee. A temporal consistent mono-lithic approach to fluid-structure interaction enabling single field predictors.
SIAM J.Scientific Computing , 37, 2015.[35] D. Mok and W. Wall. Partitioned analysis schemes for the transient interaction ofincompressible flows and nonlinear flexible structures.
Trends in Computational Struc-tural Mechanics , 01 2001.[36] A. Montlaur, S. Fernandez-Mendez, and A. Huerta. Discontinuous galerkin methodsfor the stokes equations using divergence-free approximations.
International Journalfor Numerical Methods in Fluids , 57(9):1071–1092, 2008.[37] N. C. Nguyen, J. Peraire, and B. Cockburn. An implicit high-order hybridizablediscontinuous galerkin method for linear convection–diffusion equations.
Journal ofComputational Physics , 228(9):3232–3254, 2009.[38] N. C. Nguyen, J. Peraire, and B. Cockburn. An implicit high-order hybridizablediscontinuous galerkin method for nonlinear convection–diffusion equations.
Journalof Computational Physics , 228(23):8841–8855, 2009.[39] N. C. Nguyen, J. Peraire, and B. Cockburn. A hybridizable discontinuous galerkinmethod for stokes flow.
Computer Methods in Applied Mechanics and Engineering ,199(9):582–597, 2010. 4740] N. C. Nguyen, J. Peraire, and B. Cockburn. An implicit high-order hybridizablediscontinuous galerkin method for the incompressible navier–stokes equations.
Journalof Computational Physics , 230(4):1147–1170, 2011.[41] M. Paipuri, C. Tiago, and S. Fern´andez-M´endez. Coupling of continuous and hybridiz-able discontinuous galerkin methods: Application to conjugate heat transfer problem.
Journal of Scientific Computing , 78(1):321–350, Jan 2019.[42] J. Peraire, N. Nguyen, and B. Cockburn.
A Hybridizable Discontinuous GalerkinMethod for the Compressible Euler and Navier-Stokes Equations . American Instituteof Aeronautics and Astronautics, 2018/07/05 2010.[43] P.-O. Persson, J. Bonet, and J. Peraire. Discontinuous galerkin solution of the navier–stokes equations on deformable domains.
Computer Methods in Applied Mechanicsand Engineering , 198(17):1585 – 1595, 2009.[44] B. Schott, C. Ager, and W. A. Wall. A monolithic approach to fluid-structure interac-tion based on a hybrid eulerian-ale fluid domain decomposition involving cut elements.
CoRR , abs/1808.00343, 2018.[45] R. Sevilla, M. Giacomini, A. Karkoulias, and A. Huerta. A superconvergent hybridis-able discontinuous galerkin method for linear elasticity.
International Journal forNumerical Methods in Engineering , 116(2):91–116, 2018.[46] R. Sevilla and A. Huerta. Hdg-nefem with degree adaptivity for stokes flows.
Journalof Scientific Computing , 77(3):1953–1980, 2018.[47] J. P. Sheldon, S. T. Miller, and J. S. Pitt. A hybridizable discontinuous galerkinmethod for modeling fluid–structure interaction.
Journal of Computational Physics ,326:91 – 114, 2016.[48] J. P. Sheldon, S. T. Miller, and J. S. Pitt. An improved formulation for hybridizablediscontinuous galerkin fluid-structure interaction modeling with reduced computa-tional expense.
Communications in Computational Physics , 24(5), 6 2018.[49] S. C. Soon, B. Cockburn, and H. K. Stolarski. A hybridizable discontinuous galerkinmethod for linear elasticity.
International Journal for Numerical Methods in Engi-neering , 80(8):1058–1092, 2018/07/05 2009.[50] E. Taliadorou, G. C. Georgiou, and E. Mitsoulis. Numerical simulation of the extrusionof strongly compressible newtonian liquids.
Rheologica Acta , 47(1):49–62, Jan 2008.[51] D. C. Venerus. Laminar capillary flow of compressible viscous fluids.
Journal of FluidMechanics , 555:59–80, 2006. 4852] G. Vinay, A. Wachs, and J.-F. Agassant. Numerical simulation of weakly compressiblebingham flows: The restart of pipeline flows of waxy crude oils.
Journal of Non-Newtonian Fluid Mechanics , 136(2):93 – 105, 2006.[53] W. A. Wall, A. Gerstenberger, P. Gamnitzer, C. F¨orster, and E. Ramm. Large de-formation fluid-structure interaction – advances in ale methods and new fixed gridapproaches. In H.-J. Bungartz and M. Sch¨afer, editors,
Fluid-Structure Interaction ,pages 195–232, Berlin, Heidelberg, 2006. Springer Berlin Heidelberg.[54] O. Zienkiewicz and R. Taylor.