Moving Least-Squares Aided Finite Element method (MLS-FEM): A Powerful Means to Evaluate Distributive Mixing

1. Introduction

In view of the growing demand for polymeric materials with tailored properties, efforts are being made to develop specialized materials that meet different requirements. Among the various approaches, compounding and blending have proven to be the most cost-effective and viable methods for producing tailored polymeric materials [1]. The quality of the final products obtained through these processes is significantly influenced by the effectiveness of the mixing stage [2,3].

During the polymer mixing process, at least two distinct phases are incorporated through the bulk diffusion mechanism, wherein the second phase distributes and disperses within the main polymeric phase [4]. Explicit models can be employed to investigate the development of morphology or microstructural changes that transpire during the mixing process [5]. Alternatively, implicit parameters offer another approach to assessing the extent of mixing. These implicit parameters encompass various factors such as residence time [6,7], total strain [6,7,8], interfacial area [6,9] in passive mixing [5,10], the concentration distribution [11,12], the ratio of the deformation rate to the vorticity magnitude [6,11,12,13,14] and other parameters [6,14,15] such as frame invariant. These parameters enable the evaluation of mixing quality based on the flow field of the single-phase system.

While the single-phase assumption simplifies the numerical assessment of mixing extent, it presents a significant challenge due to moving parts in mixing devices. These moving parts alter the flow domain throughout the mixing process, introducing complexities in the analysis. An approach for such a problem is to define a boundary-fitted mesh and to evaluate the flow field parameters (e.g., the components of the velocity vector and pressure values) as the positions of the internal parts change. This approach has initiated extensive application in evaluating implicit mixing parameters in various mixing devices, including twin-screw extruders [9,16,17,18,19] and internal mixers [12,20].

Given the intricate geometries of polymer mixers, the mesh generation process in numerical solutions can be particularly time-consuming. To overcome this challenge, researchers have made significant efforts to develop techniques that simplify the mesh generation step in such complex scenarios. For instance, using mesh generation based on boundary-fitted curvilinear coordinates [21,22] facilitates grid generation, albeit necessitating adjustments when internal parts' positions change. The moving mesh technique [23], initially developed for the integral form of governing equations, employs a fixed mesh. Though, in cases where grids are highly deformed or distorted, mesh regeneration becomes necessary. The sliding mesh approach [24] uses two complementary and fixed grids that slide against each other. Nevertheless, this method requires data interpolation between two separate grids during each time step, rendering it less practical. Additionally, it is not applicable in situations where the moving components of the geometry overlap, such as in intermeshing twin-screw extruders.

Lawal et al. [25] considered a background mesh covering fluid and solid parts (the rotors in a polymer mixing device). They enforced the velocity of the rotor in the finite element formulation via a penalty term to avoid any mesh regeneration during the mixing process. The fictitious domain method [26] also uses a static or background mesh covering fluid and solid parts. Within this method, the solid parts are kinematic constraints with known velocity values. Also, some controlling points describe the solid part's shape and velocity. Bertrand et al. [27] coupled this method with a mesh refinement technique to increase the accuracy of the flow field calculations. Fard et al. [5] use the fictitious domain method to evaluate the velocity field in an internal mixer. Then Fard and Anderson [28] employed the extended finite element method (XFEM) to improve the accuracy of the flow field calculations. Mesh refinement and the XFEM introduce new unknowns to the finite element formulation and increase the computation time.

Other methods, such as the cutFEM [29] and cut cell method [30,31,32], employ a background mesh in a flow field containing a solid. Nevertheless, these method has not yet been used to evaluate a mixer's flow field.

Our goal in this research is to implicitly evaluate the distributive mixing in an internal mixer in two-dimension. We assume a filled chamber with two revolving rotors (solid parts). We employ the finite element method (FEM) [33,34] to evaluate the flow field parameters. Also, to increase the accuracy and avoid remeshing during the rotors' rotations, we have used the moving least-squares (MLS) interpolation function, as described in [35,36], to enhance the classic finite element shape functions. This method neither introduces new unknowns to the finite element formulation nor requires a local mesh refinement.

The paper encompasses the following content. The second section delves into the governing equations of a two-phase flow problem while providing a concise overview of our solution strategy. The internal mixer's geometry is presented in Section 3, followed by Section 4, where we employ the MLS-FEM to compute flow field parameters in the internal mixer's cross-section. We compare these results with those obtained from boundary-fitted or classic FEM. Moving forward, Section 5 evaluates the mean total strain function as an implicit parameter for distributive mixing. Both MLS-aided and boundary-fitted FEM approaches are employed, showcasing their comparable outcomes. Ultimately, Section 6 encapsulates our concluding remarks.

2. The governing equations

The creeping motion of an incompressible fluid and its continuity condition can be mathematically represented by Equations (1-a) and (1-b), respectively. These equations employ Einstein's summation notation, utilizing a common index.

$$\frac{\partial {\sigma }_{ij}}{\partial x_j}=F_{{\sigma }_i}$$

(1-a)

$$\frac{\partial u_i}{\partial x_i}=0$$

(1-b)

$${\sigma }_{ij}={\tau }_{ij}+p{\delta }_{ij}$$

(1-c)

$${\tau }_{ij}=-2\mu D_{ij}$$

(1-d)

$$D_{ij}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)$$

(1-e)

Equations (1-a) and (1-b) introduce the coordinates $x_i$ ($i=\mathrm{1,2}$), which represents an arbitrary point ($\overrightarrow{\mathit{x}}$) in a two-dimensional space. The velocity vector ($\overrightarrow{\mathit{u}}$), expressed in Cartesian coordinates, is represented by the components $u_i$ ($i=\mathrm{1,2}$). The total stress tensor ($\stackrel{̿}{\mathit{\sigma }}$) comprises the scalar values ${\sigma }_{ij}$ ($i,j=\mathrm{1,2}$) and is evaluated using equation (1-c). Equations (1-c) and (1-d) involve ${\tau }_{ij}$ and $D_{ij}$, which denote the components of the viscous stress ($\stackrel{̿}{\mathit{\tau }}$) and deformation rate ($\stackrel{̿}{\mathit{D}}$) tensors, respectively. The parameters $P$ and $\mu$ represent the hydrostatic pressure and Newtonian viscosity, respectively. The unit tensor ($I$), denoted by ${\delta }_{ij}$, has a value of 1 when $i=j$, and 0 otherwise (Kronecker delta). The interfacial tension force, represented by ${\overrightarrow{\mathit{F}}}_{\sigma }$, is described in detail in [35] through $F_{{\sigma }_i}$ ($i=\mathrm{1,2}$). In our specific problem, the surface tension force is disregarded as the second phase is a solid body, and we assume a no-slip boundary condition at the solid-fluid interface.

In our analysis, we utilize the Galerkin finite element formulation to assess the weighted residual forms of the governing equations, as outlined in equations (2-a) and (2-b).

$$\left[{\int }_{{\mathsf{\Omega }}^e}^{}\mu \left(\frac{\partial \left\{\mathit{\psi }\right\}}{\partial x_j}\frac{\partial 〈\mathit{\psi }〉}{\partial x_j}\right)d\mathit{x}\right]\left\{{\mathit{u}}_i\right\}+\left[{\int }_{{\mathsf{\Omega }}^e}^{}\mu \left(\frac{\partial \left\{\mathit{\psi }\right\}}{\partial x_j}\frac{\partial 〈\mathit{\psi }〉}{\partial x_i}\right)d\mathit{x}\right]\left\{{\mathit{u}}_j\right\}-\left[{\int }_{{\mathsf{\Omega }}^e}^{}\frac{\partial \left\{\mathit{\psi }\right\}}{\partial x_i}〈\mathit{\phi }〉d\mathit{x}\right]\left\{\mathit{P}\right\}=\left\{{\oint }_{{\mathsf{\Gamma }}^e}^{}\left\{\mathit{\psi }\right\}\left(\stackrel{̿}{\mathit{\sigma }}.{\overrightarrow{\mathit{n}}}_{\mathit{e}}\right).{\overrightarrow{\mathit{n}}}_{\mathit{i}}d\mathit{s}\right\}$$

(2-a)

$$\left[{\int }_{{\mathsf{\Omega }}^e}^{}\left\{\mathit{\phi }\right\}\frac{\partial 〈\mathit{\psi }〉}{\partial x_i}d\mathit{x}\right]\left\{{\mathit{u}}_i\right\}=0$$

(2-b)

The equations (2-a) and (2-b) comprise one-dimensional matrices $\mathit{\phi }$ and $\mathit{\psi }$, which represent bi-linear and bi-quadratic interpolation functions for pressure (P) and velocity (${\mathit{u}}_i$,$i=\mathrm{1,2}$), respectively. The symbols

and

denote the row and column vectors, correspondingly. We use ${\mathsf{\Omega }}^e$ and ${\mathsf{\Gamma }}^e$ to represent the element domain and its boundary, respectively. The vector ${\overrightarrow{\mathit{n}}}_{\mathit{i}}$ denotes the unique base vector in the $x_i$ direction, while ${\overrightarrow{\mathit{n}}}_{\mathit{e}}$ symbolizes the normal vector to the element boundary.

In the classic FEM, the elements comply with the boundary of the fluid domain. Therefore, we name it boundary-fitted FEM, as shown in Figure 1a. However, to conduct our finite element analysis, as shown in Figure 1b, we utilize a background mesh encompassing the entire domain, including solid and fluid regions. Therefore, this strategy, known as the moving least-squares aided finite element method (MLS-FEM), eliminates the requirement for remeshing, even when the solid parts are displaced or moved.

By implementing the MLS-FEM, we effectively exclude the solid part by dividing active elements along the solid/fluid boundary. In this context, active elements refer to elements that contain both solid and fluid components. Figure 1b illustrates a representative active element, denoted as $A_0{A}_1{A}_2{A}_3$, and demonstrates its splitting process. Element $A_0{A}_1{A}_2{A}_3$ consists of two sections, one occupied by fluid. We must divide the fluid-filled portion into rectangular elements to calculate the finite element. As illustrated in Figure 1b, new elements are created when the splitting process occurs, referred to as sub-elements. These sub-elements consist of nodes, known as quasi-nodes, not all of which are part of the background mesh.

We establish a relationship between the parameter's value (such as the pressure values) at each quasi-node and the surrounding nodal values using the moving least-squares technique to eliminate any extra unknown variables in the finite element formulation. The row-wise vector $〈{\mathit{M}\mathit{L}\mathit{S}}_i〉$ in equation (3-a) reads this correlation. It is important to note that our calculations do not consider nodes within the solid region. Equations (3-b) to (3-f) outline the procedure for computing the correlating vector ($〈{\mathit{M}\mathit{L}\mathit{S}}_i〉$).

$$P_i=〈{\mathit{M}\mathit{L}\mathit{S}}_i〉\left\{P_j\right\}, i=q_1,q_2,q_3,q_4, j=A_1, A_2, \dots , A_s$$

(3-a)

$${\overrightarrow{\mathit{\beta }}}_{i,j}=\left({\overrightarrow{\mathit{x}}}_j-{\overrightarrow{\mathit{x}}}_i\right)/c, {\overrightarrow{\mathit{\beta }}}_{i,j}=\left({\beta }_{{i,j}_1},{\beta }_{i,j_2}\right)$$

(3-b)

$$〈{\mathit{h}}_i\left({\overrightarrow{\mathit{\beta }}}_j\right)〉=〈1,{\beta }_{j_1},{\beta }_{j_2},{\beta }_{j_1}{\beta }_{j_2}〉$$

(3-c)

$$\left[{\mathit{H}}_i\right]=\left[\begin{array}{c}〈{\mathit{h}}_i\left({\overrightarrow{\mathit{\beta }}}_1\right)〉\\ \begin{array}{c}〈{\mathit{h}}_i\left({\overrightarrow{\mathit{\beta }}}_2\right)〉\\ ⋮\end{array}\\ 〈{\mathit{h}}_i\left({\overrightarrow{\mathit{\beta }}}_s\right)〉\end{array}\right]$$

(3-d)

$$w_i\left({\overrightarrow{\mathit{\beta }}}_{i,j}\right)=exp\left(-\sqrt{{\beta }_{{i,j}_1}^2+{\beta }_{{i,j}_2}^2}/c\right)$$

(3-e)

$$〈{\mathit{M}\mathit{L}\mathit{S}}_i〉=〈\mathrm{1,0},\mathrm{0,0}〉{\left({\left[{\mathit{H}}_i\right]}^{†}\left[{\mathit{W}}_i\right]\left[{\mathit{H}}_i\right]\right)}^{-1}{\left[{\mathit{H}}_i\right]}^{†}\left[{\mathit{W}}_i\right]$$

(3-f)

Equation (3-b) normalizes the coordinates of the neighboring nodes with respect to each quasi-node. The scalar value $c$ is in the same order as the element size. Notably, after normalization, the coordinates of the quasi-nodes become ${\overrightarrow{\mathit{\beta }}}_{i,i}=\left(0,0\right)$. Following this, we construct the coefficient vector $〈{\mathit{h}}_i\left({\overrightarrow{\mathit{\beta }}}_j\right)〉$ and the coefficient matrix $\left[{\mathit{H}}_i\right]$ using equations (3-c) and (3-d) correspondingly. Next, utilizing the diagonal weighting matrix $\left[{\mathit{W}}_i\right]$, with its components introduced in equation (3-e), we proceed to calculate the correlating vector $〈{\mathit{M}\mathit{L}\mathit{S}}_i〉$, as outlined in equation (3-f).

Utilizing equation (3-a), we can calculate the pressure value at the quasi-nodes of each sub-element, as described in equation (4-a). Using the classic shape function, equation (4-b) establishes a correlation between the pressure value within a quasi-element (e.g., $x^{\mathrm{*}}$ in Figure 1b) and the values at the corresponding quasi-nodes. By combining equations (4-a) and (4-b), we can introduce the enhanced shape function (${\mathit{\phi }}^{\mathrm{*}}$) as expressed in equations (4-c). This enhanced shape function allows a more accurate representation of the pressure distribution within the quasi-element.

$$\left\{\begin{array}{c}{P}_{q_1}\\ \begin{array}{c}{P}_{q_2}\\ P_{q_3}\end{array}\\ P_{q_4}\end{array}\right\}=\left[\begin{array}{c}〈{\mathit{M}\mathit{L}\mathit{S}}_{q_1}〉\\ \begin{array}{c}〈{\mathit{M}\mathit{L}\mathit{S}}_{q_2}〉\\ 〈{\mathit{M}\mathit{L}\mathit{S}}_{q_3}〉\end{array}\\ 〈{\mathit{M}\mathit{L}\mathit{S}}_{q_4}〉\end{array}\right]\left\{\begin{array}{c}{P}_{A_1}\\ \begin{array}{c}{P}_{A_2}\\ ⋮\end{array}\\ P_{A_s}\end{array}\right\}$$

(4-a)

$$P_{x^{\mathrm{*}}}=〈\mathit{\phi }〉\left\{\begin{array}{c}{P}_{q_1}\\ \begin{array}{c}{P}_{q_2}\\ P_{q_3}\end{array}\\ P_{q_4}\end{array}\right\}$$

(4-b)

$$P_{x^{\mathrm{*}}}=〈\mathit{\phi }〉\left[\begin{array}{c}〈{\mathit{M}\mathit{L}\mathit{S}}_{q_1}〉\\ \begin{array}{c}〈{\mathit{M}\mathit{L}\mathit{S}}_{q_2}〉\\ 〈{\mathit{M}\mathit{L}\mathit{S}}_{q_3}〉\end{array}\\ 〈{\mathit{M}\mathit{L}\mathit{S}}_{q_4}〉\end{array}\right]\left\{\begin{array}{c}{P}_{A_1}\\ \begin{array}{c}{P}_{A_2}\\ ⋮\end{array}\\ P_{A_s}\end{array}\right\}=〈{\mathit{\phi }}^{\mathrm{*}}〉\left\{\begin{array}{c}{P}_{A_1}\\ \begin{array}{c}{P}_{A_2}\\ ⋮\end{array}\\ P_{A_s}\end{array}\right\}$$

(4-c)

In equations (2-a) and (2-b), we substitute the classic bilinear shape function ($\mathit{\phi }$) with an enhanced shape function (${\mathit{\phi }}^{\mathrm{*}}$) for active elements. This implies that we utilize the classic FEM for non-active elements or the elements occupied by the fluid.

The calculation of the enhanced bi-quadratic shape function (${\mathit{\psi }}^{\mathrm{*}}$) for velocity values follows a similar approach. However, it incorporates two additional considerations to achieve more accurate results. The first point is the coefficient vector (as described in equation 3-c) replacement with a higher-order polynomial. Secondly, we correlate the velocity of each quasi-node to not only the surrounding nodes but also the quasi-nodes lying on the solid boundary (e.g., $q_0,q_1,$ and $q_2$ in Figure 1b). Considering that the quasi-nodes situated on the solid boundary have known velocities, we treat the contribution of these nodes as a source term in the finite element formulation.

3. The geometry of the mixer

Our study analyzes flow patterns and mixing within a two-dimensional internal mixer, as illustrated in Figure 2. We name this mixer's left and right parts as the left and right sub-chamber, respectively, as written in Figure 2. It's worth noting that similar analyses are possible for mixers of varying sizes, geometries, and configurations.

4. Validating the MLS-FEM

The validation process for methods employing a background mesh involves comparing the flow field parameters, such as velocity components and pressure values, with the boundary-fitted mesh FEM [36,37]. Here, we have used a similar procedure. On the other hand, we calculate flow field parameters (the velocity component and pressure values) using the MLS-FEM and boundary-fitted mesh FEM and compare the results. If the results of the two methods comply, it confirms that the MLS-FEM approach is a viable alternative that can yield accurate results without mesh regeneration.

In this section, we have provided two cases for comparison: The first comprises counter-rotating symmetric rotors, and the second includes co-rotating asymmetric rotors. We have conducted several comparisons for various rotor arrangments and rotation directions; however, these data are not presented here for conciseness.

Figures 3a and 4a illustrate the generated mesh for the MLS-FEM (Moving Least Squares Finite Element Method). This mesh consists of a background mesh that covers both the fluid-occupied and solid-occupied regions. Figures 3b and 4b show the grid for the boundary-fitted mesh FEM, which complies with the solid boundary and requires remeshing as the rotors rotate. It is worth noting that in the MLS-FEM, there is no need to generate a new mesh when the positions of the rotors change.

Figure 3c–h and Figure 4c–h compare the flow field parameters achieved by the MLS-FEM and the boundary-fitted mesh FEM. The similarity between the results of the two methods demonstrates the MLS-FEM’s capability in evaluating the flow field parameters.

5. Distributive mixing:

This section uses MLS-FEM to qualitatively and quantitatively evaluate distributive mixing in the two-dimensional internal mixer and compares it with those of boundary-fitted mesh FEM. For this purpose, we consider the following assumptions:

We arbitrarily choose an asymmetrical arrangement and a counter-rotating operating condition.

We calculate the flow field parameters every $∆\theta =5$ degrees of rotation. Therefore, it is necessary to generate mesh after each 5-degree rotation for the boundary-fitted mesh FEM; nevertheless, the M-PVMLS method uses a fixed background mesh.

We use a Newtonian fluid assumption in our calculation; however, extending to generalize Newtonian with shear thinning behavior is not a big challenge.

The rotor speed will not affect the results since our finite element calculations use creeping flow assumption. So, we choose the rotor angular velocity to be $\omega =10Rounds/min$. By this assumption, every 5-degree rotation will take $∆t=1/12$ seconds.

We introduce four thousand (4000) inert particles to the system and trace them during the first five revolutions of the rotors.

Also, two colors are chosen for each sub-chambers particle to observe the particle exchange between sub-chambers better, as shown in Figure 5a,g.

For particle tracking, we employ equation (5):

$${\overrightarrow{\mathit{x}}}_{p_i}^{\mathrm{*}\mathrm{*}}={\overrightarrow{\mathit{x}}}_{p_i}^{\mathrm{*}}+{\overrightarrow{\mathit{u}}}_{{\mathit{p}}_{\mathit{i}}}\times \delta t,i=\mathrm{1,2},\dots ,n_{particle}$$

(5)

In equation (5), ${\overrightarrow{\mathit{x}}}_{p_i}^{\mathrm{*}}$ and ${\overrightarrow{\mathit{x}}}_{p_i}^{\mathrm{*}\mathrm{*}}$ are the current and the new spatial coordinates of the $i^{th}$ particle. The parameter ${\overrightarrow{\mathit{u}}}_{{\mathit{p}}_{\mathit{i}}}$ is the velocity vector of the particle at ${\mathit{x}}_{p_i}^{\mathrm{*}}$, calculated by interpolating among the velocity nodal values inside each element. We perform this calculation for all of the randomly distributed particles ($n_{particle}$) in the mixing chamber. The parameter $\delta t$ is the time increment in the particle tracking procedure.

We choose the time increment $\delta t$ so that the particle movement in each step is less than the element length. Disregarding the time constraint can lead the particles beyond the flow domain, interpreted as particle loss [38]. Therefore the mentioned $∆t$ is divided into smaller $\delta t$ to maintain the time increment constraint; however, particle loss may be observed, especially for the flow fields with moving parts.

Figure 5a–f show the distribution of the particles from the startup to the fifth revolution, using The MLS-FEM. A comparison of these results with those obtained using the boundary-fitted mesh FEM, as shown in Figure 5g–l, reveals qualitative similarity in the mixing patterns between the two methods.

We compare two different mixing parameters to provide a quantitative comparison between the MLS-FEM and the boundary-fitted mesh FEM results. The first parameter is the particle exchange between two sub-chambers. Secondly, the mean average of the total strain ($\overline{\gamma }$), as a criterion for distributive mixing, is also calculated based on equations (6-a) and (6-b).

$$\overline{\gamma }=mean\left({\gamma }_{p_i}\right)$$

(6-a)

$${\gamma }_{p_i}=\sum _{j=1}^{N_{track}}{\dot{\gamma }}_{p_i}^j\times \delta t_j$$

(6-b)

In equation (6-a), ${\gamma }_{p_i}$ is the total deformation of the $i^{th}$ particle after $N_{track}$ steps of tracking, as introduced in equation (6-b). The parameter ${\dot{\gamma }}_{p_i}^j$ is the deformation rate of the particle at its current position, calculated by interpolating among the nodal values of the deformation rate magnitude inside each element, as read in equations (7-a) to (7-c).

$$\dot{\gamma }=\sqrt{\frac{1}{2}{\mathsf{{\rm I}}\mathsf{{\rm I}}}_{\mathsf{\Delta }}}$$

(7-a)

$${\mathsf{{\rm I}}\mathsf{{\rm I}}}_{\mathsf{\Delta }}=\sum _{i=1}^2\sum _{j=1}^2{\dot{\gamma }}_{ij}{\dot{\gamma }}_{ij}$$

(7-b)

$${\dot{\gamma }}_{ij}=\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}$$

(7-c)

Figure 6a,b illustrate the particle exchange between two sub-chambers of the mixer and the development of the distributive mixing by an increase of the mean average total strain, respectively. In these figures, we can quantitively observe the agreement between the results of two numerical strategies. Notably, particle exchange is not always an appropriate criterion to evaluate distributive mixing. For example, for a symmetric arrangement of the rotor, its value will be zero because there would be no particle exchange between two sub-chambers.

Figure 6a,b show that increasing the particle number ($n_{particle}=40000$) and decreasing the rotation increment ($∆\theta =1$) do not affect the results; however, reducing $∆\theta$ improves the extent of the particle loss, as shown in Figure 7.

Figure 7 demonstrates the particle loss for two numerical approaches, namely MLS-FEM and boundary-fitted mesh FEM. Ignoring the initial particle loss related to our particle distributing algorithm and not the numerical method, the MLS-FEM and the boundary-fitted mesh FEM exhibit approximately 15 percent particle loss. Decreasing the rotation increment (∆θ) from 5 to 1 significantly reduces the extent of particle loss, as depicted by the solid line in Figure 7.

Reproducing these results (∆θ=1) using the boundary-fitted mesh FEM would be time-consuming since it requires much more mesh regeneration as the rotation increment decreases. However, we have employed the MLS-FEM, which requires a fixed background mesh and makes it easy to reduce rotation increment. This advantage makes the MLS-FEM more efficient in lowering the computation time and accuracy (e.g., particle loss) than the boundary-fitted FEM.

6. Conclusions

We addressed the regeneration of meshes as a challenge in the numerical analysis of the mixing state of polymer compounds in mixers with moving parts (e.g., rotors). In this context, we presented different approaches using a fixed background mesh instead of multiple boundary-fitted meshes with drawbacks, such as being time-consuming. This point inspired us to introduce MLS-FEM, which uses a fixed background mesh to analyze the flow domain without introducing new degrees of freedom in the finite element formulation.

The moving least-squares aided finite element method (MLS-FEM) uses the MLS technique to enhance the existing shape functions of the FEM. Accordingly, the MLS-FEM considers solid components part of the flow domain and does not need a new mesh as the rotors' position, shape, and orientation change.

To validate our method, we comprehensively compared MLS-FEM's results with that of the boundary-fitted mesh FEM. The qualitative and quantitative similarities between the two approaches revealed our new method's excellent performance.