A physically consistent particle method for high-viscous free-surface flow calculation

Particle methods for high-viscous free-surface flows are of great use to capture flow behaviors which are intermediate between solid and liquid. In general, it is important for numerical methods to satisfy the fundamental laws of physics such as the conservation laws of mass and momentum and the thermodynamic laws. Especially, the angular momentum conservation is necessary to calculate rotational motion of high-viscous objects. However, most of the particle methods do not satisfy the physical laws in their spatially discretized system. The angular momentum conservation law is broken mostly because of the viscosity models, which may result in physically strange behavior when high-viscous free-surface flow is calculated. In this study, a physically consistent particle method for high-viscous free-surface flows is developed. The present method was verified, and its performance was shown with calculating flow in a rotating circular pipe, high-viscous Taylor–Couette flow, and offset collision of a high-viscous object.


Introduction
Particle methods are the numerical methods which have an advantage in handling large deformations of free surface flow. The typical ones are the weakly compressible Smoothed Particle Hydrodynamics (SPH) method developed by Monaghan [1] and the strictly incompressible Moving Particle Semi-implicit (MPS) method developed by Koshizuka and Oka [2]. Particle methods were mainly applied for the low viscosity flows such as breaking waves [3]. Besides, it is expected that high viscosity fluid calculations will lead to the broader application of particle methods. It is because a technique for high-viscous calculation can be applied not only to the problems with high-viscous fluid [4][5][6][7][8] but also to the following ones where solid and liquid coexist: (1) Problems with melting and solidification [9][10][11][12][13][14][15][16][17][18].
Since these solid-liquid coexistences are often expressed by the change in viscosity, various expression in such material flows will be enabled if the broad range of viscosity can be handled in a unified manner.
For high-viscous fluid calculation using particle methods, numerical stability against high viscosity and physical consistency such as conservation of angular momentum are important. The former is effectively achieved by implicit velocity calculation with respect to viscosity term. When the particle velocities are updated explicitly, the time step width for the numerical stability will be smaller and it results in an inefficient calculation. In such cases, it is difficult to calculate very large viscosity equivalent to solid state. In fact, the implicit calculations were adopted for high-viscous calculation in the previous studies [5,6,[9][10][11][12][13][14][15][16]20], and recently, Monaghan [7] proposed pair by pair implicit calculation algorithm which can avoid solving large matrix equations.

3
The other thing is physical consistency, which is generally important for numerical method. Specializing for particle methods, it is called physically consistent when the system obtained from the continuum governing equations via space discretization satisfies the fundamental laws of physics (i.e. mass conservation, momentum conservation, angular momentum conservation, and the thermodynamic laws).
To satisfy the thermodynamic laws, analytical mechanical frameworks are helpful. In fact, some particle methods were constructed based on the frameworks. For example, the general equation for non-equilibrium reversible-irreversible coupling (GENERIC) framework [32] was used by Español et al. [33] in developing smoothed dissipative particle dynamics (SDPD). The classic analytical mechanics for energy conserving system was applied by Ellero et al. [34] and Suzuki et al. [35] in developing incompressible particle methods. Specifically, they are incompressible smoothed particle hydrodynamics (ISPH) and Hamiltonian moving particle semi-implicit (HMPS), respectively. The energy conserving framework was also used by Price [36] in proposing a particle method with variable smoothing length. Besides, the extended Lagrange mechanics for the dissipating system [37] was adopted by Kondo [38][39][40] in developing a physically consistent particle method including dissipation. In this study, the method in the previous studies [38][39][40] is termed moving particle hydrodynamics (MPH) method since its formulation is intermediate of the SPH [1] and MPS [2] methods.
The analytical mechanical frameworks can take the thermodynamic laws into consideration, but it does not always ensure satisfying the mechanical conservation laws with respect to linear momentum and angular momentum. As Weiler et al. [16] pointed out, the angular momentum conservation plays an important role for calculating rotational motion of highviscous objects with free-surface. However, when the viscosity term was discretized using the difference-based Laplacian models [2,9,18,[41][42][43], which are often used in SPH and MPS, the rotational motion is unphysically suppressed because of the torque against rotation. In some researches, the gradient model was applied twice [6, 19, 20, 22-24, 26, 27, 30] instead of directly applying the difference-based Laplacian model. The velocity gradient tensor was calculated by the first gradient model, and the divergence of the tensor was calculated by the second one. Such models could improve the calculation of rotation, and moreover, the angular momentum could be conserved when the gradient correction [44,45] was applied to each gradient model [25,28,29,31]. However, the gradient correction in addition to the nest of the gradient models made the formulation complexed, and it hindered the implicit velocity calculation with respect to the viscosity term. On the other hand, Kondo et al. [15,46] added a new unknown parameter, i.e. angular velocity, to conserve angular momentum by subtracting rotational element from the difference-based Laplacian model. The formulation was simpler than that with the nested gradient models, and it could adopt the implicit velocity calculation [15]. However, the additional parameter was not favorable because it increased the degree of freedom in the implicit calculation.
Compared to the above two methodologies, there is a simpler way for angular momentum conservation. It is just introducing damping force between each particle pair. This pairwise duping was originally used as artificial viscosity [1,47], but was diverted to calculate practical viscosity of the fluid [48][49][50] and was also applied in high-viscous calculations [4,7,16,17,21,25]. Since the formulation of the pairwise damping model is simple, it is also suitable for the implicit calculation. In fact, Weiler et al. [16] and Monaghan [7] applied the implicit viscosity calculation with the pairwise damping model. However, the pairwise damping model was not clearly related to the velocity gradient and stress tensor which are often used in solid analysis.
From the viewpoint of numerical stability against high viscosity and the physical consistency such as angular momentum conservation, the high viscosity fluid calculation of Weiler et al. [16] is currently advantageous. However, their model is not clearly related to an analytical mechanical framework. Moreover, they iteratively solved the strict incompressible condition, but it is not always required for the high viscosity calculation. The strict incompressibility is needed when high pressure is expected as in molding simulation [4,20], but the weakly compressible approach can save numerical resources, for example, in the cases where the force acting on the fluid is mainly just gravity. The MPH method [38][39][40], which is physically consistent, has both strictly incompressible version (MPH-I) [38] and weakly compressible version (MPH-WC) [39,40], and the latter can be efficient in such cases.
In this study, a physically consistent particle method for high-viscous fluid flow is developed by introducing the pairwise damping viscosity model [48][49][50] to the MPH-WC method [39,40]. To avoid the restriction of the time step width due to high viscosity, the implicit velocity calculation is adopted. Moreover, the stress tensor evaluation procedure using the particle interaction force is proposed, and the stress tensor corresponding to the pairwise damping model is clarified. To verify the present method and to show its performance in low Reynolds number condition, flow in a rotating pipe, high-viscous Taylor-Couette flow and offset collision of a high-viscous object are calculated, with comparing the results to those obtained using the difference-based viscosity model, which does not conserve angular momentum.

Discretization of governing equation
Governing equations adopted in this study are Navier-Stokes equation and an equation for pressure calculation where ρ, u, Ψ, μ, g, λ and κ are density, velocity, pressure, shear viscosity, gravity, bulk viscosity and bulk modulus, respectively. In this study, the weakly compressible approach [39,40] was adopted, and the incompressible flow was simulated by using large values for bulk viscosity λ and bulk modulus κ. In the MPH method [38][39][40], the governing equations are discretized using effective radius and weight function in a similar way as in the SPH [1] and MPS [2] methods. In this study, the weight function is given as where r e is the effective radius, and r ij is the relative position vector between particle i and j. Particle interaction models for gradient, divergence and Laplacian operators are given as respectively. The characters ϕ and A are arbitrary scalar and vector, whose subscripts i and j indicate that the parameters are of particle i and j. The vector e ij is a unit vector from particle i to particle j, which is given by The parameters w ′ ij and S w ′ are the slop of weight function with negative sign defined as and the normalization parameter with respect to w ′ ij , respectively. The gradient model in Eq. (4) will be consistent to the pressure evaluation based on the virial theorem [51] when S w ′ is given by using a well-arranged uniform particle distribution, where d is the special dimension. To avoid arbitrariness of particle arrangement in calculating S w ′ , analytical integration was used instead in this study. Specifically, in two dimensional case, S w ′ is given by The interaction models (Eq. (4)) adopted in the MPH method [38][39][40] are basically equivalent to those in the SPH [1] and MPS [2] methods. Although they have zeroth order accuracy [38], it is expected that they can approximate the differential operators in the governing equations.
The Navier-Stokes equation (Eq. (1)) is discretized with the interaction models (Eq. (4)) as where u ij = u j − u i is the relative velocity from particle i to particle j. The second term on the right-hand side in Eq. (10) is a discretized version of the shear viscosity term. However, this type of difference-based Laplacian model [2,9,18,[41][42][43] does not conserve angular momentum and unphysically hinders rotational motion. It is because the direct use of the relative velocity u ij yields the torque against rotation. One of the alternatives which conserve angular momentum conservation is the pairwise damping model [48][49][50]. The model is originally introduced for artificial viscosity [1,47] in the SPH method, but it is also applied to calculate practical viscosity [48][49][50]. Moreover, Hu et al. [50] theoretically showed the relation between the difference-based model and the pairwise damping model. In this study, the differencebased model in Eq. (10) was converted to a pairwise damping model for the angular momentum conservation as based on the relation shown by Hu et al. [50], and the two models (Eqs. (10) and (11)) were compared. The other governing equation (Eq. (2)), which is for pressure calculation, is discretized as In discretizing the second term on the right-hand side of Eq. (2), the continuum equation is used as where n i is particle number density defined as and n 0 is a base value calculated using a well-arranged uniform particle distribution. To avoid tensile instability [52], the second term on the right-hand side of Eq. (12) was ignored when n i < n 0 in this study.

Stress evaluation based on virial theorem [51]
Since the system in the MPH method [38][39][40] is physically consistent, the pressure in the system can be evaluated via the virial theorem [51]. Specifically, the virial pressure P Φ in the test region Φ is evaluated using the interaction force F ij and the relative position r ij as Assuming the test region having single particle size ΔV around particle i, Eq. (15) can be re-written as where P i is the virial pressure at particle i. For evaluating stress as well as pressure, Eq. (16) is extended as where σ i is the stress tensor at particle i. This is a natural extension of Eq. (16) because its trace divided by the special dimension d will get back to Eq. (16). The interaction forces with respect to the difference-based viscosity model in Eq. (10) and the pairwise damping viscosity model in Eq. (11) are and respectively. Here, the stress tensor regarding each viscosity model is going to be derived by substituting the interaction forces (Eqs. (18) and (19)) into Eq. (17), assuming that the particle distribution is uniform and the velocity gradient is constant. The derivation in detail is shown in "Appendix". When the difference-based model is adopted, the relation between the stress tensor σ and the velocity gradient L is On the other hand, when the pairwise damping model is adopted, the σ-L relation is The stress tensor regarding the difference-based model (Eq. (21)) is not always symmetric, but that of the pairwise damping model (Eq. (22)) is symmetric. This is corresponding to the property of the models with respect to angular momentum conservation. Although the stress tensor derived with the pairwise damping model has a bulk viscosity term, which is the second term on the right-hand side of Eq. (22), it is negligible when the fluid is nearly incompressible and tr(L) ≈ 0.

Physical consistency and analytical mechanics with dissipation [37]
For the reliability of the numerical method, it is important to satisfy the fundamental laws of physics also in its discrete system. The analytical mechanical frameworks are useful in satisfying the mass conservation law and the thermodynamic laws in a discrete calculation of particle methods. Since the method in this study includes energy dissipation term, the extended Lagrange mechanics for the system with dissipation [37] is applied. In the framework, the Lagrange equation is given as where L is the Lagrangian defined by the difference of kinetic energy T and potential energy V as and D is the Rayleigh dissipation function [37]. When the Lagrangian L and the dissipation function D are given as and respectively, the Lagrange equation (Eq. (23)) agrees with the discrete governing equations (Eqs. (11) and (12)). Therefore, it is assured that the mechanical energy in this method monotonically decreases following the second law of thermodynamics.

Time integration
To calculate motion of the particles following the spacediscretized governing equations (Eqs. (11) and (12)), time integration is to be considered. In this study, after explicitly calculating the pressure based on the weakly compressible approach [39,40] as the velocity was implicitly calculated to handle large shear viscosity as Then, the position was updated as where upper index k indicates the time steps, and Δt is the time step width. Although the velocity in the shear viscosity term (the second term on the right-hand side of Eq. (28)) was implicitly calculated, the velocity in the bulk viscosity term (the first term on the right-hand side of Eq. (27)) was explicitly calculated. This is because the large pressure load is not expected in this study and the explicit calculation is numerically efficient in such cases.

Flow in a rotating circular pipe
Flow in a rotating circular pipe is calculated with the difference-based viscosity model (Eq. (10)) and the pairwise damping viscosity model (Eq. (11)), and the results in the steady state are compared. The calculation geometry is shown in Fig. 1, and the calculation condition is shown in Table 1. The angular velocity around the pipe center with the two model is shown in Fig. 2. Since there was almost no velocity in radial direction, the angular velocity was evaluated using the position from the pipe center x i and the velocity u i as Theoretically, the angular velocity of the fluid agrees with that of the pipe. However, it spuriously exceeded the theoretical value in the region close to the inner free surface when the difference-based model was adopted (Fig. 2a). On the other hand, when the pairwise damping model was adopted, the theoretical value appeared in the whole domain (Fig. 2b).
With each model, the circumferential stress acting on the surface with radial normal σ rθ , and the radial stress acting on the surface with circumferential normal σ θr were evaluated. The stresses σ rθ , and σ θr at particle i were calculated with the stress tensor (Eq. (17)) as and respectively, where e r and e θ are the unit vectors in radial and circumferential directions. The stresses σ rθ , and σ θr in the case with the difference-based viscosity model were shown in Fig. 3a. Theoretically, the two stresses, σ rθ , and σ θr , are the same due to angular momentum conservation law, and there is no stress in this case where the fluid is to rotate like a rigid body following the pipe rotation. However, σ rθ had positive distribution and σ θr had negative distribution (Fig. 3a). These non-zero different distributions indicate that unphysical stress against the angular momentum conservation law may emerge when the difference-based viscosity model is adopted. On the other hand, with the pairwise damping model, the stresses σ rθ , and σ θr were almost zero in the whole calculation domain (Fig. 3b), which indicates that the fluid just rotated like a rigid body as in the theory.

High-viscous Taylor-Couette flow
The Taylor-Couette flow problem is a flow between two rotating pipes that share a center. Since the angular velocity distribution and stress distribution are theoretically given when laminar flow is assumed, it can be used to verify the viscosity model related to rotation. Specifically, when the radii and angular velocities of the outer and inner circles are denoted by r 1 , ω 1 , r 2 , ω 2 , respectively, the analytical solution for the distribution of angular velocity and stress are given as and respectively. The geometry for the Taylor-Couette flow calculation is shown in Fig. 4, and the calculation condition is shown in Table 2. With this condition, the Reynolds number is low enough to be laminar. The difference-based viscosity model (Eq. (10)) and the pairwise damping viscosity model (Eq. (11)) are compared. The radial distribution of the angular velocity obtained by the two model are shown in Fig. 5. Regarding the angular velocity distribution, the calculated solution was close to the theory with both viscosity models. Besides, the radial distribution of the stresses σ rθ , and σ θr obtained with the two viscosity models are shown in  Fig. 6. With the difference-based viscosity model (Fig. 6a), the distributions of stresses σ rθ , and σ θr were different from each other and they were apart from the theoretical solution. On the other hand, with the pairwise damping model (Fig. 6b), the distribution of the stresses σ rθ , and σ θr were the same and mostly followed the theory. However, the deviation remained even in the case with the smaller particle size. This is possibly because of the non-uniform particle distribution in a practical calculation. In deriving the σ-L relation of Eq. (22), it is assumed that the particle distribution is uniform, but the radial distribution is not uniform especially when the particle distance r is smaller than the particle size l because particles do not get very close to each other. Since the region having non-uniform distribution will be relatively smaller with the larger effective radius, the deviation due to   the radial distribution becomes the smaller with the larger radius. The stress distribution with varying the effective radius, r e = 0.025, 0.035, 0.045 m, is shown in Fig. 7, where the particle size l = 0.01 m is applied. It is confirmed that the closer results to the theoretical solution was obtained with the larger radius.

Offset collision of a high-viscous object
When an object having large viscosity collides with an offset, it usually rotates due to the inertial force. This example is for investigating the ability of the numerical method to capture such rotational behavior. The initial state and calculation condition are shown in Fig. 8 and Table 3, respectively. For the comparison, the calculation is conducted using both the difference-based (Eq. (10)) and pairwise damping (Eq. (11)) viscosity models. The results using the two models are shown in Fig. 9, respectively. In the figures, the color contour shows the Mises stress calculated by where σ xx , σ xy , σ yx and σ yy , are the elements of the stress tensor (Eq. (17)). Since the calculation is two dimensional, the z element is assumed zero. When the difference-based viscosity model is applied, the rotation after colliding the obstacle was unphysically suppressed (Fig. 9a). It is because the torque against rotation emerged in the model. On the other hand, the natural rotation was observed with the pairwise damping model (Fig. 9b).
The calculation results in Sects. 3.1, 3.2 and 3.3 indicate that the angular momentum conservation is of importance in capturing rotational behaviors of high-viscous material. In this study, the pairwise damping model was adopted representing the models conserving angular momentum, but there are other viscosity models with such property. For example, the models which approximate strain rate tensor and divergence of stress tensor with using nested gradient models [6, 19, 20, 22-24, 26, 27, 30] can conserve angular momentum when gradient correction [44,45] is properly applied [25,28,29,31]. Another example is an extension of the difference-based model [15,46], where the rotational element is subtracted using an additional parameter, i.e. angular velocity. It is expected that these alternative models conserving angular momentum also perform well in calculating the rotation of high-viscous flows. However, they would need more calculation cost than the simple pairwise damping model especially in the implicit velocity calculation to cope with very high viscosity.

Conclusion
In this study, a physically consistent particle method for high-viscous fluid flow was developed by introducing the pairwise damping viscosity model [48][49][50] to the MPH-WC method [39,40]. To avoid the restriction of the time step width due to high viscosity, the implicit velocity calculation was adopted. Moreover, the stress tensor evaluation equation using the particle interaction force was proposed. With the equation, the relation between stress and velocity gradient were derived. The stress tensor derived with the difference-based viscosity model was not symmetric, but the one derived with the pairwise damping model was symmetric, reflecting the property of each model with respect to angular momentum conservation.
To verify the present method and show its performance in low Reynolds number condition, flow in a rotating pipe, high-viscous Taylor-Couette flow and offset collision of a high-viscous object were calculated, with comparing to the difference-based viscosity model which does not conserve angular momentum. In the calculation of flow in a rotating pipe, the difference-based model suffered spurious velocity and unphysical stress, but the pairwise damping model could calculate the fluid rotation following the pipe as in the theory. In the Taylor-Couette calculation, the velocity distribution agreed with the theory using both models. However, the stress distribution was apart from the theory with the difference-based model, while it agreed well with the theory with the pairwise damping model. In the calculation of offset collision of a high-viscous object, the rotational behavior was unphysically hindered by the

Appendix: Relation between stress tensor and velocity gradient
Here, the relation between the stress tensor σ and the velocity gradient L is to be derived based on the stress evaluation equation (Eq. (17)). The derivation is shown in two dimensions, but the derived relations (Eqs. (21) and (22)) will be the same in three dimensions. When the velocity gradient L is assumed constant, the relative velocity u ij is related to the relative position r ij as Substituting the interaction force with respect to the difference-based model (Eq. (18)) into Eq. (17) and applying Eq. (36), is obtained, where e ijx and e ijy are the x and y elements of e ij . Assuming the uniform particle distribution, the summation in Eq. (37) can be replaced with the integration as L 00 e ijx e ijx + L 01 e ijx e ijy L 00 e ijx e ijy + L 01 e ijy e ijy L 10 e ijx e ijx + L 11 e ijx e ijy L 10 e ijx e ijy + L 11 e ijy e ijy Substituting the normalization parameter S w ′ in Eq. (9) for Eq. (38), the σ-L relation (Eq. (21)) is obtained. The similar derivation can be done with respect to the pairwise damping model. Replacing the interaction force in Eqs. (17) to (19), and using Eq. (36), the stress is expressed as With assuming the uniform particle distribution, the expression can be replaced with the form with integral as L 00 cos 2 + L 01 cos sin L 00 cos sin + L 01 sin 2 (L 00 e ijx e ijx + L 01 e ijx e ijy + L 10 e ijx e ijy + L 11 e ijy e ijy ) e ijx e ijx e ijx e ijy e ijx e ijy e ijy e ijy . Finally, by substituting Eq. (9) into Eq. (41), the σ-L relation (Eq. (22)) is derived.
Funding Not applicable.

Availability of data and material
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Conflict of interest
The author reports no conflict of interest in this study.

Code availability
The code that supports the findings of this study is available from the corresponding author upon reasonable request and licensing.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.