GSPH with viscosity

GodunovSPH with shear viscosity : implementation and tests


The acceleration and energy dissipation terms due to the shear viscosity have been implemented and tested in GodunovSPH. The double summation method has been employed to avoid the well known numerical noise of the second derivative in particle based codes. The plane Couette flow with various initial and boundary conditions have been used as tests, and the numerical and analytical results show a good agreement. Not only the viscosity–only calculation, but the full hydrodynamics simulations have been performed, and they show expected results as well. The very low kinematic viscosity simulations show a turbulent pattern when the Reynolds number exceeds . The critical value of the Reynolds number at the transition point of the laminar and turbulent flows coincides with the previous works approximately. A smoothed dynamic viscosity has been suggested to describe the individual kinematic viscosity of particles. The infinitely extended Couette flow which has two layers of different viscosities has been simulated to check the smoothed dynamic viscosity, and the result agrees well with the analytic solution. In order to compare the standard SPH and GodunovSPH, the two layers test has been performed again with a density contrast. GodunovSPH shows less dispersion than the standard SPH, but there is no significant difference in the results. The results of the viscous ring evolution has also been presented as well, and the numerical results agrees with the analytic solution.

Stars – stars: novae, cataclysmic variables, methods: numerical

1 Introduction

Smoothed Particle Hydrodynamics (hereafter SPH) is a widely used method in numerical astronomy since its first announcement (Gingold & Monaghan, 1977; Lucy, 1977), because it is practically the only multidimensional Lagrangian technique. The Eulerian remap process is not necessary in SPH, so every particle can remember its own evolution history. Furthermore, by using a kernel with compact support and the tree algorithm (Barnes & Hut, 1986), the SPH method is able to achieve a reasonable speed in finding neighbors and the calculation of self-gravity.

In spite of its various advantages, SPH has some problems. The long standing one may be the side effect of the artificial viscosity. The artificial viscosity implemented in the standard SPH (Monaghan, 1992) is very effective to capture shockwaves and to prevent the particle penetration in the complicated fluid motion, but has a side effect in a velocity shear. The problem is that the bulk and shear components cannot be decoupled in the artificial viscosity of SPH, and the exact contribution of the shear viscosity in the numerical calculation is hard to manage. As a result, SPH has difficulty modeling a Keplerian accretion disk around a compact object or a protostar. The physical shear viscosity plays an important role in the evolution of an accretion disk (Shakura & Sunyaev, 1973; Lynden-Bell & Pringle, 1974; Pringle, 1981; Frank et al., 2002), so accurate handling the physical and numerical viscosities in a numerical code is very important in these simulations.

Several approaches have been tried to address this problem. One of them is the derivation of an empirical viscosity by the modification of the coefficients appearing in the artificial viscosity (Lubow, 1991; Artymowicz & Lubow, 1994). Another method is the direct implementation of the physical viscosity in an (inviscid) SPH code (Flebbe et al., 1994; Takeda et al., 1994; Watkins et al., 1996; Speith & Riffert, 1999; Lanzafame, 2003). The latter showed good results especially in cataclysmic variable (hereafter CV) systems typically consisting of a Roche lobe filling low-mass main sequence donor star transferring mass to a white dwarf primary star by way of an accretion disk.

Another approach to solve the side effect of the artificial viscosity is the use of a Riemann solver instead of the artificial viscosity (Inutsuka, 2002; Cha & Whitworth, 2003). This method is called GodunovSPH (hereafter GSPH), because Godunov (1959) firstly employed the Riemann solver between the adjacent cells in the finite difference method. A one–dimensional Riemann solver has been used on the line joining and particles in GSPH, and provides the time-independent pressure and velocity at the contact discontinuity as a result. The pressure and velocity are used in the momentum and energy equations to update the velocity and the specific internal energy of the particle in the time domain. The Riemann solver generates a small but sufficient dissipation to describe shockwaves. Cha & Whitworth (2003) proved that the Riemann solver can capture high Mach number shocks, and effectively prevents the particle penetration even in a head–on collision of fluids without the addition of any extra numerical dissipation, such as artificial viscosity.

Another problem of SPH is the numerical surface tension (Agertz et al., 2007). The momentum equation of standard SPH shows an unphysical force across the density contrast even in pressure equilibrium, and the unphysical force acts like a surface tension on the contact discontinuity (Cha et al., 2010). The numerical surface tension suppresses the initial perturbation effectively, so any hydrodynamic instability, such as the Kelvin-Helmholtz instability cannot grow in the standard SPH simulation. Price (2008) suggested artificial conduction for the thermal communications of adjacent particles to reduce the unphysical force, and showed that standard SPH employing the artificial conduction can describe the Kelvin-Helmholtz instability across a density contrast.

A different approach to emphasize the consistency of a particle code has been suggested by Inutsuka (2002). He suggests mutually smoothed equations of motion and energy in GSPH. The mutually smoothed equations consider the host and neighbor particles as a smoothed one, while standard SPH or the method of Cha & Whitworth (2003)3 considers the neighbor particles as point masses. The mutually smoothed momentum equation is able to get the –order consistency, and does not show any numerical surface tension (Cha et al., 2010). Eventually, GSPH can describe the Kelvin–Helmholtz instability across the density contrast without any numerical dissipation term of energy such as the artificial conduction. Fuerthermore, GSPH can hold magnetohydrodynamics with riemann solvers of a magnetized fluid (Iwasaki & Inutsuka, 2011; Tsukamoto et al., 2013, 2015)

We have developed a parallel GSPH code with the Message Passing Interface (MPI4). With our parallel GSPH code, we are working to simulate the dynamical evolution of a CV accretion disk. CVs provide a useful testbed for exploring astrophysical viscosity. They are physically compact and hence vary on timescales short enough that hundreds to thousands of orbital cycles of time series data can be collected in a single season from a variety of systems (Wood & Burke, 2007; Wood et al., 2009). The substantial and growing database of observations provide extremely useful reality check for the theoretical calculations. While we are still some years away from fully-physical radiation-magnetohydrodynamical models of accretion disks, advances in computational hardware and software continue to close that gap. Without a doubt the question of the nature of astrophysical viscosity is central to the time evolution of astrophysical accretion disks.

As the first step of this plan, we implemented several numerical routines to handle the shear viscosity stress tensor in the momentum and energy equations of GSPH. Furthermore, the implementation has been tested in the context of the Couette flow (Couette, 1890). By comparison of the numerical results with the analytic solution, we can demonstrate that our implementation works correctly.

In the practical calculation, not only the viscous acceleration, but the hydrodynamics acceleration should be calculated simultaneously. We have performed the full hydrodynamics simulations with the shear viscosity calculation in Couette flow as well. All results show the expected result based on the analytic solution of the genuine Couette flow.

The kinematic viscosity of a fluid is not a constant in general. In order to manage the individual kinematic viscosity, a smoothed dynamic viscosity has been suggested. We have confirmed that the smoothed dynamic viscosity is able to describe the viscous acceleration when the particles have a different kinematic viscosity.

We have performed the simulations of the Couette flow with a density contrast to see the difference between the standard SPH and GSPH. GSPH shows a less dispersive result, but the difference is not meaningful.

The Couette flow is good to check the correctness of the simulations with viscosity. However, we need a more realistic test problem to check our implementation. The evolution of the viscous ring around a central gravity (Lynden-Bell & Pringle, 1974) has been performed as another test. The results also show an agreement to the analytic solution.

Sections 2 and 3 briefly announce the equations of GSPH and the implementation of the shear viscosity, respectively. Sections 4 and 5 contain the simulations of the plane Couette flow and the viscous ring evolution, respectively. Finally, section 6 presents the conclusion and discussion.

2 Fundamentals of GSPH

The Lagrangian hydrodynamics of the inviscid fluid is described by


where is the specific total energy, and the other variables have their usual meanings. Here, the dot above physical quantities means the total derivative (i.e., Lagrangian derivative with respect to time). Instead of the total energy equation, the specific internal energy, has been used more frequently in the particle codes, and the specific internal energy equation becomes


The equation of state,


is needed to close the system. Here, is the specific heat ratio.

With the kernel convolution or the Euler-Lagrange equation, Eqs. (2) and (4) are converted to the momentum and energy equations of GSPH, respectively (Inutsuka, 2002),


where the starred quantities, and are the results of the Riemann solver. Note that we have used lowercase Greek characters (, , …) as the indices of the direction, and and as the particle indices in this paper. We have used the iterative Riemann solver suggested by van Leer (1997) to get the starred quantities. GSPH describes shockwaves without any numerical dissipation because of the Riemann solver, and is free from the side effect of the artificial viscosity. Here, is the updated velocity of the host particle by , and is the integral of the squared particle volume. in Eqs (6) and (7) is the gradient of kernel function, given by


All details of Eqs. (6) and (7) are described in Eqs. (61) – (64) of Inutsuka (2002).

3 Implementation of the physical viscosity

3.1 Viscous stress tensor

The viscous stress tensor, acting on a fluid element is described by


where and are the shear and bulk viscosity stress tensors, given by


component–wise, respectively. Here, is the Kronecker delta. Furthermore, and are the dynamic viscosities of the shear and bulk components, respectively. We will omit the bulk viscosity because the idealized monoatomic gas has been assumed in the simulations.

The kinematic shear viscosity, has a relation with the dynamic shear viscosity ,


3.2 Acceleration and energy dissipation due to the viscosity

The acceleration, due to the viscosity is given by


The viscosity street tensor contains the gradient of velocity already, so the viscous acceleration has the second derivative of velocity components. It is well known that the simple treatment of the second derivative by the standard SPH formulation is very sensitive to the particle disorder, and therefore too noisy. It happens not only in the viscous fluid calculation, but in any kind of simulation including the second derivative, e.g., radiative transfer with the diffusion approximation (Viau et al., 2006). Several previous lines of research have been explored to solve this problem, and most of them can be categorized into three groups. Fatehi et al. (2009) presented a nice review and comparison on this subject. We will summarize Fatehi et al. (2009) below for the convenience of the reader.

The first method is the double summation (Flebbe et al., 1994; Watkins et al., 1996; Speith & Riffert, 1999; Lanzafame, 2003). Perhaps, it is the most widely used method in the physical viscosity implementation of SPH. The second derivative of a physical quantity, is given by




and is a coefficient of the physical process, e.g., viscosity. Eq. (15) uses a symmetric form of particles, but a non–symmetric form,


can be considered as well. We have used he symmetric form (Eq. (15)) all the test calculations below.

In the double summation method, every particle should calculate the derivative with its own neighbors by Eq. (15) firstly. Then, the second derivative of the quantity is calculated by Eq. (14). The essence of the double summation method is to use the information of the “neighbors of a neighbor.” With Eqs. (14) and (15), particle can use not only the information of particle , but also the information of the neighbors of particle . Another good property of the double summation method is that the boundary particles can be managed easily. Obviously, this method has a fatal disadvantage in the calculation speed. However, it is inevitable to consider the information the neighbors of a neighbor. We have used this double summation method in our implementation.

Another widely used method is the difference scheme (Brookshaw, 1985; Cleary & Monaghan, 1999; Monaghan, 2005; Viau et al., 2006). It uses the single summation to calculate the second derivative, so should be more efficient in calculation speed. Furthermore, this method shows a excellent result when particles have a different , and the contrast of them is very steep (Cleary & Monaghan, 1999; Viau et al., 2006). In contrast to the difference scheme, The double summation method has difficulty when the difference between and is large. So we will suggest a smoothed dynamic viscosity, and show that it is effective in the practical calculations in Sec. 4.4.

The last method is to use the second derivative of the kernel function (Chaniotis et al., 2002). We do not consider this method here.

We have employed the double summation method in the viscosity implementation. First of all, the viscous stress tensor has been evaluated by Eq. (10), and then the viscous acceleration is calculated by the GSPH momentum equation,




Here, is the unit vector in -direction.

The viscous coupling between the fluid elements is a typical dissipative mechanism, so, the effect of the viscosity should be considered in the energy equation. The equation of the total specific energy with the viscosity dissipative term becomes


By comparison with Eq. (3), the pure contribution of the viscosity on the internal energy is given by


With the similar way for the momentum equation, the viscous energy equation becomes,




There is no ambiguity in Eqs. (18) and (22) because is a symmetric tensor. A full description on the viscous energy dissipation can be found in Flebbe et al. (1994).

The viscous acceleration and energy dissipation are added to Eqs. (6) and (7) to make the final acceleration and internal energy increment by the operator splitting method, respectively.

4 Plane Couette flow

The plane Couette flow (hereafter Couette flow) (Couette, 1890) describes the motion of a viscous flow between two moving parallel plates. One plate is at the top, and the other is at the bottom of the fluid. The plates at the top is moving with a constant velocity along the direction of the fluid extension, and the bottom one is moving by a different velocity or at rest. With the non-slip boundary condition between the plates and the fluid, the fluid gains momentum from the plate movement due to the viscous coupling. Although it is the simplest shear flow which follows Newton’s law of viscosity, actual experiments have been rarely done because of practical difficulties (Kitoh et al., 2005).

We assume a two–dimensional situation. The fluid is at rest initially, and the plate at the top of the fluid () is moving right with a constant velocity, , while the plate at the fluid bottom () is at rest. The Couette flow under this circumstance is described by


where is the –component of the fluid velocity. Here, the kinematic viscosity has been assumed as a constant all over the fluid. It is the simplest diffusion equation with an inhomogeneous boundary condition,


and the initial velocity profile,


where is the calculation domain size. The analytic solution is to be found easily by the separation of variables, and is given by


4.1 Problem setup

We have used two–dimensional GSPH in the Couette flow simulations, and all physical quantities in the code are dimensionless. A unity box () with grids in the –plane has been used in the calculation. Initially, the particles are uniformly positioned at the center of the grids to minimize the numerical noise.

The moving blocks are treated as the boundary zone in the simulations, so they are beyond the calculation domain. The thickness of the blocks in –direction is 0.2, and is enough to cover the several times of the smoothing length of the fluid particles.

The upper moving block should have a constant velocity, , and the bottom is at rest in the Couette problem. However, the moving blocks are considered as a boundary zone in the simulations, so the upper and lower blocks should have a variable velocity. It is because the the non–slip boundary condition has been employed (i.e., extremely strong viscous coupling) in the tests. The value of is set to in the simulation, but the actual velocity of the moving blocks has been determined by the extrapolation of the analytic solution (Eq. (26)) to the boundary zone at every evolution time.

Another practical problem is the startup velocity of the fluid. A sudden movement of the top block with respect to the fluid may cause a slip, and a noise in the velocity profile appears at the very early stage. The noise won’t decay even in the later evolution. So, we have set an the initial time, ,


and determine a corresponding velocity of the fluid at by Eq. (26) as the startup velocity. Here, is the viscous timescale, given by


For the startup velocity of the fluid, one can use


instead of Eq. (26). Here, erfc() is the complementary error function. Eq. (29) is more convenient in the simulations.

Two cases of simulations have been performed. One is the genuine Couette flow. The fluid particles are moved by the viscosity acceleration only in this case. Therefore, any hydrodynamic acceleration or energy calculation have not been involved in this case. However, note that neighbor finding, the density update and the individual smoothing length determination based on the updated density have been calculated. So, one can see the numerical dispersion in this case if it exists.

The other case simulated here is the full hydrodynamics calculation with the viscosity. In this case, hydrodynamics acceleration and energy calculation by Eqs. (6) and (7) have been added to the viscous acceleration and energy dissipation. There is no analytical description in this case, however, one may expect the almost same result of the former case. It is because the pressure equilibrium has been assumed, and the main driving force of the fluid motion is the viscous coupling with the moving plates rather than hydrodynamics. We expect that the latter may be useful to check the influence of the numerical dispersion in our GSPH code.

The density of the fluid is set to the unity initially in the whole calculation domain, and the pressure is set to


where, is the specific heat ratio, and is set to in the simulations with hydrodynamics. This pressure value makes the sound speed unity.

The calculations with hydrodynamics have been performed with the second-order GSPH code. Refer to Inutsuka (2002) and van Leer (1997) to see the details of the second-order scheme.

4.2 Results of the Couette flow

We have performed several simulations with various kinematic viscosities, from to , and all simulations show the same result at each normalized evolution time by the viscosity timescale.

Figure 1: Simulations of the Couette flow with and in left and right panels, respectively. These are the velocity–position plots of the viscosity–only calculations. The red solid lines are the analytic solution given by Eq. (26) at from the left to the right, and the black dots are the numerical solution at the same time. Around , the fluid nearly achieved the equilibrium solution.

Fig. 1 shows the position-velocity plots of the simulations of and cases without hydrodynamics. The red solid lines are the analytic solution calculated by Eq. (26), and the black dots are the numerical results. The numerical and analytic solutions show a good agreement, and the equilibrium solution,


has been almost obtained around . These are two–dimensional simulations, so one black dot in the figures represents a group of fluid particles rather than a single fluid particle. One can see the suppressed numerical dispersion. Note that these simulations calculate the viscosity force only, but the other GSPH routines, e.g., finding neighbors, the density estimation and the smoothing length determination are used in the calculation.

The viscosity force calculation by Eq. (21) has been confirmed successfully by the tests above, but we have to check the overall performance of our viscous GSPH code. It is because the hydrodynamic acceleration and adiabatic changes of the internal energy are also important in the viscous fluid motion. So, we have performed the same simulations again with a full steps of hydrodynamics. Fig. 2 is the results with hydrodynamics. The kinematic viscosity and other configurations of the simulations are the same as the viscosity–only simulations.

Figure 2: This is the same simulations as Fig. 1, but hydrodynamics included. The kinematic viscosity and other configurations of the simulations of left and right panels are the same as Fig. 1. The contribution of the pressure gradient are negligible, so the results are nearly the same as the simulations without hydrodynamics. One can see that the numerical dispersion is suppressed.

Although there is no analytic solution in this case, one may expect nearly the same results as the viscosity–only calculation, because the contribution of the pressure gradient is negligible. Essentially, there is no significant difference in Figs. 1 and 2, as expected.

4.3 High Reynolds number simulations

The Reynolds number defined by


has been frequently used as a parameter of the possibility of turbulence. The simulations in the previous section have the kinematic viscosity range in , and the corresponding values are in the range . No dispersion appears in this range.

We have performed some simulations with a higher number to define the possibility of the occurrence of turbulence. The configuration of the simulations has been changed in these tests. Specifically, the two blocks at the top and bottom of the fluid are moving in opposite directions with the same speed. There is an analytic solution to this boundary condition, given by


The system size in Eq. (33) is set to , and the plate moving velocity, is set to .

Figure 3: The simulations with a higher value. The value of is set to , and the corresponding kinematic viscosity is . The left and right panels are without and with hydrodynamics calculations, respectively. Note that the output times of the left and right panels are different. It is because the viscosity–only calculation with a high value becomes unstable in the later time. In the left panel, , and the right panel, . They are more dispersive than the smaller cases.

Fig. 3 shows the results of the high simulations. The left and right panels show the results without and with hydrodynamics, respectively. Note that the output times are different in the two panels because the result of the viscosity–only calculation becomes very unstable in the later times. The kinematic viscosity used in the simulations is , and the corresponding value becomes 100. The fluid becomes turbulent around . Although it is still controversial, the Couette flow is expected to be turbulent around (Schneider et al., 2010). Our results show the occurrence of turbulence at a smaller value more or less. Another interesting point is that the full hydrodynamics calculation is more stable (i.e., less dispersive in the final snapshot) than the pure viscosity calculation.

Fig. 4 shows the velocity dispersion along – and –directions. The red big and black dots are the results of and , respectively. The both are the results of the full hydrodynamics calculations. First of all, the low case show very clean profiles in –direction without any dispersion. The curved features of the low case in along –direction is due the imperfect settling to the equilibrium solution. Fig. 4 is at , and one may expect a less curved pattern around . One can see a very dispersive pattern in the high simulation. The viscosity coupling between the moving blocks and the fluid acts on in –direction only, but all velocity components in every direction gain an acceleration, and the fluid shows the turbulent pattern eventually.

Figure 4: The dispersion of velocities in – and –directions. The red thick and black thin dots are the results of and cases, respectively. The corresponding values are and , respectively. More dispersive pattern arises in the high case, and the severe dispersion may cause a turbulence. In the low case, the curved feature of in –direction is due the the imperfect settling to the equilibrium solution at .

4.4 A smoothed dynamic viscosity : Implementation and test

We have tested our viscosity managing routines with a universal kinematic viscosity for all particles so far. However, the kinematic viscosity of a particle is not constant over the calculation domain in general, and is a function of the physical quantities of the fluid (e.g., density and temperature). We have found that Eqs. (17), (21) are not effective when the the host and neighbor particles have a different kinematic viscosity. A smoothed dynamic viscosity has been introduced to solve this problem. The smoothed dynamic viscosity is determined by


where is approximated by


Note that the kinematic viscosity of the individual particle is not a smoothed quantity, and should be determined in advance by the physical conditions of the particle.

The same momentum and energy equations, Eqs. (17) and (21) are used, but the dynamic viscosities, and of the equations should be replaced by the smoothed dynamic viscosities, and .

Two layers of a viscous flow have been considered as a test problem. The fluid is assumed to be extended infinitely in the – and –directions, and has different kinematic viscosities in the upper and lower layers. Contrary to the genuine Couette flow, there is no moving solid block in this test, but the upper layer is moving by a initial velocity, to the right, while the lower layer is at rest initially. The viscous coupling between the two layers will generate an acceleration in the layers. The kinematic viscosity is a function of –poisition only, so they are given by


There is an analytic solution for this infinitely extended viscous flow (Carslaw & Jaeger, 1959), and the –component of the fluid velocity is give by


for the upper layer, and


for the lower layer. Here, , and erf is the error function. The value of is set to 0.1 in the test.

Figure 5: Two layers simulation results at . The steeper profile is the later result. The left and right panels are without and with hydrodynamics, respectively. The red solid lines and black dots are the analytic and numerical solutions, respectively. The two solutions agrees with each other well. The kinematic viscosity of the upper and lower layers are 0.1 and 0.2, respectively.

Fig. 5 is the results of the two–layered viscous flow. The kinematic viscosities of the upper and lower layers are and , respectively. The output times of the figure are scaled by the viscous time scale, defined by


The numerical solution (black dots) coincides with the analytic solution (red solid lines) very well, and numerical dispersion does not appear at all.

Figure 6: Two layers simulations of a higher contrast in the kinematic viscosity. The kinematic viscosities of the upper and lower layers are 0.5 and 0.1, respectively. The agreement between the numerical and analytic solutions are still good. The output times are the same as Fig. 5.

Fig. 6 is the results of the higher contrast in the kinematic viscosities. The kinematic viscosity of the upper layer is 5 times bigger than the lower layer. Still, the numerical and analytic solutions show a good agreement.

4.5 Test with a density contrast : Comparison with the standard SPH

Perhaps one interesting subject is the comparison between the standard SPH and GSPH. However, the Couette flow should not be a good example to show the difference of the two numerical methods, because there is no hydrodynamics in the Couette flow. Furthermore, it is hard to find a significant difference in the two methods even in the full hydrodynamics calculation. It is obvious because the Couette flow assumes the pressure equilibrium, so the hydrodynamic force is negligible compare to the viscous one as presented in Fig. 2.

Another consideration is necessary to see the difference. So we have performed the hydrodynamic Couette flow tests with a density contrast. As we have described in Sec. 1, the standard SPH shows the numerical surface tension across the density contrast, so one may expect to see a difference in the results of the standard SPH and GSPH.

Essentially, all configurations of the simulation are the same as the two layers simulations in the previous section except the density contrast. The lower layer is denser and more viscous than the upper one by 4 and 2 times, and the density and viscosity of the lower layer are set to 4 and 0.2, respectively. Note that the contrast of the dynamic viscosity becomes 8 across the two layers.

Initially, the upper layer moves to the right by the velocity of 0.1. Without the hydrodynamics, the analytic solution given in Eqs. (37) - (38) can be used by the modified and . The evolution time of the simulations are scaled by the viscous time scale, given in Eq. (39).

Figure 7: Simulations of the Couette flow with a density contrast. Each panel shows the result at from the upper left to the lower right. The left- and right-hand-side are the results of GSPH and the standard SPH, respectively. The embedded small boxes are the closeup at the boundary of the two different density layers. The red solid lines are the analytic solution given in Eqs. (37) - (38), and the solid circles and crosses are the upper and lower fluid, respectively. No significant difference appears, but the standard SPH simulation shows a more dispersive pattern at the boundary of the layers. The lower layer is denser and more viscous, and the density and kinematic viscosity of the lower layer are set to 4 and 0.2 in the code unit, respectively.

Fig. 7 shows the results. The left- and right-hand-side are the results of GSPH and the standard SPH, respectively. In the SPH simulation, the coefficients of the artificial viscosity (they are called in general) are set to 1 and 2, respectively. Any other kind of treatment, e.g., the artificial conduction (Price, 2008) or Balsala switch (Balsara, 1995) has not been used. The red solid lines are the analytic solution without hydrodynamics given by Eqs. (37) - (38). The solid circles and crosses are the upper and lower layer material, respectively. Apparently, there is no significant difference, but the numerical dispersion of GSPH looks less serious, especially, at the border of the two layers in the first and last snapshots. It may be understandable, because the two layers are parallel, so the curvature becomes 0. The numerical surface tension is very similar to the physical one, so the higher curvature will make the stronger surface tension (Cha et al., 2010). One may see a clearer difference in the results, if any other curved geometry (e.g. cylindrical fluid) has been employed.

5 Viscous ring evolution

The Couette flow is a good test problem for the viscous fluid simulations due to its simplicity and clear analytic description. However, we need a more realistic test problem to check our implementation of the physical viscosity. Lynden-Bell & Pringle (1974) derived the self-similar solution for the dynamical evolution of a viscous gas ring under a central gravity, and it becomes a good test for a numerical code of the shear viscosity (e.g. Flebbe et al., 1994). Only a brief description of the problem is presented here, because there are many good sources that explain this problem (e.g. Pringle, 1981; Frank et al., 2002).

A ring of gas is located around the central gravitational object. The total mass and the initial radius of the ring are set to and , respectively. The ring rotates around the central object, and spreads due to the kinematic viscous, . A dimensionless time,


and radius,


are useful to describe the evolution of the ring. The time evolution of the surface density of the ring is given by


where and are the modified Bessel function and , respectively. Furthermore, the radial velocity of the ring is given by


The initial distribution of the ring material is set to the delta function,


in the analytic work of Lynden-Bell & Pringle (1974), but we choose the solution at as the initial condition. Note that the radial velocity of our initial condition should be defined by Eq. (43) rather than . Initially, the mass and radius of the ring are set to 1, and the kinematic viscosity, is set to or . The two–dimensional GSPH code has been used, and the total number of particles in the simulations is 10410.

Figure 8: The surface density of viscous ring evolution test. The left and right panels are the result with and , respectively, The evolution times are at . The red solid lines are the analytic solution, and the dots are the surface density of the ring. The overall agreement are good, but the evolution is little bit faster than the analytic solution. Furthermore, the smaller kinematic viscosity shows the more dispersive pattern.

Fig. 8 shows the results. The left and right panels are and , respectively. The red solid lines are the analytic solution at , and the black dots are the surface density of the numerical calculations. Note that the snapshot at is the initial condition. The results shows an agreement to the analytic solution, but a dispersive pattern happens in the late evolution stage of the small kinematic viscosity case (Speith & Riffert, 1999; Iwasaki, 2015; Sugiura & Inutsuka, 2016). Furthermore, the numerical results both are faster more or less than the analytic solution (Watkins et al., 1996).

6 Conlusion and discussion

Although there are many advantages of SPH as a multidimensional Lagrangian code, the side effect of the artificial viscosity is a great obstacle to simulate the accretion disks in CVs and around a protostar. The exact contribution of the physical shear viscosity in the standard SPH is hardly managed due to the interference of the artificial viscosity.

We have implemented several numerical routines to manage the viscous acceleration and energy dissipation in our GSPH code. Contrary to the standard SPH, GSPH uses the Riemann solver to describe the fluid motion and energy dissipation, therefore, it is free from the side effect of the artificial viscosity. The double summation method has been used in the treatment of the viscosity stress tensor and its derivatives. The second derivative of the velocity can be described successfully, and the numerical noise suppressed by this method.

The plane Couette flow has been chosen as a test problem. Our implementation shows a good agreement with the analytic solution. We have performed not only the viscosity–only simulations, but also the full hydrodynamics simulations of the viscous fluid. The hydrodynamics simulations also show an expected results.

Simulations with a high value have been performed as well. We can see a dispersive pattern around , and it is thought to be an occurrence of the turbulence. An experimental and analytic expectation of value for the turbulence occurrence in the Couette flow is approximately, so our simulation looks more sensitive against the occurrence of the turbulence.

In general, the kinematic viscosity is not a constant all over the calculation domain, and is a function of the physical properties of a fluid element. The smoothed dynamic viscosity has been introduced to describe the individual viscosity of particles. Simulations with two layers which have different kinematic viscosities have been performed, and the results show a good agreement to the analytic solution, even in a high contrast of the kinetic viscosity.

The same two layers simulation has been formed with a density contrast as well. The result of GSPH shows a less dispersive pattern at the boundary of the different density layers, but the results of the standard SPH and GSPH does not show any meaningful difference. It is because the numerical surface tension is not very critical due to the plane geometry.

The viscous gas ring simulation has been performed. The numerical solution can follow the analytic solution well, but one can see a dispersive pattern at the later stage of the small viscosity case. Furthermore, the evolution shows a faster evolution than the analytic solution more or less.

We do not claim that GSPH is a complete inviscid code. The Riemann solver generates a non–zero dissipation to describe shockwaves. However, it is hard to estimate how much the dissipation is in the Godunov method (Dunhill et al., 2013; Puri & Ramachandran, 2014). Perhaps, the intrinsic numerical viscosity of GSPH should be added to the physical viscosity somehow. However, we do not expect the mixing between the numerical and physical viscosities to be a simple arithmetic sum. We guess it may be the reason why the critical value between the laminar and turbulent flows in our simulations is different to the previous study. The exact contribution of the intrinsic viscosity of GSPH, and a better understanding for the occurrence of the turbulence are left for a future study.


We thank Anthony Whitworth and Andrew King for useful discussion on the selection of the test problem. This material is based upon work supported by the National Science Foundation under Grant No. AST-1305799 to Texas A&M University-Commerce and by NASA under grant 11-KEPLER11-0038.


  1. pagerange: GodunovSPH with shear viscosity : implementation and testsGodunovSPH with shear viscosity : implementation and tests
  2. pubyear: 2015
  3. There is an important difference between Inutsuka (2002) and Cha & Whitworth (2003). Inutsuka (2002) uses the mutually smoothed motion and energy equations, while Cha & Whitworth (2003) uses the same equations as the standard SPH. Specifically, Inutsuka (2002) considers both of the host and neighbor particles as smoothed, while Cha & Whitworth (2003) considers only the host particle as smoothed. The neighbor particles in the standard SPH and Cha & Whitworth (2003) are considered as point masses. Therefore, we chose that GSPH uses the method of Inutsuka (2002). In order to distinguish the two methods, Cha & Whitworth (2003) is called GPH sometimes. If the kernel function of neighbor particles is assumed to be a delta function, then GSPH is reduced to GPH.
  4. For the standard documents of MPI, refer to


  1. Agertz O., Moore B., Stadel J., Potter D., Miniati F., Read J., Mayer L., Gawryszczak A., Kravtsov A., Nordlund Å., Pearce F., Quilis V., Rudd D., Springel V., Stone J., Tasker E., Teyssier R., Wadsley J., Walder R., 2007, MNRAS, 380, 963
  2. Artymowicz P., Lubow S. H., 1994, ApJ, 421, 651
  3. Balsara D. S., 1995, J. Comp. Phys., 121, 357
  4. Barnes J., Hut P., 1986, Nature, 324, 446
  5. Brookshaw L., 1985, PASA, 6, 207
  6. Carslaw H. S., Jaeger J. C., 1959, Conduction of heat in solids. Oxford University Press
  7. Cha S.-H., Inutsuka S.-I., Nayakshin S., 2010, MNRAS, 403, 1165
  8. Cha S.-H., Whitworth A. P., 2003, MNRAS, 340, 73
  9. Chaniotis A. K., Poulikakos D., Koumoutsakos P., 2002, J. Comp. Phys., 182, 67
  10. Cleary P. W., Monaghan J. J., 1999, J. Comp. Phys., 148, 227
  11. Couette M., 1890, Ann. de Chimie et de Physique, 21, 433
  12. Dunhill A. C., Alexander R. D., Armitage P. J., 2013, MNRAS, 428, 3072
  13. Fatehi R., Fayazbakhsh M., Manzari M. T., 2009, Int. Jour. Nat. and App. Sci, 3, 50
  14. Flebbe O., Muenzel S., Herold H., Riffert H., Ruder H., 1994, ApJ, 431, 754
  15. Frank J., King A., Raine D. J., 2002, Accretion Power in Astrophysics: Third Edition. Cambridge University Press
  16. Gingold R. A., Monaghan J. J., 1977, MNRAS, 181, 375
  17. Godunov S. K., 1959, Math. Sbornik, 47, 271
  18. Inutsuka S.-I., 2002, J. Comp. Phys., 179, 238
  19. Iwasaki K., 2015, J. Comp. Phys., 302, 359
  20. Iwasaki K., Inutsuka S.-I., 2011, MNRAS, 418, 1668
  21. Kitoh O., Nakabyashi K., Nishimura F., 2005, Jour. Fluid Mech., 539, 199
  22. Lanzafame G., 2003, A&A, 403, 593
  23. Lubow S. H., 1991, ApJ, 381, 268
  24. Lucy L. B., 1977, AJ, 82, 1013
  25. Lynden-Bell D., Pringle J. E., 1974, MNRAS, 168, 603
  26. Monaghan J. J., 1992, ARA&A, 30, 543
  27. Monaghan J. J., 2005, Rep. Prog. Phys., 68, 1703
  28. Price D. J., 2008, J. Comp. Phys., 227, 10040
  29. Pringle J. E., 1981, ARA&A, 19, 137
  30. Puri K., Ramachandran P., 2014, J. Comp. Phys., 270, 432
  31. Schneider T. M., de Lillo F., Buehrle J., Eckhardt B., Dörnemann T., Dörnemann K., Freisleben B., 2010, Phys. Rev. E, 81, 015301
  32. Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337
  33. Speith R., Riffert H., 1999, J. Comp. Appl. Math., 109, 231
  34. Sugiura K., Inutsuka S.-i., 2016, J. Comp. Phys., 308, 171
  35. Takeda H., Miyama S. M., Sekiya M., 1994, Prog. Theor. Phys, 92, 939
  36. Tsukamoto Y., Iwasaki K., Inutsuka S.-i., 2013, MNRAS, 434, 2593
  37. Tsukamoto Y., Iwasaki K., Okuzumi S., Machida M. N., Inutsuka S., 2015, MNRAS, 452, 278
  38. van Leer B., 1997, J. Comp. Phys., 135, 229
  39. Viau S., Bastien P., Cha S.-H., 2006, ApJ, 639, 559
  40. Watkins S. J., Bhattal A. S., Francis N., Turner J. A., Whitworth A. P., 1996, A&AS, 119, 177
  41. Wood M. A., Burke C. J., 2007, ApJ, 661, 1042
  42. Wood M. A., Thomas D. M., Simpson J. C., 2009, MNRAS, 398, 2110
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Comments 0
Request comment
The feedback must be of minumum 40 characters
Add comment
Loading ...

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description