tnew time step formula - sphere - GPU-based 3D discrete element method algorithm with optional fluid coupling
 (HTM) git clone git://src.adamsgaard.dk/sphere
 (DIR) Log
 (DIR) Files
 (DIR) Refs
 (DIR) LICENSE
       ---
 (DIR) commit f6c8394ac2e8e2c86191521693237b39387e8e29
 (DIR) parent c187341cc547490c3de9c3397af82b2893e0d1d4
 (HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
       Date:   Thu, 19 Jun 2014 14:41:49 +0200
       
       new time step formula
       
       Diffstat:
         M python/sphere.py                    |      31 ++++++++++++++++++++-----------
       
       1 file changed, 20 insertions(+), 11 deletions(-)
       ---
 (DIR) diff --git a/python/sphere.py b/python/sphere.py
       t@@ -2312,7 +2312,8 @@ class sim:
                    current = 0.0,
                    file_dt = 0.05,
                    step_count = 0,
       -            dt = -1):
       +            dt = -1,
       +            epsilon = 0.01):
                '''
                Set temporal parameters for the simulation. *Important*: Particle radii,
                physical parameters, and the optional fluid grid need to be set prior to
       t@@ -2332,24 +2333,32 @@ class sim:
                :type step_count: int
                :param dt: The computational time step length [s]
                :type total: float
       +        :param epsilon: Time step multiplier (default = 0.01)
       +        :type epsilon: float
                '''
        
       -        # Computational time step (O'Sullivan et al, 2003)
       -        #self.time_dt[0] = 0.17 * \
       -                #math.sqrt((4.0/3.0 * math.pi * r_min**3 * self.rho[0]) \
       -                #/ numpy.amax([self.k_n[:], self.k_t[:]]) )
       -        # Computational time step (Zhang and Campbell, 1992)
                if dt > 0:
                    self.time_dt[0] = dt
                    if (self.np[0] > 0):
                        print("Warning: Manually specifying the time step length when "
                        + "simulating particles may produce instabilities.")
                else:
       -            r_min = numpy.amin(self.radius)
       -            self.time_dt[0] = 0.075 *\
       -                    math.sqrt((V_sphere(r_min) * self.rho[0]) \
       -                    / numpy.amax([self.k_n[:], self.k_t[:]]) )
       -
       +            r_min = numpy.min(self.radius)
       +
       +            # Radjaii et al 2011
       +            m_min = self.rho * 4.0/3.0*numpy.pi*r_min**3
       +            k_max = numpy.max([self.k_n[:], self.k_t[:]]))
       +            self.time_dt[0] = epsilon/(numpy.sqrt(k_max/m_min))
       +
       +            # Zhang and Campbell, 1992
       +            #self.time_dt[0] = 0.075 *\
       +                    #math.sqrt((V_sphere(r_min) * self.rho[0]) \
       +                    #/ numpy.amax([self.k_n[:], self.k_t[:]]) )
       +
       +            # Computational time step (O'Sullivan et al, 2003)
       +            #self.time_dt[0] = 0.17 * \
       +                    #math.sqrt((4.0/3.0 * math.pi * r_min**3 * self.rho[0]) \
       +                    #/ numpy.amax([self.k_n[:], self.k_t[:]]) )
        
                # Check numerical stability of the fluid phase, by criteria derived by
                # von Neumann stability analysis of the diffusion and advection terms