Lagrangian Differencing Dynamics for Time-Independent Non-Newtonian Materials

This paper introduces a novel meshless and Lagrangian approach for simulating non-Newtonian flows, named Lagrangian Differencing Dynamics (LDD). Second-order-consistent spatial operators are used to directly discretize and solve generalized Navier–Stokes equations in a strong formulation. The solution is obtained using a split-step scheme, i.e., by decoupling the solutions of the pressure and velocity. The pressure is obtained by solving a Poisson equation, and the velocity is solved in a semi-implicit formulation. The matrix-free solution to the equations, and Lagrangian advection of mesh-free nodes allowed for a fully parallelized implementation on the CPU and GPU, which ensured an affordable computing time and large time steps. A set of four benchmarks are presented to demonstrate the robustness and accuracy of the proposed formulation. The tested two- and three-dimensional simulations used Power Law, Casson and Bingham models. An Abram slump test and a dam break test were performed using the Bingham model, yielding visual and numerical results in accordance with the experimental data. A square lid-driven cavity was tested using the Casson model, while the Power Law model was used for a skewed lid-driven cavity test. The simulation results of the lid-driven cavity tests are in good agreement with velocity profiles and stream lines of published reports. A fully implicit scheme will be introduced in future work. As the method precisely reproduces the pressure field, non-Newtonian models that strongly depend on the pressure will be validated.


Introduction
Various types of fluids used in industry and in nature do not act in accordance with Newton's law, assuming that the fluid viscosity is constant. Non-Newtonian flows include mud, paint, blood, lymph fluid, cell fluid, milk, chocolate, edible oil, etc. Non-Newtonian fluids exhibit viscous properties on a regular basis, and it is important for the designer or engineer to be familiar with the flow behavior of such fluids, be able to identify the fluids' physical properties, and know how to use these properties to predict flow behavior in the industrial process. Viscoplastic and viscous fluids are both time-independent fluid categories, and novel numerical models belonging to these classes of non-Newtonian fluids will be discussed in this paper.
Purely viscous models simulate shear-thickening and shear-thinning fluids. Shearthinning fluids can be described by the Casson model [1], while the Power Law model [2] is a more generalized model. It is assumed that the Casson model has infinite viscosity for a zero shear rate, a yield stress below which the fluid does not flow, and zero value of the viscosity at an infinite shear rate. Shear-thinning behavior describes the melting of polymers, polymer solutions, biological fluids, mud and mayonnaise. The Power Law is a widely used, mathematically simple model that can approximately simulate the behavior of a non-Newtonian fluid. Depending on the flow-behavior index n, it can be classified solution is known, and for laminar flow inside a square cavity with a constant velocity lid, for which the analytical solution is unknown.
In Bognár and Csáti [21] paper, a spectral approach is used to address the lid-driven cavity flow problem and develop a MATLAB software to solve the Navier-Stokes equation. Bruneau and Saad [22] researched Hopf bifurcation, and the solution's behavior for intermediate and high Reynolds numbers, based on a finite-differences (FD) discretization and a multigrid solver with cell-by-cell relaxation. The findings of Coclite et al. [23] and Coclite [24], using time LBM in the rheology of elliptical particles, can provide specific insight into the rational design of micro-and nano-particles as drug carriers.
The Smoothed Particle Hydrodynamics (SPH) method is a Lagrangian and meshless method, which is widely used for a variety of fluid flow simulations [25][26][27]. Researchers and scientists are continually developing the SPH method, solvers and models according to the type of flow, although its traditional formulation suffers from an inaccurate pressure field [28]. The SPH has also been applied to non-Newtonian free-surface problems, such as mud and molding flows [29,30]. Hosseini et al. [31] presented a GPU implementation to achieve a better performance. Shao and Lo [30] simulated a dam-break problem and discussed the flow features of Newtonian and non-Newtonian flows using the Incompressible SPH (ISPH). Its advantage lies in the ease of tracking the free surface using a similar procedure to that employed in the the moving particle semi-implicit (MPS) method, and the findings correlated well with the experiments. Fan et al. [32] developed a matrix-free, implicit SPH solver for highly viscous non-Newtonian flows that contain high-pressure areas. The standard, explicit SPH method was not feasible due to the need for a very small time step to achieve a stable simulation. Zhu et al. [33] assessed how well the plastic viscosity can be determined using the SPH approach. Papanastasiou's Bingham constitutive model [34] was implemented in the SPH model and tested against the results of published data. Xu et al. [35] enhanced the SPH approach to 3D non-Newtonian flows with complex free-surface shapes and the Casson model. An artificial stress-term is inserted into the momentum equation to prevent tensile instability, which leads to the clustering of particles and non-physical defects in fluid stretching. Xenakis et al. [28] used a diffusion-based ISPH method to describe free-surface flows. The approach was developed to solve inelastic non-Newtonian flows by introducing a new viscous term that is more suitable for such flows.
The objective of the present work is to introduce a Lagrangian Differencing Dynamics (LDD) method for the simulation of non-Newtonian fluids, which combines the advantages of the meshless Lagrangian methods, such as SPH and MPS, with the advantages of implicit schemes, such as FVM or FEM. Peng et al. [36] applied LDD to granular flow modeling using a Drucker-Prager model, and found that its accuracy in terms of pressure and impact force is higher than SPH. This success promoted the idea of applying the LDD to simulate a wide range of non-Newtonian materials in this work. The LDD method relies on second-order consistent operators introduced by Basic et al. [37], and extends the scheme introduced by Bašić et al. [38]. The scheme solves the Navier-Stokes equations (NSE) using the split-step approach, where the pressure field is solved implicitly, and the velocity is solved explicitly. On the other hand, in this work, the scheme is extended to solve generalized NSE with variable viscosity, while solving the momentum equation in a semi-implicit manner. After obtaining the pressure and velocity fields, the LDD method relies on Lagrangian advection, which is shown to enable large time steps.
In this paper, the authors reproduced four benchmarks from different authors to validate the proposed methodology. The original benchmarks reproduced in this paper used different numerical methods and mathematical approaches to solve a specific set of problems for non-Newtonian and Newtonian fluids. The validation tests replicated in this study using the LDD method are as follows: dam break case originally validated using SPH and Bingham model [28], square lid-driven cavity simulated using the FVM and Casson model [11], skewed lid-driven cavity flow simulated using the FVM method and Power Law model [39] and fresh concrete Abram slump test problem validated using the PFEM and Bingham model [18].
The remainder of the paper is organized as follows. In Section 2, the mathematical model for simulating various viscoplastic materials is presented. Section 3.1 explains the numerical procedure for solving the introduced mathematical system. Section 4 presents the validation cases and discusses the numerically obtained results. Finally, the conclusions are drawn in Section 5.

Governing Equations
The generalized Navier-Stokes equations (NSE), in vector form, for the non-Newtonian incompressible fluid flow are described and solved. The conservation of momentum and mass is given as follows: where the advective derivative is expressed as D/Dt, the velocity vector as u, the fluid density as ρ, the fluid pressure as p, the stress tensor as τ, and F ext as the vector of external forces. Equations (1) and (2) imply temporal and spatial dependency. The stress tensor τ is defined for an incompressible fluid using the following expression: where µ is the dynamic viscosity of the fluid, and E is the strain rate tensor, which is defined as: where ∇u indicates the velocity-gradient tensor of the flowing material. The shear rate is defined as: where the colon (or double-dot) operator is defined as E : E T ≡ trace(EE T ), and trace represents the sum of the matrix diagonal elements.

Bingham Model
The Bingham model is a widely used, viscoplastic non-Newtonian model. This viscoplastic model is commonly used in engineering due to its mathematical simplicity, i.e., it is a two-parameter model. It is used in food, drilling, oil and gas, chemical and many other industries, as the majority of industrial fluids comply with Bingham's law. According to the Papanastasiou [34] model, the stress tensor is calculated as: where τ 0 is the yield stress, µ ∞ is the dynamic viscosity at infinite shear rate,γ = 2E, and m is the regularization parameter. The effective viscosity is calculated as follows:

Casson Model
The Casson model is a rheological model that is used to describe viscoelastic flow. It is expressed in accordance with the Papanastasiou [34] regularization: where τ 0 is the yield stress, µ ∞ is the dynamic viscosity at infinite shear rate and m is the regularization parameter. The effective viscosity is obtained as follows:

Power Law Model
Power Law is a widely used and mathematically simple model that can approximately simulate the behavior of a non-Newtonian fluid. In this generalized model, for purely viscous fluids, shear stress tensor is calculated as: the generalized model of the Power Law is defined by the effective viscosity in function of the shear rate, as follows: where µ 0 represents the flow consistency index, and n is the flow behavior index. Depending on the flow-behavior index n, it can mathematically model three types of fluids. The fluid describes shear-thinning flows for 0 < n < 1, Newtonian behavior for n = 1, and shear-thickening flows for n > 1. The zero-shear viscosity is approached at very low shear rates, while the infinite shear viscosity is approached at very high shear rates.

Lagrangian Differencing
The defined fluid domain is represented without topological information, using a cloud of points (i.e., nodes) described by the set Ω. At a given time, each node i ∈ Ω at location x i (t) in a specific time instant t carries some fluid properties. Each node i interacts with a set of neighbor nodes, N , which are found in the sphere with radius h, as shown in Figure 1. A symmetric and positive weighting function W(r, h 0 ), is used to evaluate the strength of the interaction between nodes, based on the distance between them, where r = x i − x j . The weighting function radius is represented by h 0 , and when r ≥ h 0 , the weighting function equals zero. The weighting function is defined as W(r, h 0 ) = (1 − r/h 0 ) 3 , where 0 ≤ r ≤ h 0 . Unlike in the SPH method, where the particles are constrained to a specific volume and the smoothing function must fulfill the quadrature requirements [27], in this study the weighting function is used only for the weighting of interactions in the finite-difference manner; therefore, it can take an arbitrary shape. This so-called Lagrangian differencing only accounts for a small number of immediate neighboring nodes [36], because a high accuracy is guaranteed due to the renormalization, which is introduced below. Consequently, the introduced method is more efficient when compared to classical SPH and MPS methods [36,38], which often employ a large compact radius to encapsulate a large number of neighbors because the convergence requirements of SPH approximation include particle spacing ∆ → 0, and ratio h/∆ → 0 [36]. In the LDD method, only the nearest neighbors need to be accounted for; therefore, h/∆ is taken as 1.6, significantly below the SPH standards. Values of some continuous function f (x) may be approximated anywhere between nodes by using Shepard's formula: where indicates a discrete version of some expression, f j ≡ f x j is introduced for compactness, and W j ≡ W x − x j , h 0 indicates the weight of the neighbor node j in the neighborhood. A renormalization tensor B i for the discrete gradient is defined as [37]: where (13), the discrete approximation of the gradient is obtained using the expression: The derivation of this expression for the gradient is given in [37], and the convergence theorems are presented in [40]. The method was validated to be consistent up to the second order, and robust for highly irregular arrangements of nodes.
An accurate discrete expression for the Laplacian operator used in this study belongs to a class of renormalized Laplacian operators introduced by Basic et al. [37], which were validated on particularly irregular neighborhood arrangements. The discrete Laplacian is defined as: where o i is the offset vector of the node i: which points from x i to the point where the arrangement of the neighborhood dominates. Bašić et al. [38] have shown that the renormalization process enhances the operator when approximating a scalar field Laplacian, and that it is responsible for reaching second-order accuracy when solving Poisson problems, while the original SPH and MPS formulations yield first-order accuracy. This operator is a crucial ingredient in discretizing the pressure and velocity equations.

Numerical Solver
In the split-step scheme, the pressure and velocity steps are solved separately. The pressure is obtained by solving the Poisson equation within the proper boundary conditions, and the solution must satisfy the continuity constraint. The pressure Poisson equation (PPE) is defined as: where n is the normal boundary. In Equation (17) the no-slip boundary condition along the Γ w interface is obtained by dotting Equation (1) with n. The external acceleration is taken as a constant; therefore, ∇ · F ext = 0 is excluded from Equation (17). Furthermore, ∇ · ∇ · τ is also taken as zero in Equation (17), since it was empirically tested and did not significantly affect the simulated problems. In other words, the spatial change in the stress divergence vector field does not generate a significant pressure gradient or affect the flowing material. This does not hold in general, but holds for the types of flows simulated in this study. The general impact of the term ∇ · ∇ · τ on the pressure and velocity will be analyzed in future work. The PPE was solved using the BiCGStab iterative solver in the matrix-free manner, as described by Bašić et al. [38]. After solving the PPE (17), the pressure gradient may be explicitly calculated using Equation (14), and the results may be substituted into the momentum Equation (1). The shear rate was calculated using Equation (5), also based on explicitly calculated velocity gradients for the strain rate Equation (4) using Equation (14). At this point, the velocity vector is the only unknown in the momentum Equation (1). Since the stress tensor τ is defined by Equation (3) through the strain rate E (Equation (4)), the stress tensor divergence may be written as ∇ · τ = µ ∇ 2 u + 2E ∇µ, and substituted back into Equation (1), which can be stated as: In this study, the momentum equation is solved for velocity in the semi-implicit formulation. The term with the velocity Laplacian is placed on the left-hand-side of the equation, while the term with the viscosity gradient is explicitly evaluated and placed on the right-hand-side. Furthermore, the Lagrangian acceleration term is discretized by the Euler approximation. The problem is, therefore, defined by the following expression: where n denotes the step with unknown variables, and n − 1 denotes the previous timestep with known values. The numerical method is built by directly discretizing equations using the discrete spatial operators described in Section 3.1. The analyses given in Basic et al. [37] have shown that the types of problems described using Equations (17) and (19) reach second-order convergence for Dirichlet, Neumann and mixed-boundary conditions. Equation (19), along with the no-slip boundary condition, can be straightforwardly solved using an iterative matrix-free solver.
After obtaining the velocity field by solving u n in Equation (19), mesh-free nodes are advected in the Lagrangian manner. Second-order Euler approximation is applied to advance positions of mesh-free nodes, using the velocities of the previous and currently obtained time-steps. Therefore, the newly advected position of a node i is explicitly obtained by calculating the expression: The Lagrangian solution dictates that each node advects along its streamline. Since the size of the time-step is finite and the flow is often complex, advecting along the streamline tangent can lead to the eventual separation or collision of nodes' trajectories. The possibility of distorted node arrangements and errors in conservation of the volume are resolved by introducing a reordering step after the advection of mesh-free node neighbourhoods.
The reordering step employs the Position-Based Dynamics (PBD) method introduced by Macklin and Müller [41], which iteratively displaces node neighborhoods to enforce constant density in the fluid. This results in equally spaced nodes in the point cloud, thereby enforcing incompressibility and allowing the utilization of larger time-steps. The velocities at reordered locations are interpolated using Equation (12) only for excessive time-step sizes and significant reorderings, as described by Bašić et al. [38]. Generally, using reasonable time-step sizes does not result in significant reordering, but a minor movement of the node position within its radius to maintain the equidistant property.
In should be noted that the boundary surfaces of the domain and body boundary surfaces are required when imposing boundary conditions in Equations (17) and (19), but no preparation of any kind of mesh is required. Each time-step fluid nodes that close to the walls is projected onto the geometry, as shown in Figure 2. The newly projected nodes impose Neumann (pressure gradient) boundary conditions in the PPE or Dirichlet (no-slip) boundary condition in the momentum equation. The projections guarantee that the closest fluid nodes are aligned along their (geometry) normals. This property of the LDD method allows for the accurate imposition of the Neumann boundary conditions.

Square Lid-Driven Cavity Flow
Lid-driven cavity tests are popular benchmarks for the newly introduced CFD solvers and methods [42]. This test's simple geometry yields complex and distinct types of flows; therefore, it can provide a better understanding of the industrial complex processes in closed recirculating regions. While a strong extension occurs near the lid edges, the rotational flow can be seen in the center of the recirculating region. The ideally moving lid meets the stationary wall due to idealization, which results in a discontinuous velocity. It is hard to obtain the full range of kinematics, as well as rapid changes in pressure and stress near the corners, than it first appears. It is important to have a suitable method of calculating convection-dominated momentum transfer. For the above reasons, lid-driven flow in a cavity is a good initial experiment to validate flows for various Reynolds numbers and flow properties.
Casson fluid flow was simulated by the two-dimensional, steady-lid-driven cavity problem. The LDD numerical simulations were validated against the simulations conducted by [11]. The cavity space in dimensionless units is 1 × 1, and the moving lid had a steady, dimensionless velocity of u lid = 1. The no-slip condition was applied to the lid and wall boundaries. Two point cloud resolutions were tested. The time-step used for both resolutions was δt = 10 −3 , and the average time of the time-step calculation was 12 ms on a modern GPU. 10 s of actual physical time was simulated. The Reynolds number considered was Re CA = 100 and the Bingham number Bn = 0.01, with the Reynolds and Bingham numbers, respectively, defined as follows: where ρ is the fluid density, U ∞ is the velocity at infinite shear rate, l is the reference length, µ ∞ is the dynamic viscosity at infinite shear rate. Figure 3 renders the velocity magnitude and streamlines. The obtained results agree very well with those presented by Neofytou [11]. The plotted streamline pattern in Figure 3 showed that the center of the primary vortex was successfully and accurately captured, while the secondary vortices were also clearly reproduced in the bottom corners of the cavity box. Figure 4 plots the velocity profiles perpendicular to the vertical and horizontal centerlines of the cavity. Two simulations with different initial resolutions, 100 × 100 and 200 × 200, were performed to prove the convergence of the method. The u-component of the velocity for the finer resolution is in excellent agreement with the FEM method. The coarser resolution had slightly underestimated values, but overall agreed well with the referent results. The v-velocity plotted in the horizontal centerline of the cavity for both resolutions is in very good agreement with FEM. The only discrepancy occured at its highest and lowest peak, where the v-component of the velocity showed slightly lower absolute values, even though it agrees very well in all other areas.

Skewed Lid-Driven Cavity Flow
The Power Law viscosity model was tested for a fluid circulating in a lid-driven skewed cavity flow. Demirdžić et al. [43] was the first to publish results for skewed angles, α = 30 • and 45 • for Newtonian fluids. Thohura et al. [39] compared their results with the published results and extended the investigation to various skew angles. An experiment of the skewed cavity using the Power Law, reported by Thohura et al. [39], is reproduced in this paper using the LDD method, and the results are compared.
The circulation pattern and vortex formation are highly dependent on the Reynolds number for any rheological behavior. Therefore, a relatively high Reynolds number was chosen as the simulation, Re = 500, to show the stability and robustness of the LDD method. The Power Law index was n = 1.5. The cavity in dimensionless units had a size of 1 × 1, and the angle of the side wall to the baseline was α = 60 • . The lid moved with a steady, dimensionless velocity of u lid = 1. The no-slip condition was applied to the lid and wall boundaries. Three initial resolutions of 50 × 50, 100 × 100 and 200 × 200 nodes were tested. The time step used for simulations was δt = 10 −3 , and the calculation of a step took 26 ms in average. 20 s of actual physical time was simulated until steady simulation was reached.
In Figure 5, plotted streamlines correspond very well to the reference data provided by Thohura et al. [39]. The positions of the vortices are correctly captured, and vortex in the lower-right corner is more clearly represented in the LDD method than in the FVM. Figure 6 compares the simulated results of the u and v components of the velocity along the vertical and horizontal centerlines of the cavity with the referent data. For the u-component of the velocity, the plotted curves match almost perfectly, while v-component of the velocity has small discrepancy up to the height of Y = 0.65. Overall, a good match between the obtained results is achieved, which proves that the method is capable of simulating non-Newtonian Power Law fluids.

Dam Break of a Bingham Fluid
A dam break test was performed in accordance with the experimental parameters specified by Komatina and Jovanovíc [44]. A non-Newtonian fluid was represented by a mixture of water and clay in a reservoir released onto a channel with an initial length of 2 m and initial height of 0.1 m. The mixture density of water and mud was ρ = 1200 kg/m 3 , the Bingham yield stress was τ B = 25.0 Pa and the viscosity was µ B = 25.0 Ns/m 2 . The no-slip boundary condition was used along the channel walls. The simulation snapshots were taken at times t = 0.1, 0.3, 0.6, 1.0 s, and rendered in Figure 7. The results of the LDD method were compared to the experimental data and other numerical methods in Figure 8. Shao and Lo [30] conducted the dam-break test by employing the Incompressible SPH (ISPH) method, while Cremonesi et al. [16] tested a Lagrangian formulation for weakly compressible fluids by employing an explicit solver based on PFEM. The simulated flow visually and numerically corresponded to the experimental data. The distinctive characteristics of the flow could be recognized. At the start of the flow, the free-standing end of the fluid column flowed, and the upper free corner began to collapse. Subsequently, the fluid tended to decrease in height as it flowed, obtaining peak velocity at the flow front. Towards the end, the flow behaved as a creeping flow and the surface profile was almost unchanged, which corresponds with the experimental findings. From the time instance t = 0.0 ÷ 0.2 s all the numerical methods showed a discrepancy that can be explained by the lack of a vertical wall in the channel that releases the flow. The numerical results of the LDD method are closer to the experimental results between t = 0.2 ÷ 0.8 s times than the other numerical methods. From t = 0.8 ÷ 1.2 s, the FEM method [16] and the LDD method both show a good fit with the experimental data. The results obtained by the LDD method clearly show the ability of the novel method to simulate Bingham-type flows and reproduce realistic results. A convergence study was performed using two simulations with different resolutions and time-step sizes. The first simulation used 200,000 nodes (∆ = 1 mm) and the second simulation used 50,000 nodes (∆ = 2 mm). For both resolutions, the tested time step sizes were δt = 10 −3 and δt = 10 −4 . Table 1 shows the comparison of the tracked mass front position, obtained numerically and experimentally. The results are plotted in Figure 9. For the domain of 50,000 nodes, the average calculation time for one time-step δt = 10 −3 took 35 ms, while for the domain of 200,000 nodes and time-step δt = 10 −4 , it took 200 ms. The results of the convergence study show that the domain of 200,000 nodes and the time-step δt = 10 −4 gave the best results, and this was, therefore, used for the dam break comparison.

Fresh Concrete Slump Test
As a standard laboratory experiment, the slump test was used to determine the workability of fresh concrete. The slump test will be reproduced for the purpose of validating the method in three dimensions. Franci and Zhang [18] provided a detailed comparison of the experimental data and their numerical simulations.
The conical container was filled with concrete and the evolution of the form was measured after the container was removed. The spread and the slump of the evolving concrete were measured, i.e., fluid height and diameter. The test was completed once there was no movement in the fluid. This so-called Abram's test [45] was computed using a single-phase Bingham model, with the regularization parameter m = 1. Geometrical and material data are given in Table 2. The computation was made using 50,000 nodes in the domain and an average time-step took 55 ms on the GPU. Time-step value was set to δt = 10 −3 s, and 40 s of physical time was simulated. To test the stability of the solver, the experiment was tested using larger time-steps, up to δt = 10 −2 s, which remained stable.
The spread of cement was computed with good accuracy, as shown in Figure 10. During the time interval t = 0.0 ÷ 10.0 s, both PFEM and LDD numerical simulations resulted in approximately the same evolution of the diameter. The diameter values for those time instances are higher than in the experiment. After t = 10.0 s the results of the LDD and PFEM began to differ. They both showed a smaller diameter than the experiment, but the LDD method asymptotically reached the experimental results. At the time instance t = 40.0 s, the LDD method reached the experimentally obtained diameter, while the PFEM method still showed a lower value. The simulation snapshots presented in Figure 11 show the flow results for three time instances, t = 0.5, 5.0, 40.0 s. The axisymmetric flow reached maximum velocity at the top of the free surface, immediately after the container was removed. As the top slowly drops down, and the diameter evolves, the velocity decreases. After the cone completely collapses, the diameter slowly increases, until the final timepoint, where fluid has almost zero movement. These results are in agreement with those presented by Franci and Zhang [18].

Conclusions
In this paper, a novel meshless and Lagrangian method for 2D and 3D non-Newtonian problems with free-surface is introduced, named Lagrangian Differencing Dynamics (LDD). The method uses Lagrangian differences derived from finite differences to discretize and solve generalized Navier-Stokes equations. The fluid domain is represented by a cloud of points, in which each point advects in Lagrangian manner, while points near walls are projected to impose boundary conditions. Therefore, it is possible to simulate complex flows and complex shapes using relatively large time-steps. Due to the second-order consistent operators, and implicit solving of the pressure and velocity fields, the time-steps used were relatively large. In addition, parallelization on the CPU and GPU, as well the Lagrangian movement, allowed for simulations to be carried out relatively quickly.
Four tests with non-Newtonian fluids were validated using the Bingham, Casson and Power Law models. The results of the simulations are in good agreement with data from the corresponding experiments. The accuracy, high speed and robustness of the method are shown through different tests, parameters and flow types. For all of the performed benchmarks, the convergent behavior of the solver was verified. The velocity and pressure solvers showed stable and convergent behavior. It was noted that the regularization parameter m has an impact on the results, which will be investigated in future work.
Square and skewed cavity simulations were conducted in accordance with experiments and numerical simulations from the literature. The square lid-driven cavity test was simulated using the Casson model for Re = 100, while the skewed lid-driven cavity was simulated using the Power Law model for Re = 500. The problem of discontinuous velocity where the moving lid meets the stationary wall was determined, as well as the recirculation in the cavity. The full spectrum of kinematics and pressure fields was achieved. The distinctive characteristics of the flow were recognized in the dam break test simulated with the Bingham model. The Abram slump test was successfully reproduced and the evolution of the form agreed well with the experimental data.
Instead of the presented semi-implicit scheme, a completely implicit scheme for the velocity will be introduced in future work. Since the method is capable of reproducing accurate pressure fields, it will be validated by non-Newtonian models that significantly depend on the pressure. Future work will also address the validation of problems of Fluid-Structure Interaction (FSI). The current implementation of the introduced methodology is based on uniformly spaced nodes in the point-cloud, but the introduced operators handle non-uniform arrangements. Therefore, in future work, a local refining of the point cloud will be introduced to obtain more accurate solutions within turbulent areas of flow.