Moving least‐squares aided finite element method: A powerful means to predict flow fields in the presence of a solid part

With the assistance of the moving least‐squares (MLS) interpolation functions, a two‐dimensional finite element code is developed to consider the effects of a stationary or moving solid body in a flow domain. At the same time, the mesh or grid is independent of the shape of the solid body. We achieve this goal in two steps. In the first step, we use MLS interpolants to enhance the pressure (P) and velocity (V) shape functions. By this means, we capture different discontinuities in a flow domain. In our previous publications, we have named this technique the PVMLS method (pressure and velocity shape functions enhanced by the MLS interpolants) and described it thoroughly. In the second step, we modify the PVMLS method (the M‐PVMLS method) to consider the effect of a solid part(s) in a flow domain. To evaluate the new method's performance, we compare the results of the M‐PVMLS method with a finite element code that uses boundary‐fitted meshes.


INTRODUCTION
Many physical and industrial problems involve a moving or stationary solid in a flow domain, for example, gear lubrication, [1][2][3] continuous stirred-tank reactor (CSTR) mixers, [4][5][6] static mixers, 7 and polymer mixing equipment and processing devices, such as single-screw extruders, 8,9 twin-screw extruders 10,11 and batch mixers. 12,13 popular way to evaluate the flow field parameters (the velocity components and the pressure values) in a solid-fluid flow problem is to generate a boundary-fitted mesh or grid in a Lagrangian framework for the flow domain and update the grid by the solid part movement.This strategy requires a mesh-regeneration step as the solid part moves.This paper uses the classical method to refer to this method.
Since mesh generation is time-consuming, using a background mesh for calculations is desirable.The background mesh comprises elements that cover the whole domain, including the fluid and solid part(s), independent of the shape and orientation of the undeformable part(s).For example, Fard et al. 14 adopted an existing mapping method 15 to evaluate the flow field in the cross-section of an internal mixer.In their research, the solid parts (the rotors of the mixer) are assumed as a fluid with very high viscosity compared to the first phase and behave like a solid body as a result.They extended their method to calculate the flow parameters in a 3D geometry and employed an extended finite element (XFEM) to evaluate them 16 more accurately.The XFEM enriches the Galerkin finite element method's regular or standard shape functions by adding a new term, enabling it to predict different discontinuities. 17There are various suggestions for enriched shape function 18 ; however, they all lead to an unsymmetrical stiffness matrix, which limits our choices for direct and iterative solvers.Notably, the XFEM introduces new unknowns or degrees of freedom to the finite element formulation.
The cut FEM [19][20][21][22] is another alternative to consider discontinuities, such as solid parts, in a flow field.This method uses a Nitsche formulation to calculate the integration of the weighted residual form, which is based on shape functions that are continuous in space but discontinuous in each time step.This method employs a penalty parameter to calculate or enforce the interface velocity (e.g., no slip at the solid and fluid interface).This method leads to an unsymmetric stiffness matrix, destroying the solution for some direct and iterative solvers.
The moving least-squares technique already has applications in solving discontinuous problems with meshless methods, 23,24 and the MLS interpolation functions are used to develop new shape functions. 25The intrinsic XFEM 26,27 uses the MLS function to create different shape functions for each node, depending on its position according to the interface.The developed shape functions do not have the Kronecker- properties (despite finite element standard shape functions), which makes it impossible to apply the Driklet boundary condition directly.Also, the support domain of each node, which includes the nodes of the neighboring and the next-neighboring elements, increases the bandwidth of the stiffness matrix and the calculation volume.
Our main target in this research is to introduce a new method to evaluate the velocity and pressure values in flow fields containing a constant or stationary solid body.This method, which is based on the FEM and employs the MLS interpolation functions, 25,27 requires a fixed background mesh independent of the shape and the orientation of the solid part.Also, it introduces no new unknowns (degree of freedom) into the finite element formulation and applies no-slip boundary conditions directly on the solid and fluid interface.It preserves the stiffness matrix's symmetry, giving us more direct or iterative solver choices.
We achieve this target in two steps.In the first step, we employ the moving least-squares (MLS) interpolation functions to enhance the pressure (P) and the velocity (V) shape functions.Therefore, it would be possible to consider the effect of the discontinuous parameters (i.e., the velocity and pressure values) in the finite element calculations.We have already described up to this point in 28,29 under the name of the PVMLS method (pressure and velocity shape function enhanced by the MLS interpolants).The PVMLS method calculates the pressure and velocity values in two-phase flow fields more accurately, especially in the vicinity of the interface, compared to the classical FEM.However, the PVMLS method is not a proper choice when the second phase is a stationary or moving solid body.Therefore, as the second step, we modify the PVMLS method to evaluate the flow field parameters in a domain containing both fluid and solid parts (the M-PVMLS method).
The content of this paper is as follows.After this section, the next part briefly discusses the governing equations of a two-phase flow problem and our previous solution strategies.In the third section, we suggest a method for applying the Dirichlet boundary condition on the solid-liquid interface and introduce the M-PVMLS method.Then, in the fourth and fifth parts, the performance of the M-PVMLS method in predicting the flow field parameters for Newtonian and generalized Newtonian fluids is investigated.Finally, we will conclude.

THE GOVERNING EQUATIONS AND SOLUTION STRATEGIES
Equations (1-a) and (1-b) mathematically express the creeping motion of an incompressible fluid and its continuity condition, respectively.The equations used in Einstein's summation notation over a similar index.The expanded formulation can be found in. 16 ij x j = F  i (1a) F I G U R E 1 An example of an active element and its division into sub-elements.Solid points q 1 q 2 q 3 q 4 show the quasi-nodes of a sub-element.[Colour figure can be viewed at wileyonlinelibrary.com] In equations (1-a) and (1-b), x i (i = 1, 2) are the coordinates of any arbitrary point (⃗ x) in a two-dimensional space.The components of the velocity vector (⃗ u) in a Cartesian coordinate are denoted by u i (i = 1, 2).The total stress tensor () comprises the scalar values of  ij (i, j = 1, 2) and is evaluated using equation (1-c).In equations (1-c) and (1-d),  ij and D ij are components of viscous stress () and deformation rate (D) tensors, respectively.The parameters P and  are hydrostatic pressure and Newtonian viscosity, correspondingly. ij is the unit tensor (I) component, such that  ij is one if i = j and otherwise 0 (Kronecker delta).The interfacial tension force ( ⃗ F  ) is shown by , which is thoroughly described in. 30In our problem, the surface tension force is zero, because the second phase is a solid body, and we assume a no-slip boundary condition at the solid-fluid interface.
Our finite element analysis requires a background mesh covering the whole domain, including solid and fluid portions, as shown in Figure 1.We employ a Galerkin finite element formulation to evaluate the weighted residual forms of the governing equations, as read in equations (2-a) to (2-e).

[ K jj
] The one-dimensional matrices  and  in equations (2-a) to (2-d) are bi-linear and bi-quadratic standard continuous finite element shape functions, respectively.The row and column vectors are shown by ⟨⟩ and {} symbols, respectively.The terms {u i } (i = 1, 2) and {P} are 9 × 1 and 4 × 1 matrices containing the velocity vector components and pressure values on each element's nodes, respectively.In these equations, the element domain and its boundary are expressed by Ω e and Γ e , correspondingly.⃗ n i is the unique base vector in x i direction and ⃗ n e is the normal vector to the element boundary.The calculations mentioned in equations (2-a) to (2-e) are performed for each element separately, and the integration is estimated using a third-order Gaussian quadrature method.Then, the outcomes are assembled to form the equation set of the mixed formulation, 31 by which one simultaneously calculates velocity components and pressure values.The equation set is written in the matrix form, in which the coefficient matrix (the stiffness matrix) is symmetric. 16,31e used the  and  interpolation functions for the whole domain; however, they need some enhancement on active elements or the elements containing two phases (i.e., solid and fluid).Our previous publications 28,29 describe enhancing the active elements' pressure and velocity shape functions; however, one can find a brief description here.
Figure 1 shows a typical active element (ABCD), which is divided into several minor elements, namely sub-elements (e.g., q 1 q 2 q 3 q 4 ).Sub-elements have three properties.They cover the active element area and do not overlap each other.Also, they are occupied by one phase only.
As we have shown in Figure 1 by solid points, a sub-element introduces some new nodes (quasi-node) to the problem (e.g., q 1 q 2 q 3 q 4 ).To avoid any new degree of freedom in our finite element formulation, we employ the moving least-squares (MLS) interpolation function to correlate the value of a parameter (e.g., velocity or pressure) on a quasi-node to the nodes in its vicinity.We have described the MLS calculations in Appendix A.
To provide an example, we have clarified the pressure nodes in the neighborhood of quasi-nodes q i ( i = 1, 2, 3, 4) in Figure 1 by hollow circles.Equation (3-a) reveals that the coefficient matrix [ MLS p ] 4×s , calculated using the MLS method, correlates the pressure values ) on the quasi-nodes (q i , i = 1, 2, 3, 4) of a typical sub-element with their surrounding nodes ) .
To evaluate the pressure inside each active element (e.g., x * in Figure 1), one can use continuous shape functions , as mentioned in equation (3-b).By considering equations (3-a) and (3-b), we can calculate the enhanced shape functions  * , as stated in equation (3-c).We employ a similar procedure to evaluate the enhanced shape function of the velocity ( * ), which is not mentioned here for conciseness.To our knowledge,  * and  * are continuous shape functions in each sub-element. [ As depicted in equations (4-a) to (4-e), the enhanced shape functions  * and  * take the places of  and  (just for active elements) in equations (2-a) to (2-e), respectively.Therefore, it is possible to consider discontinuities of pressure and velocity components in two-phase flow problems, as described under the PVMLS method in. 28The terms { u * i } and { P * } are the velocity vector components and pressure values on the nodes in the vicinity of each active element, respectively.
Since the enhanced shape functions  * and  * are not polynomials anymore, integrating equations (4-c) to (4-e) requires special precautions.Using the second-order Gaussian quadrature integration, we observed results with checkered plate marks (especially for velocity contours), which vanish by employing the Gaussian quadrature integration with third-order and higher.Replacement of the Gaussian quadrature integration with other integration techniques, such as kernel gradient smoothing integration, 32 could be an alternative to reduce calculation time.
Although the PVMLS method mainly deals with two-phase flow systems containing two different fluids, it is applicable for solid-fluid problems if one phase is assumed to have a very high viscosity value and, consequently, undeformable.In this research, we modify the PVMLS method (M-PVMLS) to apply the Dirichlet boundary condition directly to the solid-fluid interface to avoid unnecessary finite element calculations for the elements and sub-elements occupied by a solid phase.

APPLYING THE DIRICHLET BOUNDARY CONDITIONS ON THE SOLID-FLUID INTERFACE
In the following, we will apply no-slip conditions on the solid-liquid boundaries.When the no-slip condition, as a Dirichlet boundary condition, is imposed on the weighted residual form of the governing equations, it specifies the velocity value a node alongside the boundary must take. 31It means that the velocity of each boundary node is zero relative to the solid boundary.The no-slip boundary condition readily applies to the classic FEM, 33 as their shape functions have the Kronecker- properties.Since such a property cannot always be provided for the moving least-squares interpolation functions, applying Dirichlet boundary condition is a challenge for the meshless methods 27 or the intrinsic extended finite element methods. 25he M-PVMLS method directly enforces the no-slip boundary condition; however, the difficulty of imposing the Dirichlet boundary condition in the M-PVMLS method arises because the nodes on the boundaries (fluid-solid interfaces), as shown by triangles in (Figure 1), do not fit in the background meshing system.In other words, they are quasi-nodes and not the main nodes.
We describe the incorporation of the Dirichlet boundary condition for a typical term of . One can apply a similar procedure for the remaining terms in equation (4-a).
As with the PVMLS method, it turns out that the vector  * (the enhanced shape function) correlates the velocity value of a quasi-node to the velocity of the surrounding nodes (say n vel nodes).However, to apply the boundary condition on the quasi-nodes, the mentioned correlation is done for n vel surrounding nodes plus three quasi-nodes on the interface, as shown schematically by triangles in (Figure 1).In this situation, the term can be written in the following form: In equation ( 5), the terms u i,m l (l = 1, 2, 3) are the velocity components of the quasi-nodes on the fluid-solid boundary.Although we calculate the matrix according to equation (4-c), we cannot assemble it to form the equation set of the mixed finite element formulation; while the nodes m l (m l = n vel + l, l = 1, 2, 3) do not belong to the background meshing system.
Therefore, we apply the Dirichlet boundary condition on these quasi-nodes to remove the related data from the local stiffness matrix . In other words, since the velocity values of the quasi-nodes on the interface are known, the boundary condition can be directly applied to the weighted residual form, as stated in equation (6-a).To have a symmetrical coefficient matrix of unknowns (u i,1 , u i,2 … and u i,n vel ) we decompose into a symmetric and a non-symmetric matrix, as read in equation (6-b).
The three last lines of the matrixes in the first term of equation (6-b) lead to an apparent equality (u i,m l = u i,m l ) and do not take part in assembling the mixed finite element method equation set.Therefore, the contribution of the term on the stiffness matrix, and the right-hand side of the equation set is re-written as equations (7-a) to (7-c): The matrix [A] can be used in assembling the mixed finite element method equation set and the vector {B} is considered and treated a source term on the right-hand side of the equation set.Notably, the matrix [A] is symmetric, and applying the boundary condition will not ruin the symmetric nature of the total stiffness matrix.

EVALUATING M-PVMLS METHOD IN DIFFERENT GEOMETRIES
As stated, the M-PVMLS method employs a fixed background mesh covering the whole domain, including fluid and solid parts.Similar to the methods that comprise a background mesh, 34,35 we compare the results of the M-PVMLS method with an analysis using a boundary-fitted mesh. 36,37The selected geometries and the flow conditions for the comparison are as follows: • We consider a parallel channel as a flow domain.We arbitrarily choose the channel's depth and length to be 2h and 6h (h = 0.5).
• The viscosity of the fluid is unity.• The solid part is supposed to have circular or rectangular shapes.
• The solid part can be fixed or rotating • The driving force for the fluid to flow is the movement of the channel's wall (similar to the simple shear flow), the pressure difference at the channel's two ends, or the solid part's rotation.
First, we consider a domain similar to the simple shear flow with a circular obstacle in the middle of the channel.As we stated before, for the M-PVMLS method, the mesh is generated for the whole domain, while the classic method excludes the solid part.Such a meshing system is an advantage of the M-PVMLS method since changing the obstacle shape and position requires no new mesh.(Figure 2A,B) show the background mesh (suitable for the M-PVMLS method) and the classic method, respectively.
The flow field parameters, including the velocity in the channel length direction, the transverse flow direction, and the pressure values, are compared in (Figure 2C-H).Both methods reveal a similar deviation from the fully developed flow (according to the velocity components values and the streamlines) in the vicinity of the obstacle and a high-pressure gradient in the same area.The similarity between the results and the streamlines, as shown in (Figure 2I,J), demonstrates the capability of the M-PVMLS method in evaluating the flow field parameters.
We have similarly studied the pressure flow, as shown in (Figure 3A-J), and the rotating solid part in (Figure 4A-J).
In (Figures 5-7), we have replaced the obstacle with a square-shaped block with a dimension of h × h.As (Figures 5A, 6A, and 7A) show, the mesh of the M-PVMLS method is the same as the cases with the circular obstacle, while the classic method requires a new mesh.
A close survey of (Figures 2-7) demonstrates that the results of the M-PVMLS method comply with the outcomes of the boundary-fitted mesh, independent of the solid part's shape and the driving force of the flow field.We have used the comparison with the classic FEM to validate the M-PVMLS method.
Notably, this section's flow problems are typical examples of our numerous numerical experiments.However, for the sake of brevity, they are not all listed here.

EVALUATING THE M-PVMLS METHOD FOR GENERALIZED NEWTONIAN FLUIDS
This section uses the M-PVMLS method to calculate the pressure and velocity values in a flow field comprising a shear-thining fluid.For this purpose, we replace the constant viscosity () in equation (1-d) with a shear-rate dependant viscosity ().In our study, we have used a power-law model, as stated in equations (8-a) to (8-d)

APPENDIX A
The moving least-squares method correlates the value of a parameter on a point to the scattered points in its neighborhood.
A comprehensive literature review on this subject in 38 is available.We have also used this technique to correlate the value of a parameter (e.g., pressure) on the sub-element's quasi-nodes to their surrounding nodes (i.e., background mesh's nodes in their vicinity).
Figure 1A shows a typical sub-element (q 1 q 2 q 3 q 4 ) in the ABCD element, and solid circles highlight all surrounding nodes.We describe the evaluation of the MLS coefficients, which correlate the quantity of the pressure on a typical quasi-node (e.g., q 1 ) to its surrounding nodes in six steps.In our calculation, we neglect the surrounding nodes, which are not in the same phase as the quasi-nodes, 39 and consider the ones (P * j , j = 1, 2, … , s) distinguished by hollow squares in Figure 1A. ) ) and the interface radius (R int ), as stated in equations (a) and (b).We have already defined the interface radius in, 30 which is proportional to the element size.Notably, in the new coordinate, the spatial position of the quasi-node is (0, 0).
F I G U R E 1A An example of a sub-element and its surrounding nodes.Hollow points q 1 q 2 q 3 q 4 show the quasi-nodes of a sub-element.
Step 2: Select a suitable polynomial basis vector.A linear basis vector is ideal for pressure, as is shown in equation (c).We recommend a higher-order polynomial for the velocity basis vector. 28b Step3: Form the following matrix: Step 4: Select a suitable weight function, which decreases as the distance from the quasi-node increases, as shown in equation (e). w In equation (e), the parameter c is a real and positive constant, which is proportional to the element size.
Step 5: Form the following diagonal matrix.

E 2
Flow field parameters (the velocity components, the pressure values, and the streamlines) in a geometry similar to the simple shear flow with a circular obstacle.(a) and (b): generated mesh for the M_PVMLS and the classic methods, respectively; (c) and (d): The contour of the velocity in the channel direction for the M_PVMLS and the classic methods, respectively; (e) and (f): The contour of the velocity in the transverse channel direction for the M_PVMLS and the classic methods, respectively; (g) and (h): The pressure contour for the M_PVMLS and the classic methods, respectively; (i) and (j): The streamlines for the M_PVMLS and the classic methods, respectively.[Colour figure can be viewed at wileyonlinelibrary.com]

3
Flow field parameters (the velocity components, the pressure values, and the streamlines) in a geometry similar to pressure-flow with a circular obstacle.(a) and (b): generated mesh for the M_PVMLS and the classic methods, respectively; (c) and (d): The contour of the velocity in the channel direction for the M_PVMLS and the classic methods, respectively; (e) and (f): The contour of the velocity in the transverse channel direction for the M_PVMLS and the classic methods, respectively; (g) and (h): The pressure contour for the M_PVMLS and the classic methods, respectively; (i) and (j): The streamlines for the M_PVMLS and the classic methods, respectively.[Colour figure can be viewed at wileyonlinelibrary.com]

4
Flow field parameters (the velocity components, the pressure values, and the streamlines) in a channel with a rotating circle in the middle.(a) and (b): generated mesh for the M_PVMLS and the classic methods, respectively; (c) and (d): The contour of the velocity in the channel direction for the M_PVMLS and the classic methods, respectively; (e) and (f): The contour of the velocity in the transverse channel direction for the M_PVMLS and the classic methods, respectively; (g) and (h): The pressure contour for the M_PVMLS and the classic methods, respectively; (i) and (j): The streamlines for the M_PVMLS and the classic methods, respectively.[Colour figure can be viewed at wileyonlinelibrary.com]

5
Flow field parameters (the velocity components, the pressure values, and the streamlines) in a geometry similar to the simple shear flow with a square obstacle.(a) and (b): generated mesh for the M_PVMLS and the classic methods, respectively; (c) and (d): The contour of the velocity in the channel direction for the M_PVMLS and the classic methods, respectively; (e) and (f): The contour of the velocity in the transverse channel direction for the M_PVMLS and the classic methods, respectively; (g) and (h): The pressure contour for the M_PVMLS and the classic methods, respectively; (i) and (j): The streamlines for the M_PVMLS and the classic methods, respectively.[Colour figure can be viewed at wileyonlinelibrary.com]

E 6
Flow field parameters (the velocity components, the pressure values, and the streamlines) in a geometry similar to pressure-flow with a square obstacle.(a) and (b): generated mesh for the M_PVMLS and the classic methods, respectively; (c) and (d): The contour of the velocity in the channel direction for the M_PVMLS and the classic methods, respectively; (e) and (f): The contour of the velocity in the transverse channel direction for the M_PVMLS and the classic methods, respectively; (g) and (h): The pressure contour for the M_PVMLS and the classic methods, respectively; (i) and (j): The streamlines for the M_PVMLS and the classic methods, respectively.[Colour figure can be viewed at wileyonlinelibrary.com]

7
Flow field parameters (the velocity components, the pressure values, and the streamlines) in a channel with a rotating square in the middle.(a) and (b): generated mesh for the M_PVMLS and the classic methods, respectively; (c) and (d): The contour of the velocity in the channel direction for the M_PVMLS and the classic methods, respectively; (e) and (f): The contour of the velocity in the transverse channel direction for the M_PVMLS and the classic methods, respectively; (g) and (h): The pressure contour for the M_PVMLS and the classic methods, respectively; (i) and (j): The streamlines for the M_PVMLS and the classic methods, respectively.[Colour figure can be viewed at wileyonlinelibrary.com]

8
Flow field parameters (the velocity components, the pressure values, and the streamlines) for a fluid obeying the power-law model in a geometry similar to pressure-flow with a circular obstacle.(a) and (b): generated mesh for the M_PVMLS and the classic methods, respectively; (c) and (d): The contour of the velocity in the channel direction for the M_PVMLS and the classic methods, respectively; (e) and (f): The contour of the velocity in the transverse channel direction for the M_PVMLS and the classic methods, respectively; (g) and (h): The pressure contour for the M_PVMLS and the classic methods, respectively; (i) and (j): The streamlines for the M_PVMLS and the classic methods, respectively.[Colour figure can be viewed at wileyonlinelibrary.com]

9
Flow field parameters (the velocity components, the pressure values, and the streamlines) for a fluid obeying the power-law model in a geometry similar to the pressure flow with a rectangular obstacle.(a) and (b): generated mesh for the M_PVMLS and the classic methods, respectively; (c) and (d): The contour of the velocity in the channel direction for the M_PVMLS and the classic methods, respectively; (e) and (f): The contour of the velocity in the transverse channel direction for the M_PVMLS and the classic methods, respectively; (g) and (h): The pressure contour for the M_PVMLS and the classic methods, respectively; (i) and (j): The streamlines for the M_PVMLS and the classic methods, respectively.[Colour figure can be viewed at wileyonlinelibrary.com]