tadaptive-grid-shear-test.py - 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
       ---
       tadaptive-grid-shear-test.py (3992B)
       ---
            1 #!/usr/bin/env python
            2 
            3 # Import sphere functionality
            4 import sphere
            5 
            6 ### EXPERIMENT SETUP ###
            7 initialization = False
            8 consolidation  = True
            9 shearing       = True
           10 fluid          = True
           11 rendering      = False
           12 plots          = True
           13 
           14 # Number of particles
           15 np = 1e4
           16 
           17 # Common simulation id
           18 sim_id = "adapt-grid"
           19 
           20 # Deviatoric stress [Pa]
           21 devslist = [80e3]
           22 #devs = 0
           23 
           24 ### INITIALIZATION ###
           25 
           26 # New class
           27 init = sphere.sim(np = np, nd = 3, nw = 0, sid = sim_id + "-init")
           28 
           29 # Save radii
           30 init.generateRadii(mean = 0.02)
           31 
           32 # Use default params
           33 init.defaultParams(gamma_n = 100.0, mu_s = 0.6, mu_d = 0.6)
           34 init.setYoungsModulus(7.0e9)
           35 
           36 # Add gravity
           37 init.g[2] = -9.81
           38 
           39 # Periodic x and y boundaries
           40 init.periodicBoundariesXY()
           41 
           42 # Initialize positions in random grid (also sets world size)
           43 hcells = np**(1.0/3.0)
           44 init.initRandomGridPos(gridnum = [hcells, hcells, 1e9])
           45 
           46 # Set duration of simulation
           47 init.initTemporal(total = 5.0)
           48 
           49 if (initialization == True):
           50 
           51     # Run sphere
           52     init.run(dry = True)
           53     init.run()
           54 
           55     if (plots == True):
           56         # Make a graph of energies
           57         init.visualize('energy')
           58 
           59     init.writeVTKall()
           60 
           61     if (rendering == True):
           62         # Render images with raytracer
           63         init.render(method = "angvel", max_val = 0.3, verbose = False)
           64 
           65 
           66 
           67 # For each normal stress, consolidate and subsequently shear the material
           68 for devs in devslist:
           69 
           70     ### CONSOLIDATION ###
           71 
           72     # New class
           73     cons = sphere.sim(np = init.np, nw = 1, sid = sim_id + "-cons-devs{}".format(devs))
           74 
           75     # Read last output file of initialization step
           76     lastf = init.status()
           77     cons.readbin("../output/" + sim_id + "-init.output{:0=5}.bin".format(lastf), verbose=False)
           78 
           79     # Periodic x and y boundaries
           80     cons.periodicBoundariesXY()
           81 
           82     if fluid:
           83         # set fluid and solver properties
           84         cons.initFluid(mu=2.080e-7, p=0.0, cfd_solver=1)
           85         cons.setFluidTopFixedPressure()
           86         cons.setFluidBottomNoFlow()
           87         cons.setMaxIterations(2e5)
           88         cons.setPermeabilityPrefactor(2.0e-16)
           89         cons.setFluidCompressibility(1/2.2e9)
           90 
           91     # Setup consolidation experiment
           92     cons.consolidate(normal_stress = devs)
           93     cons.adaptiveGrid()
           94 
           95     # Set duration of simulation
           96     cons.initTemporal(total = 1.5)
           97 
           98     if (consolidation == True):
           99 
          100         # Run sphere
          101         cons.run(dry = True) # show values, don't run
          102         cons.run() # run
          103 
          104         if (plots == True):
          105             # Make a graph of energies
          106             cons.visualize('energy')
          107             cons.visualize('walls')
          108 
          109         cons.writeVTKall()
          110 
          111         if (rendering == True):
          112             # Render images with raytracer
          113             cons.render(method = "pres", max_val = 2.0*devs, verbose = False)
          114 
          115 
          116     ### SHEARING ###
          117 
          118     # New class
          119     shear = sphere.sim(np = cons.np, nw = cons.nw, sid = sim_id + "-shear-devs{}".format(devs))
          120 
          121     # Read last output file of initialization step
          122     lastf = cons.status()
          123     shear.readbin("../output/" + sim_id + "-cons-devs{}.output{:0=5}.bin".format(devs, lastf), verbose = False)
          124 
          125     # Periodic x and y boundaries
          126     shear.periodicBoundariesXY()
          127 
          128     if fluid:
          129         # set fluid and solver properties
          130         shear.initFluid(mu=2.080e-7, p=0.0, cfd_solver=1)
          131         shear.setFluidTopFixedPressure()
          132         shear.setFluidBottomNoFlow()
          133         shear.setMaxIterations(2e5)
          134         shear.setPermeabilityPrefactor(2.0e-16)
          135         shear.setFluidCompressibility(1/2.2e9)
          136 
          137     # Setup shear experiment
          138     shear.shear(shear_strain_rate = 0.05)
          139     shear.adaptiveGrid()
          140 
          141     # Set duration of simulation
          142     shear.initTemporal(total = 20.0)
          143 
          144     if (shearing == True):
          145 
          146         # Run sphere
          147         shear.run(dry = True)
          148         shear.run()
          149 
          150         if (plots == True):
          151             # Make a graph of energies
          152             shear.visualize('energy')
          153             shear.visualize('shear')
          154 
          155         shear.writeVTKall()
          156 
          157         if (rendering == True):
          158             # Render images with raytracer
          159             shear.render(method = "pres", max_val = 2.0*devs, verbose = False)