A smoothed finite element approach for computational fluid dynamics: applications to incompressible flows and fluid–structure interaction

In this paper the cell-based smoothed finite element method (CS-FEM) is introduced into two mainstream aspects of computational fluid dynamics: incompressible flows and fluid–structure interaction (FSI). The emphasis is placed on the fluid gradient smoothing which simply requires equal numbers of Gaussian points and smoothing cells in each four-node quadrilateral element. The second-order, smoothed characteristic-based split scheme in conjunction with a pressure stabilization is then presented to settle the incompressible Navier–Stokes equations. As for FSI, CS-FEM is applied to the geometrically nonlinear solid as usual. Following an efficient mesh deformation strategy, block-Gauss–Seidel procedure is adopted to couple all individual fields under the arbitrary Lagriangian–Eulerian description. The proposed solvers are carefully validated against the previously published data for several benchmarks, revealing visible improvements in computed results.


Introduction
Computational fluid dynamics (CFD) is a modern discipline concerned with mathematical modeling, numerical methods and software tools of fluid dynamics. The principle of the majority of CFD problems relies on the Navier-Stokes (NS) equations which describe how different variants (e.g. the velocity, pressure, temperature and density) of a moving fluid are related and predict what is the flow state physically. As a matter of fact, the numerical resolution B Tao He taohe@shnu.edu.cn 1 of these balance equations (with moving and deformable boundaries) has drawn a substantial amount of endeavors from both scientific and engineering communities. Under the umbrella of this major branch, finite element method (FEM) is widely utilized for analysis and design of a variety of specializations like incompressible flows past bluff/streamlined bodies [63,78] and fluid-structure interaction (FSI) [4,58]. Normally, the standard Galerkin finite element procedure for incompressible flows is confronted with two sources of numerical instabilities [63]. One source is due to the presence of the convective acceleration in the NS equations, mainly causing spurious oscillations in the velocity field. The other source rests with the inappropriate usage of a pair of interpolation functions for velocity and pressure fields, which primarily poses the pressure oscillations. In the past decades, a number of stabilized FEMs have been successfully devised to prevent these potential instabilities. Popular approaches contain streamline upwind/Petrov-Galerkin (SUPG) formulation [8], Taylor-Galerkin method [18], Galerkin/least-squares (GLS) technique [35], pressurestabilized Petrov-Galerkin (PSPG) formulation [66], spacetime FEM [57,59,60,64,65], characteristic-based split (CBS) scheme [53,77], etc. Nowadays FEM has become one of stable and robust numerical methodologies for solving CFD related problems and delivering deep insights into fluid physics. For example, the complex three-dimensional turbocharger flow is handled by the space-time variational multiscale FEM incorporating isogeometric analysis [60].
Typical finite element solution results in a system of algebraic equations, which approximates the original partial differential equations (PDEs) based upon a finite element discretization. During the process, we probably observe the overly-stiff phenomenon owing to the fully compatible strain field [49]. In examining meshless and finite element methods, Liu and his colleagues [48] found gradient smoothing [9,71] an elegant remedy to the afore-mentioned overly-stiff issue, as well as a valuable alliance of these two methods.
In the seminal publication [48], the smoothed finite element method (SFEM) is proposed by incorporating gradient smoothing operation with the traditional FEM. The essential idea behind SFEM lies in modification of the compatible strain field whereby a Galerkin model may deliver some superior properties. This method is saliently featured with the softened stiffness matrix that yields more accurate solution to discrete PDEs than the standard FEM. After a decade of development, a group of SFEM models has been fostered with versatile applications in solid mechanics, heat transfer and acoustics whose governing equations perfectly suit the technique after introducing divergence theorem. The reader is referred to the monograph [49] and the review article [72] for the state of the art of SFEM.
In recent years, SFEM has seemingly been applied into a range of FSI problems as follows. The immersed SFEM is developed in [70,[73][74][75] where SFEM is responsible for nonlinear solids. Wang et al. [68] integrated SFEM in solid mechanics with a strong-form fluid solver under the arbitrary Lagrangian-Eulerian (ALE) description [34]. Similarly, the authors adopted SFEM to quantify the structural finite deformation triggered by fluidic excitation in a partitioned manner [24,27,28,30]. However, these scenarios do not provide any settlements tailored for the NS equations, but rather replicate SFEM's early success in solid mechanics.
The major dilemma that SFEM faces in CFD stems from interpolating the mixed product of a quantity and its gradient in the NS equations, such as the convective acceleration. For this reason, the underlying investments may be discouraged from CFD. Most recently, we have witnessed a joyful breakthrough that makes the cell-based smoothed FEM (CS-FEM) accessible to incompressible NS equations in terms of three schemes to interpolate nodal quantities involving the mixed product [38]. This work is motivated by the need for a broadening utilization of SFEM for CFD. A natural preference is given to the simplest CS-FEM that is initiated on the basis of bilinear four-node quadrilateral (Q4) element. The easy smoothing treatment is proposed for the fluid equations to bypass the strenuous operations. For this purpose, we set up the equal amounts of Gaussian points (GPs) and smoothing cells (SCs) in each element to compute fluid fluxes. That is to say, the contribution from all SCs is accumulated in the loop circulating the Gaussian integration. As a result of this collective effort, a second-order smoothed CBS (S-CBS) scheme with the stabilized pressure gradient projection (SPGP) technique [1,6,12,54] is utilized to decouple the fluid velocity and pressure. As to FSI, the nonlinear block-Gauss-Seidel procedure [25,26] is preferred to interconnect individual fields owing to its attractive simplicity with good convergence. In short, the marriage of SFEM and fractional-step method may soothe the pressure fluctuation on (dynamic) boundaries. This is potentially important to both incompressible flows and FSI simulations.
The remainder of this paper is organized as follows. The theory of CS-FEM is briefly recalled in Sect. 2. The ALE form of fluid governing equations is given in Sect. 3 while the structural dynamics is depicted in Sect. 4. The mesh updating method is described in Sect. 5. Subsequently, Sect. 6 explains the partitioned coupling algorithm in detail. Several benchmark examples are investigated in Sect. 7. Concluding remarks are drawn in the final section.
where n c is the number of SCs in the ith element. As illustrated in Fig. 1, the smoothed gradient of a field variable b at a point x c in an SC is approximated by where ∇ means the gradient operator and ∇ is its smoothed counterpart, Ω designates the SC and W is the Heavisidetype kernel that possesses the following properties [71] W (x − x c ) 0 and Applying Gauss theorem into the right-hand side of Eq. (1) yields where Γ is the boundary of Ω and n is the unit outward normal of Γ . A piecewise constant kernel W is now formulated in form of where A c = Ω dΩ is the area of the SC. Substituting Eq. (4) into Eq. (3), we have where the gradient of a constant vanishes automatically. The Galerkin procedure leads to the following approxi- where N I is the shape function at node I , the bar indicates a nodal quantity and Einstein summation convention is applied. With the aid of Eq. (6), Eq. (5) is rewritten as Since one-point Gaussian quadrature is sufficiently accurate for two-node line integral, the item enclosed within external brackets on the right-hand side of Eq. (7) can be transformed to its algebraic form where 4 indicates the number of segments per quadrilateral SC, x gp i is the GP on the ith segment Γ i and l i is the length of Γ i . By now, no coordinate transformation is involved in the gradient smoothing process such that a heavily distorted mesh is possibly accommodated [48]. Furthermore, the imported smoothing concept may pass on to wet boundaries the improved traction or pressure. The construction of shape functions for CS-FEM is shown in Fig. 2. A Q4 element is partitioned into four quadrilateral SCs in consideration of the stability condition [48]. Of total nine nodes, extra five nodes are generated to compute the smoothed shape functions by simply averaging those values at four corners [14,48].

Governing equations
Without loss of generality, the NS equations governing an isothermal incompressible viscous fluid flow on a timedependent domain Ω f ⊂ R 2 in a time interval (0, T ) are written in their ALE formulation of where u is the fluid velocity, ρ f is the fluid density, c = u−w is the convective velocity, w is the mesh velocity, f f is the fluid body force and σ f is the fluid stress tensor. For the Eulerian flows, we have w = 0 such that c degenerates into u.
The constitutive equation for Newtonian fluid reads as where p is the pressure, I denotes the identity tensor, μ is the dynamic viscosity, indicates the rate-of-strain tensor and superscript T means transpose. It is assumed that proper boundary conditions (abbreviated to BCs) are imposed at different segments of the domain boundary Γ f below where Γ f d and Γ f n designate the Dirichlet and Neumann subboundaries, respectively, and n f is the unit outward normal of Γ f n . The fluid problem is initiated by prescribing initial conditions as follows where x and t, of course, represent the spatial and temporal coordinates, respectively. The coupling conditions on fluid-structure interface Σ will be presented in a separate subsection.
In view of the characteristic length L and the free-stream velocity U , we define the dimensionless scaleŝ along with the constitutive relation where Re = ρ f U L/μ is the Reynolds number and all hats are dropped. The nondimensionalized BCs and initial conditions share the same form as Eqs. (12) and (13).

Solution procedure
The CBS scheme [53,77] combines the characteristic Galerkin method [50] with the fractional-step method [10,62]. The former process suppresses spurious oscillations via higher-order time stepping in the convection-dominated flows whereas the latter procedure stabilizes the pressure field. The second-order pressure splitting error can be guaranteed by inclusion of the pressure gradient, but the increased accuracy inevitably imperils the stabilizing properties of the original first-order scheme [12,54]. In what follows, the additional SPGP stabilization [12] is introduced to overcome this penalty.
An auxiliary equation in regard to the variable q is defined, whereby the continuity equation (14) is modified as with φ denoting the stabilization parameter. We will numerically discuss φ later as it is yet unclear how to exactly determine the parameter in theory [12]. The temporal discrete version of Eqs. (18), (15) and (17) may be written as where superscript n denotes the nth time slice and Δt = t n+1 −t n is the time step. The auxiliary variable q is explicitly treated in Eq. (19) while the pressue of Eq. (20) is temporally discretized in the semi-implicit manner. Following the CBS procedure, Eq. (20) admits the decomposition below where u * signifies the intermediate velocity, and the body force and the third-order terms are neglected. Taking the divergence of Eq. (23) and expanding the semidiscrete form of Eq. (19) at the next time level yield where the third-order terms are discarded as well.
With introduction of the gradient smoothing, the main steps of the stabilized second-order S-CBS scheme are arranged below Step 1: Predict the velocity field Step 2: Update the pressure field Step 3: Correct the velocity field Step 4: Renew the auxiliary variable Imposition of BCs is straightforward: the velocity BCs are prescribed for Steps 1 and 3 on Γ f u , while the pressure BC for Step 2 on Γ f p . Besides, q 2 = 0 is applied on the pressure-free outlet for Step 4.

Time-step limitations
It is of interest to remark that the semi-implicit CBS scheme is conditionally stable [78]. The general time-step limitations are recommended as [53,78] Δt Δt crit = min(Δt conv , Δt diff ), (29) where Δt crit signifies the critical time step, and Δt conv and Δt diff are the convection and diffusion limits, respectively. The latter two velocities are calculated from where h means the characteristic size of the element.
To account for the stability and convergence, Codina [12] advocated that the pressure stabilization parameter must satisfy the following relationship for viscous dominated flows on the Eulerian mesh. By inspecting Eq. (30), we may reconsider the inequality as Recalling Eq. (29), the range of the stabilization parameter is suggested as or, it may be replaced for safety by

Finite element discretization
The standard Galerkin spatial approximation is performed to discretize the fluid equations in space. Since the CBS scheme permits the equal-/low-order interpolation for both velocity and pressure, the two primitive variables are approximated as where N designates the shape function of Q4 element. Substitution of spatial approximation (35) into the semi-discrete form of Eqs. (25)-(28) entails the final matrix form where the assembled matrices and vectors are presented below At a closer observation of the above representation, the two smoothed element matrices being derivated from the second derivatives are directly finalized by assembly of all SCs of the element e, i.e., The remaining smoothed element matrices manifest themselves in the mixed product of a quantity and its first derivative, see e G for example. To handle those items, we simply dictate that the number and numbering of GPs per Q4 element exactly equal to those of SCs. Since 2 × 2 GPs in Q4-FEM and n c = 4 in CS-FEM are often adopted in practice, it is straightforward to estimate e G below As can be seen from Fig. 3, the contribution of four SCs and GPs per Q4 element successfully circulates within one recurrence. Compared to [38], we organize the smoothed element integral into a more comprehensible pattern.

Structural dynamics
Consider a structural domain Ω s ⊂ R 2 with the boundary Γ s which comprises the same three types of boundaries as well.
A structure immersed in a fluid continuously sustains the fluctuating fluid force. The equation of motion is expressed in the Lagrangian description with proper initial and boundary conditions. The isotropic assumption is made for the structural problem.

Rigid-body motion
In the case of a single rigid body undergoing both translation and rotation (see Fig. 4), the structural displacement is represented by d = {d 1 , d 2 , θ} T where subscripts 1, 2 and θ designates the horizontal, vertical and rotational components defined at the center of gravity G, respectively. The equation of structural motion is formulated by  where the dot illuminates the derivative with respect to t, m i , c i and k i (i = 1, 2 and θ ) stand for the generalized mass, damping and stiffness of the structure, R = {F d , F l , F m } T is the applied fluid force, F d , F l and F m signify the drag, lift and pitching moment, respectively. As pictured in Fig. 4, the compatibility condition must be satisfied between the surface point P and the center of gravity G [55]. Next, the dimensionless scaleŝ and the reduced parameters are computed to nondimensionalize Eq. (42), where the drag coefficient C d , the life coefficient C l and the moment coefficient C m are the dimensionless applied forces, the mass ratiô m i is the dimensionless mass, ξ i is the damping ratio, f ri is the reduced natural frequency, and f ni is the natural frequency. Therefore, the dimensionless equation of structural motion is visualized as

Flexible-body motion
For an elastic solid, the elastodynamics equation governing the conservation law of linear momentum reads as where ρ s is the structural density, f s is the structural body force, σ s is the Cauchy stress tensor and the structural damping is omitted. Other material constants contain Young's modulus E and Poisson's ratio ν. The plane stress assumption is made for the two-dimensional case.
To accommodate the geometrical nonlinearity, the Saint Venant-Kirchhoff constitutive model is assumed below where S is the second Piola-Kirchhoff stress tensor, D stands for the constitutive tensor, E means the Green-Lagrangian strain tensor, and F = I + ∇d is the deformation gradient tensor. The second Piola-Kirchhoff stress tensor, S, is related to the Cauchy stress tensor, σ s , via the geometric transformation given by where J = det(F). The initial and boundary conditions are imposed to close the system of solid equations in the following manner where n s is the unit outward normal of Γ s n . Likewise, the following dimensionless scales are defined in order to enable the nondimensionalization of Eq. (44). Discarding all superscript hats, the dimensionless version of the geometrically nonlinear elastodynamics equation is established as alongside with the given initial and boundary conditions.

Finite element discretization
Here we commence spatial discretization for the elastic solid. As usual, the standard Galerkin procedure is used with the finite element approximation to the displacement, velocity and acceleration which generates the incremental equilibrium equation for dynamic analysis below where K represents the tangent stiffness matrix, M s is the mass matrix, Δd =d n+1 −d n is the increment of nodal displacement, R is the external force and P is the internal force.
Depending upon the geometrical nonlinearity, it is necessary to iterate Eq. (50) in each load step until a required tolerance is satisfied. This linearization is carried out by the modified Newton-Raphson procedure using total Lagrangian formulation [2]. The mass of the body considered is assumed to be conserved in dynamic analysis. Hence the smoothed equilibrium iteration equation is written as where δd (k) is the incremental displacement in the kth subiteration at the current time step and the tangent stiffness matrix is decomposed into linear and nonlinear parts, namely K = K l + K nl . The resultant matrices and vectors admit the following representation The key to compute these quantities consists in the smoothed deformation gradient tensor F = I + ∇d [13,14]. Details of the modified Newton-Raphson procedure considering specific time discretization methods can be found in [2,7].

Time marching method
The widespread availability of step-by-step time integration algorithms is seen in computational analyses of structural dynamics. Here, the structural movement is integrated in time with the Generalized-α method [11] which is generally superior to the Newmark-β method [52]. To do this, the semi-discrete equation of motion is applied to a general midpoint within one time interval, implying that the following modified equation holds where M s , C and K represent the mass, damping and stiffness matrices, respectively, and we prescribë To setd n+1 as the single unknowns in Eq. (52), the Newmark approximations [52] to the acceleration and velocity at t n+1 are stated as Accordingly, the generalized midpoint acceleration and velocity are given by The time integration parameters β, γ , α m and α f are defined as functions of the spectral radius ρ ∞ [11], whose optimal expressions take the form of where 0 ρ ∞ 1 for the desired level of numerical dissipation. Here we specify ρ ∞ = 0.1 for the rigid body [16] whereas ρ ∞ = 0.5 for the elastic solid [17].
In addition, the calculation of smoothed internal force complies with the interpretation of [43] while working on the elastic solid.

Two-level mesh updating
Imposition of interface conditions in time requires that the position of moving interface is accurately captured in the ALE domain whilst maintaining the satisfactory mesh quality. Hence the mesh deformation is of cardinal significance in fluid-structure coupling. For instance, in the space-time FEM the variational formulation written over its space-time domain automatically takes into account the deformation of the spatial domain with respect to time. Such a process is particularly effective for forced motion of a cylinder where the mesh movement is known a priori [64]. For free motion of a body, a general pseudo-elasticity equation approach is proposed in association with the stabilized space-time FEM [39].
The present mesh deformation method adopts a blend of moving submesh approach (MSA) [45] and the ortho-semi-  torsional spring analogy model [51] in the ALE context. Its fundamental principle comprises two stages below -Spring analogy method assimilates the triangle submesh to the structural motion;   -MSA creates a mapping between the submesh's deformation and that of ALE mesh.
Interested readers are recommended to refer to [24,[31][32][33] for thorough implementation. Though MSA moves fluid nodes with the aid of a background mesh, this technique can reduce the expenditure on spring analogy method while preserving the mesh topology [30,33].
On the other hand, the midpoint rule is applied to the mesh velocity scheme as it automatically meets geometric conservation law for two-dimensional stabilized FEM [46] and outstrips the second-order differencing scheme [19].

Interface coupling conditions
In the partitioned scheme, the interplay between the fluid and structure is accomplished via separately enforcing the velocity continuity and traction equilibrium on Σ as follows where t f = σ f · n s and t s = σ s · n s are the fluid and structural tractions, respectively, n s represents the unit outward normal of Σ pointing from the structure to the fluid and n f = −n s . Since the external force acting on the immersed rigid body is a concentrated load vector, the stress equilibrium on Σ becomes where Δx is the distance between surface point P and center of gravity G, as shown in Fig. 4. Also, the geometrical continuity is supplemented thanks to the mesh movement Moreover, interface conditions (60)-(62) may be recast in a hybrid way to alleviate the adverse time-lag effect [29,30,36].

Block-Gauss-Seidel coupling algorithm
The FSI system constitutes a coupled set of nonlinear algebraic equations to be solved for each time step. For numerical stability, kinematic and kinetic compatibilities are compulsively imposed on Σ through block-Gauss-Seidel procedure which implicitly couples all interacting fields. Extra acceleration technique like the Aitken's Δ 2 method [44] may be adopted for faster convergence. Within one time interval, the present coupling algorithm is elaborated hereinafter.
Step 1: Initialize all variables and set k = 0 Step 2: Extrapolate the interfacē Fig. 10 Vorticity contour of the rigid circular cylinder at Re = 100 Fig. 11 Sketch of geometry and boundary conditions for the freely oscillating circular cylinder Step 3: Start fixed-point iterations and set k ← k + 1 Step 4: Rearrange the fluid mesh Ω f n+1(k) Step 5: Calculate the mesh velocity Step 6: Derive other geometrical quantities if necessary Step 7: Compute the intermediate velocity u * − u n = Δt −c n · ∇u n − ∇ p n + 1 Re ∇ 2 u n + Δt 2 c n · ∇(c n · ∇u n + ∇ p n ) Step 8: Update the pressure Step 9: Correct the velocity Step 10: Renew the auxiliary variable q n+1(k) = ∇ p n+1(k) Step 11: Deduce the fluid load and pass it to the structure/solid Step 12: Solve equation of the structural equation Δtd n − α f Kd n Step 13: Estimate the interfacial residuals Step 14: Check the convergence and the maximum number of subiterations: if not convergent, then go ahead; otherwise, proceed to the next time step Step 15: Relax the position of the interfacē Step 16: Return to Step 3 The stop criterion at the kth subiteration is simply judged with max(g n(k)  where n fs is the number of nodes on the interface, the convergence tolerance is ε = 1.0 × 10 −6 and k max = 200 is the user-defined constant that controls the maximum subiterations at each time step.
Alternatively, a variant algorithm may be acquired in case that Eq. (18) is discretized in time by

Steady cavity flow
The geometry of lid-driven cavity flow is defined in Fig. 5a. The cavity is meshed with 40 × 40 Q4 elements in Fig. 5b. Re = 100 and Δt = 1.0 × 10 −2 are chosen for this problem. The velocity components computed without the SPGP technique are severely oscillatory in Fig. 6. This is because the pressure difference p n+1 − p n will approach to zero in the CBS scheme once steady state is reached. Furthermore, in accordance with [1,12], the stability of the second-order scheme seems more sensitive to a smaller time step. Fig. 7 exhibits no oscillations at steady state since the difference between the Laplacian of p and the divergence of q multiplied by φ stabilizes the pressure variation. Besides, the curve obtained from a smaller φ is closer to [21]. Nithiarasu and Zienkiewicz [54] explained that modifying φ possibly reduces numerical oscillations but could incur accuracy deterioration elsewhere. Among all φ in Table 1, φ = 0.25Δt

Unsteady flow over a circular cylinder
The incompressible flow past a circular cylinder is attempted at Re = 100. The problem definition is plotted in Fig. 8a whereas the finite element discretization is composed of 5190 Q4 elements and 5341 nodes in Fig. 8b. The time step is set as Δt = 1.0 × 10 −2 . Table 2 lists the mean value of drag coefficient C d,mean , the root-mean-square error (RMSE) of drag coefficient C d,rmse , the amplitude of lift coefficient C l,max , the RMSE of lift coefficient C l,rmse and the Strouhal number St. The unstabilized scheme generates larger values of the first three indicators, whereas all stabilized schemes agree well with the existing data [5,6,22,32,42,56,66]. The predicted timevarying C d and C l underline the negligible deviation between φ = 0.1Δt and 1.0Δt in Fig. 9. The expenditure examined in Table 3 explains that φ = 0.25Δt consumes the least time again. Unlike the steady flow, adjusting φ does not deteriorate the accuracy. In Fig. 10, the vorticity contour using φ = 0.25Δt reflects that a repeating pattern of swirling vortices is caused by unsteady separation of the flow around the blunt body.    Table 4 analyzes the φ-sensitivity through them = 0.471 case. We see that the unstabilized scheme begets a failure whereas the stabilized schemes give nearly identical data. We choose φ = 1.0Δt for all mass ratios, given its performance. Eq. (33) still holds for FSI as the stability criterion because φ = 1.0Δt < Δt crit = 0.2828.

Vortex-induced vibration of a very light circular cylinder
The time history of aerodynamic parameters is plotted in Fig. 13 form = 0.408, at which our FSI method establishes the stable and smooth cylinder response. However, the enlarged view in Fig. 14a indicates the failure of the traditional strong staggered coupling (SSC) scheme [37]. By contrast, the present coupling scheme based upon standard fixed-point iterations agrees well with that computed by the nonlinear interface force correction (NFIC) approach [37]. Fig. 15 shows the x 1 -x 2 trajectory at variousm, illustrating that the VIV at low Re is a self-limiting process [76]. The cylinder takes on the nearly symmetrical trajectory shaping the classical Lissajous figure of "8". Vorticity fields atm = 0.393, 0.298 and 0.157 are displayed in Fig. 16 where the 2S vortex-shedding mode [69] is seen in the wake.
The unsteady periodic long-term oscillatory vibration of the tip is fairly depicted in Fig. 19a. Figure 19b demonstrates that the slightly longer time is required to reach the smaller characteristic amplitude in [28]. However, the under-estimated amplitude may be obtained even though denser spatial discretization or higher-order interpolation is used for the beam [15,20,40,41]. Three typical snapshots of vorticity and pressure fields are displayed in Fig. 20. It is seen that transient flow patterns and structural oscillations vary significantly in different phases.

Conclusions
This paper has reported the straightforward implementation of CS-FEM into two major areas of CFD. The stabilized second-order S-CBS scheme is proposed to solve incompressible NS equations. In the fluid equations GPs cooperates with equal SCs for each smoothed element integral, whereas CS-FEM works for the solid routinely. The structural equations are advanced in time by the Generalized-α method. The dynamic mesh is efficiently updated via MSA in combination with spring analogy method. Block-Gauss-Seidel procedure is adopted for the fluid-structure interplay within the ALE framework. The proposed methodologies do not only make trivial revision to available FE codes, but also exhibit outstanding performance in numerical tests. The main findings are summarized below -The SPGP technique is crucial to the second-order S-CBS scheme in incompressible flows and FSI. -The stabilization parameter has an impact on numerical accuracy and efficiency. In particular, Eq. (34) is recommended for the Eulerian flows while Eq. (33) for FSI. -The FSI solver never asks for accelerated fixed-point iterations even in the case of extremely low mass ratio.