Romain Beucher Version 0.1
romain.beucher@unimelb.edu.au
The following notebook is an implementation of the Numerical Sandbox Extension Experiment similar to Buiter et al., 2006. The test is commonly referred as one of the GEOMOD benchmarks and is used to test the large deformation viscous-plastic behaviour of geodynamic numerical codes.
The initial purpose of the GEOMOD numerical experiments was to compare results from analog and numerical experiments. It is not a numerical benchmark sensus-stricto.
The test has been implemented using Fantom in Thieulot, 2011 and more recently with Aspect in Glerum et al., 2017.
Results for Underworld from this notebook
A | B |
---|---|
Results for the extension experiment after 2cm extension (Buiter et al., 2006)
import UWGeodynamics as GEO
loaded rc file /workspace/user_data/UWGeodynamics/UWGeodynamics/uwgeo-data/uwgeodynamicsrc
u = GEO.UnitRegistry
velocity = 2.5 * u.centimeter / u.hour
model_length = 20. * u.centimeter
bodyforce = 1560 * u.kilogram / u.metre**3 * 9.81 * u.meter / u.second**2
KL = model_length
Kt = KL / velocity
KM = bodyforce * KL**2 * Kt**2
GEO.scaling_coefficients["[length]"] = KL
GEO.scaling_coefficients["[time]"] = Kt
GEO.scaling_coefficients["[mass]"]= KM
Set-up of the extension experiment. A viscous layer (PDMS, 10 x 0.5 cm) lies in the central part of the model on the base. The rest of the model consists of three ‘sand’ layers (only differing in colour). Extension is achieved by moving the right wall with the attached 10 cm long sheet outwards to the right.
We use a uniform resolution of 0.2 cm (400 x 100 elements)
Model = GEO.Model(elementRes=(200, 50),
minCoord=(0. * u.centimeter, -3.5 * u.centimeter),
maxCoord=(20. * u.centimeter, 1.5 * u.centimeter),
gravity=(0.0, -9.81 * u.meter / u.second**2))
Model.mesh_advector(axis=0)
Model.outputDir="outputs_tutorial3B"
Model.add_visugrid(elementRes=(100, 25),
minCoord=(0. * u.centimeter, -3.5 * u.centimeter),
maxCoord=(20. * u.centimeter, 0.0 * u.centimeter))
Now that we have our "universe" (box, sand pit) ready, we need to fill it with some materials. The geodynamics module is designed around that idea of materials, which are essentially a way to define physical properties across the Model domain.
A material (or a phase) is first defined by the space it takes in the box (its shape).
air = Model.add_material(name="Air", shape=GEO.shapes.Layer(top=Model.top, bottom=Model.bottom))
sand1 = Model.add_material(name="Sand1", shape=GEO.shapes.Layer(top=0.*u.centimeter, bottom=Model.bottom))
sand2 = Model.add_material(name="Sand2", shape=GEO.shapes.Layer(top=-1. * u.centimeter, bottom=-2. * u.centimeter))
vertices = [( 5.* u.centimeter, -3.0 * u.centimeter),
(15.* u.centimeter, -3.0 * u.centimeter),
(15.* u.centimeter, -3.5 * u.centimeter),
( 5.* u.centimeter, -3.5 * u.centimeter)]
silicon = Model.add_material(name="Silicon", shape=GEO.shapes.Polygon(vertices))
import numpy as np
x = np.linspace(GEO.nd(Model.minCoord[0]), GEO.nd(Model.maxCoord[0]), 1000)
interface1 = GEO.nd(sand1.top)
interface2 = GEO.nd(sand2.top)
interface3 = GEO.nd(sand2.bottom)
interface1 = Model.add_passive_tracers(name="Interface1", vertices=[x, interface1])
interface2 = Model.add_passive_tracers(name="Interface2", vertices=[x, interface2])
interface3 = Model.add_passive_tracers(name="Interface3", vertices=[x, interface3])
import glucifer
Fig = glucifer.Figure(figsize=(1200,400))
Fig.Points(interface1, pointSize=2.0)
Fig.Points(interface2, pointSize=2.0)
Fig.Points(interface3, pointSize=2.0)
Fig.Points(Model.swarm, Model.materialField, fn_size=2.0)
Fig.show()
air.density = 10. * u.kilogram / u.metre**3
sand1.density = 1560. * u.kilogram / u.metre**3
sand2.density = 1560. * u.kilogram / u.metre**3
silicon.density = 965. * u.kilogram / u.metre**3
Viscosities can be defined as a Quantity or a simple scalar value. It is also possible to load predefined rheologies from the rheology_library.
The rheology library contains some commonly used rheologies stored in a python dictionary structure. We can list the keys defining the rheologies as follows:
# Assign function to materials
air.viscosity = 1.0e2 * u.pascal * u.second
sand1.viscosity = 1.0e13 * u.pascal * u.second
sand2.viscosity = 1.0e13 * u.pascal * u.second
silicon.viscosity = 5.0e4 * u.pascal * u.second
Plastic behavior is assigned using the same approach as for viscosities.
sandPlasticity = GEO.DruckerPrager(cohesion=10. * u.pascal,
cohesionAfterSoftening=10. * u.pascal,
frictionCoefficient=0.73,
frictionAfterSoftening=0.60)
sand1.plasticity = sandPlasticity
sand2.plasticity = sandPlasticity
import underworld.function as fn
conditions = [(Model.x > GEO.nd(10.1 * u.centimetre), GEO.nd(2.5 * u.centimeter / u.hour)),
(Model.x > GEO.nd(9.9 * u.centimetre), (Model.x - GEO.nd(9.9 * u.centimetre)) * GEO.nd(velocity) / GEO.nd(0.2 * u.centimetre)),
(True, 0.0)]
fn_condition = fn.branching.conditional(conditions)
Model.set_velocityBCs(left=[0 * u.centimeter / u.hour, None],
right=[2.5 * u.centimeter / u.hour, None],
bottom=[fn_condition, 0.0])
<UWGeodynamics._velocity_boundaries.VelocityBCs at 0x7f07b453d9e8>
The model is initialize using the init_model method which will solve the initial steady state temperature field and the pressure field.
Model.minViscosity = 1.0e2 * u.pascal * u.second
Model.maxViscosity = 1.0e9 * u.pascal * u.second
Model.init_model()
FigDensity = glucifer.Figure(figsize=(1200, 400), title="Density Field (kg/m3)", quality=3)
FigDensity.Points(Model.swarm, GEO.dimensionalise(Model.densityField, u.kilogram / u.metre**3), pointSize=2.0)
FigDensity.show()
FigPressure = glucifer.Figure(figsize=(1200, 400), title="Pressure Field (Pa)", quality=3)
FigPressure.Surface(Model.mesh, GEO.dimensionalise(Model.pressureField, u.pascal))
FigPressure.show()
FigViscosity = glucifer.Figure(figsize=(1200, 400), title="Viscosity Field (Pa.s)", quality=3)
FigViscosity.Points(Model.swarm, GEO.dimensionalise(Model.viscosityField, u.pascal * u.second), pointSize=2.0)
FigViscosity.show()
Model.solver.set_inner_method("mumps")
Model.solver.set_penalty(1e6)
Model.run_for(1.0 * u.hours, checkpoint_interval=10 * u.minutes)