thalfshear-darcy-stress-mod-starter.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
       ---
       thalfshear-darcy-stress-mod-starter.py (2745B)
       ---
            1 #!/usr/bin/env python
            2 import sphere
            3 import numpy
            4 import sys
            5 
            6 # launch with:
            7 # $ ipython halfshear-darcy-stress-starter.py <device> <fluid> <c_phi> <k_c>
            8 #     <sigma_0> <mu> <shear_stress> <mod_A> <mod_f> <timefactor>
            9 
           10 device = int(sys.argv[1])
           11 wet = int(sys.argv[2])
           12 c_phi = float(sys.argv[3])
           13 k_c = float(sys.argv[4])
           14 sigma0 = float(sys.argv[5])
           15 mu = float(sys.argv[6])
           16 shear_stress = float(sys.argv[7])
           17 mod_A = float(sys.argv[8])
           18 mod_f = float(sys.argv[9])
           19 timefac = float(sys.argv[10])
           20 
           21 if wet == 1:
           22     fluid = True
           23 else:
           24     fluid = False
           25     
           26 sim = sphere.sim('halfshear-sigma0=' + str(sigma0), fluid=False)
           27 print('Input: ' + sim.sid)
           28 sim.readlast()
           29 
           30 sim.fluid = fluid
           31 if fluid:
           32     sim.id('halfshear-darcy-sigma0=' + str(sigma0) + '-k_c=' + str(k_c) +
           33             '-mu=' + str(mu) + '-ss=' + str(shear_stress) + '-A=' + str(mod_A)
           34             + '-f=' + str(mod_f))
           35 else:
           36     sim.id('halfshear-sigma0=' + str(sigma0) + '-ss=' + str(shear_stress) + \
           37             '-shear')
           38 
           39 sim.checkerboardColors(nx=6,ny=3,nz=6)
           40 sim.cleanup()
           41 sim.adjustUpperWall()
           42 sim.zeroKinematics()
           43 
           44 #sim.shear(0.0/20.0)
           45 #sim.shear(1.0/20.0 * velfac)
           46 sim.shear(shear_stress = shear_stress)
           47 K_q_real = 36.4e9
           48 K_w_real =  2.2e9
           49 K_q_sim  = 1.16e9
           50 K_w_sim  = K_w_real/K_q_real * K_q_sim
           51 
           52 if fluid:
           53     #sim.num[2] *= 2
           54     sim.num[:] /= 2
           55     #sim.L[2] *= 2.0
           56     #sim.initFluid(mu = 1.787e-6, p = 600.0e3, cfd_solver = 1)
           57     sim.initFluid(mu = mu*timefac, p = 0.0, cfd_solver = 1)
           58     sim.setFluidBottomNoFlow()
           59     sim.setFluidTopFixedPressure()
           60     #sim.setDEMstepsPerCFDstep(10)
           61     sim.setMaxIterations(2e5)
           62     sim.setPermeabilityPrefactor(k_c)
           63     sim.setFluidCompressibility(1.0/K_w_sim)
           64     # the fluid modulation should be 180 degree out of phase with the wall
           65     # modulation
           66     sim.setFluidPressureModulation(A=mod_A, f=mod_f/timefac, phi=numpy.pi)
           67 
           68 sim.w_sigma0[0] = sigma0
           69 sim.w_m[0] = numpy.abs(sigma0*sim.L[0]*sim.L[1]/sim.g[2])
           70 
           71 #sim.setStiffnessNormal(36.4e9 * 0.1 / 2.0)
           72 #sim.setStiffnessTangential(36.4e9/3.0 * 0.1 / 2.0)
           73 sim.setStiffnessNormal(K_q_sim)
           74 sim.setStiffnessTangential(K_q_sim)
           75 sim.mu_s[0] = 0.5
           76 sim.mu_d[0] = 0.5
           77 sim.setDampingNormal(0.0)
           78 sim.setDampingTangential(0.0)
           79 #sim.deleteAllParticles()
           80 #sim.fixvel[:] = -1.0
           81 
           82 sim.initTemporal(total = 20.0*timefac, file_dt = 0.01, epsilon=0.07)
           83 #sim.initTemporal(total = 20.0, file_dt = 0.00001, epsilon=0.07)
           84 #sim.time_dt[0] *= 1.0e-2
           85 #sim.initTemporal(total = 1.0e-4, file_dt = 1.0e-5, epsilon=0.07)
           86 sim.setTopWallNormalStressModulation(A=mod_A, f=mod_f/timefac)
           87 
           88 # Fix lowermost particles
           89 #dz = sim.L[2]/sim.num[2]
           90 #I = numpy.nonzero(sim.x[:,2] < 1.5*dz)
           91 #sim.fixvel[I] = 1
           92 
           93 sim.run(dry=True)
           94 sim.run(device=device)
           95 sim.writeVTKall()
           96 #sim.visualize('walls')
           97 #sim.visualize('fluid-pressure')