Difficulty using esys-particle for large-scale thrust fault model

Hi All,

I'm trying to use ESyS-Particle to model a large-scale thrust fault on the planet Mercury. I've set up a 2D simulation, with particles generated by gengeo. The model is 80 km wide and 10 km high. The fault is simulated as an inclined plane (or line in 2D) of broken bonds in the subsurface and is blind (i.e. it doesn't extend all the way to the surface of the model).

I'm using 3 scripts in total. The first one simply generates the particle assembly and writes it to a .geo file. The second script generates bonds between the particles, breaks the bonds along a line segment (the fault), tags the particles (for visualization) and writes it to another .geo file. The final script, which is shown below, loads the particle assembly with bonds, adds walls to the model, and moves the left wall incrementally to the right, in order to simulate compression.

I can make the model work using unit density and RotBondPrms with unrealistic bogus parameters and I can get it to display 'realistic' bulk rock behavior. However, when using the correct density and BrittleBeamPrms with realistic parameters for young's modulus, cohesion etc., the model has a tendency to become very unstable. Even when I decrease the time step to very small values (10e-9), I can't reproduce realistic bulk rock behavior.

I'm using the following base units:

1 length unit = 400 m
1 mass unig = 1 kg
1 time unit = 1 s

Could someone advise if I'm using ESys-Particle in the correct way? Can it be used for something large-scale like this in the first place? The third script, which executes the deformation, is pasted below. I would really appreciate any help you might be able to offer, I've been struggling with this for a long time.

Many thanks!
Rubio Vaughan

begin script:
--------------------------------------------------------------------------------
#thrust.py: Fault propagation fold simulation using ESyS-Particle
# Author: R. Vaughan
# Date: 11 January 2014
# Organisation: VU University Amsterdam
# (C) All rights reserved, 2014.
#
#
#import the appropriate ESyS-Particle modules:
from esys.lsm import *
from esys.lsm.util import *
from esys.lsm.geometry import *
from WallLoader import WallLoaderRunnable

scale = 400 # 1 model unit is 400 meters
maxRadius = 400 / scale
minRadius = 100 / scale
width = 80000 / scale # 80 km
height = 10000 / scale # 10 km
timeSteps = 100000
timeStepIncr = 1000
timeStepSize = 1.0000e-06

particleTag1 = 1
particleTag2 = 2
particleTagFault = 4

#instantiate a simulation object:
sim = LsmMpi (numWorkerProcesses = 1, mpiDimList = [1,1,1])

#initialise the neighbour search algorithm:
sim.initNeighbourSearch (
 particleType = "RotSphere",
 gridSpacing = 2.5 * maxRadius,
 verletDist = 0.2 * minRadius
)
sim.force2dComputations(True)

#set the number of timesteps and timestep increment:
sim.setNumTimeSteps (timeSteps)
sim.setTimeStepSize (timeStepSize)

#specify the spatial domain for the simulation
margin = height / 2
domain = BoundingBox(
 Vec3(-margin,-margin,0),
 Vec3(width+margin,height*3,0)
)
sim.setSpatialDomain (domain)

sim.readGeometry("geometry2.geo")

sim.setParticleDensity (
 tag = 1,
 mask = 0,
 Density = (scale^3) * 2800 / 0.83 # based on 2800 kg/m3, adjusted for scale and porosity
)

#create a wall at the bottom of the model:
sim.createWall (
 name = "bottom_wall",
 posn = Vec3(0.0, 0.0, 0.0),
 normal = Vec3(0.0, 1.0, 0.0)
)

#create a wall at the left side of the model:
sim.createWall (
 name = "left_wall",
 posn = Vec3(0.0, 0.0, 0.0),
 normal = Vec3(1.0, 0.0, 0.0)
)

#create a wall at the right side of the model:
sim.createWall (
 name = "right_wall",
 posn = Vec3(width, 0.0, 0.0),
 normal = Vec3(-1.0, 0.0, 0.0)
)

bond_param_brittle = BrittleBeamPrms(
 name="pp_bonds",
 youngsModulus=20.0e9 * scale, # 20 GPa, adjusted for scale
 poissonsRatio=0.25,
 cohesion=4.0e7 * scale, # 4 MPa, adjusted for scale
 tanAngle=0.577, # 30 degrees
 tag=1
)

pp_bonds = sim.createInteractionGroup (
 bond_param_brittle
)

# frictional interactions for regular particles
sim.createInteractionGroup (
 prms = FrictionPrms(
  name="friction",
  youngsModulus=20.0e9 * scale, # 20 GPa, adjusted for scale
  poissonsRatio=0.25,
  dynamicMu=0.4,
  staticMu=0.6
 )
)

#create an exclusion between bonded and frictional interactions:
sim.createExclusion (
 interactionName1 = "pp_bonds",
 interactionName2 = "friction"
)

# add gravity
sim.createInteractionGroup(
 GravityPrms(name="mercury-gravity", acceleration=Vec3(0,-3.7 / scale,0)) # Mercury's gravity acceleration, adjusted for scale
)

#add translational viscous damping:
sim.createInteractionGroup (
 LinDampingPrms(
  name="damping1",
  viscosity=0.1,
  maxIterations=50
 )
)

#add rotational viscous damping:
sim.createInteractionGroup (
 RotDampingPrms(
  name="damping2",
  viscosity=0.1,
  maxIterations=50
 )
)

#specify elastic repulsion from the bottom wall:
sim.createInteractionGroup (
 NRotElasticWallPrms (
  name = "bottom_wall_repel",
  wallName = "bottom_wall",
  normalK = 1.0e14 * scale
 )
)

sim.createInteractionGroup (
 NRotBondedWallPrms (
  name = "left_wall_bond",
  wallName = "left_wall",
  normalK = 1.0e14 * scale,
  particleTag = 1,
  tagMask = 0
 )
)

sim.createInteractionGroup (
 NRotElasticWallPrms (
  name = "left_wall_repel",
  wallName = "left_wall",
  normalK = 1.0e14 * scale
 )
)

sim.createInteractionGroup (
 NRotBondedWallPrms (
  name = "right_wall_bond",
  wallName = "right_wall",
  normalK = 1.0e14 * scale,
  particleTag = 1,
  tagMask = 0
 )
)

sim.createInteractionGroup (
 NRotElasticWallPrms (
  name = "right_wall_repel",
  wallName = "right_wall",
  normalK = 1.0e14 * scale
 )
)

#add a wall loader to move the left wall:
wall_loader_left = WallLoaderRunnable(
 LsmMpi = sim,
 wallName = "left_wall",
 vPlate = Vec3 (80000 / scale / timeStepSize, 0.0, 0.0),
 startTime = 0,
 rampTime = 1000
)
sim.addPreTimeStepRunnable (wall_loader_left)

#add a CheckPointer to save simulation data at regular intervals:
sim.createCheckPointer (
 CheckPointPrms (
  fileNamePrefix = "/home/rubio/esys-particle/thrust2D/checkpoints/checkpoint",
  beginTimeStep = 0,
  endTimeStep = timeSteps,
  timeStepIncr = timeStepIncr
 )
)

#execute the simulation:
sim.run()

--------------------------------------------------------------------------------