commit 7a535a7e52e53c1d5967573ab4b97a062198c120
parent 5e5b568af441c1457133afd2deb7e09867fbb31e
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Fri, 5 Apr 2019 08:47:59 +0200
Run over range of effective normal stresses
Diffstat:
| M | 1d_fd_simple_shear.jl | | | 80 | +++++++++++++++++++++++++++++++++---------------------------------------------- |
1 file changed, 33 insertions(+), 47 deletions(-)
diff --git a/1d_fd_simple_shear.jl b/1d_fd_simple_shear.jl
@@ -3,51 +3,54 @@
ENV["MPLBACKEND"] = "Agg"
import PyPlot
+let
+
# Simulation settings
## Gravitational acceleration magnitude
-const G = 9.81
+G = 9.81
## Wall parameters
# A larger normal stress deepens the deformational depth
-const P_wall = 100.0e3 # normal stress [Pa]
+P_wall_ = [10, 20, 40, 60, 80, 120] .* 1e3 # normal stress [Pa]
-const v_x_bot = 0.0 # bottom velocity along x [m/s]
-const μ_wall = 0.40 # stress ratio at top wall
+v_x_bot = 0.0 # bottom velocity along x [m/s]
+μ_wall = 0.40 # stress ratio at top wall
## Nodes along z
-const nz = 100
+nz = 100
## Material properties
# lower values of A mean that the velocity curve can have sharper curves, e.g.
# at the transition from μ ≈ μ_s
-const A = 0.48 # nonlocal amplitude [-]
+A = 0.48 # nonlocal amplitude [-]
# lower values of b mean larger shear velocity for a given stress ratio above
# yield
-const b = 0.9377 # rate dependence beyond yield [-]
+b = 0.9377 # rate dependence beyond yield [-]
-const μ_s = 0.3819 # static yield friction coefficient [-]
-const ϕ = 0.38 # porosity [-]
+μ_s = 0.3819 # static yield friction coefficient [-]
+μ_s = atan(deg2rad(22.0)) # Damsgaard et al 2013
+ϕ = 0.38 # porosity [-]
# lower values of d mean that the shear velocity curve can have sharper curves,
# e.g. at the transition from μ ≈ μ_s
-const d = 1e-3 # representative grain size [m]
+#const d = 1e-3 # representative grain size [m]
+d = 0.04 # representative grain size [m] - Damsgaard et al 2013
-const ρ_s = 2.485e3 # grain material density [kg/m^3]
+ρ_s = 2.485e3 # grain material density [kg/m^3]
## Domain size
-const origo_z = 0.0
+origo_z = 0.0
#const L_z = 20.0*d
-#const L_z = 200.0*d
-const L_z = 0.5
+L_z = 0.7
## Spatial coordiantes
-const z = collect(range(origo_z, L_z, length=nz))
-const Δz = z[2] - z[1]
+z = collect(range(origo_z, L_z, length=nz))
+Δz = z[2] - z[1]
## Allocate other arrays
μ = zero(z) # local stress ratio
@@ -61,9 +64,6 @@ g_ghost = zeros(size(z)[1]+2) # local fluidity with ghost nodes
## Shear plastic strain rate (eq 2)
γ_dot_p(g, μ) = g.*μ
-## Stress ratio at the upper wall
-#μ_wall(τ_wall, P_wall) = τ_wall.*P_wall
-
## Normal stress
p_lithostatic(P_wall, z) = P_wall .+ (1 - ϕ).*ρ_s.*G.*(L_z .- z)
@@ -104,18 +104,6 @@ end
#friction(τ, p) = τ./(p .+ 1e-16)
-## Fluidity
-#fluidity(p, μ, g) = 1.0./(ξ.(μ).^2.0) .* (g .- g_local.(p, μ))
-
-# ## set up the forcing function for the Poisson equation. The forcing is the
-# ## right-hand side of the fluidity equation
-# function forcing(p, μ, g_ghost)
-# return fluidity.(μ, g_ghost[2:end-1], p)
-# end
-
-## 1d Laplacian operator (2nd order central finite difference)
-#laplacian(phi, i, dx) = (phi[i-1] - 2.0*phi[i] + phi[i+1])/dx
-
## A single iteration for solving a Poisson equation Laplace(phi) = f on a
## Cartesian grid. The function returns the normalized residual value
function poisson_solver_1d_iteration(g_in, g_out, r_norm,
@@ -205,7 +193,6 @@ end
function plot_profile(z, v, label, filename)
PyPlot.clf()
PyPlot.plot(v, z, "+k")
- #PyPlot.gca().invert_yaxis()
PyPlot.xlabel(label)
PyPlot.ylabel("\$z\$ [m]")
PyPlot.tight_layout()
@@ -219,23 +206,22 @@ init_μ(μ_wall, ϕ, ρ_s, G, z, P_wall) =
# Main
-## Set velocity BCs
-v_x[1] = v_x_bot
+for P_wall in P_wall_
+ ## Set velocity BCs
+ v_x[1] = v_x_bot
-## Calculate stressses
-p = p_lithostatic(P_wall, z)
-μ = init_μ(μ_wall, ϕ, ρ_s, G, z, P_wall)
-#τ = shear_stress(v_x, Δz)
+ ## Calculate stressses
+ p = p_lithostatic(P_wall, z)
+ μ = init_μ(μ_wall, ϕ, ρ_s, G, z, P_wall)
-#f = forcing(p, μ, g_ghost)
-#g_ghost[2:end-1] = f
+ implicit_1d_jacobian_poisson_solver(g_ghost, p, μ, Δz)
-implicit_1d_jacobian_poisson_solver(g_ghost, p, μ, Δz)
+ v_x = γ_dot_p(g_ghost[2:end-1], μ)
-v_x = γ_dot_p(g_ghost[2:end-1], μ)
+ plot_profile(z, v_x, "Shear velocity, \$v_x\$ [m/s]", "1d_fd_simple_shear_v_x_P$(P_wall).png")
+ plot_profile(z, μ, "Stress ratio, μ [-]", "1d_fd_simple_shear_mu_P$(P_wall).png")
+ plot_profile(z, p, "Normal stress, \$p\$ [Pa]", "1d_fd_simple_shear_p_P$(P_wall).png")
+ plot_profile(z, g_ghost[2:end-1], "Fluidity, \$g\$", "1d_fd_simple_shear_g_P$(P_wall).png")
+end
-plot_profile(z, v_x, "Shear velocity, \$v_x\$ [m/s]", "1d_fd_simple_shear_v_x.png")
-plot_profile(z, μ, "Stress ratio, μ [-]", "1d_fd_simple_shear_mu.png")
-plot_profile(z, p, "Normal stress, \$p\$ [Pa]", "1d_fd_simple_shear_p.png")
-#plot_profile(z, f, "Forcing, \$f\$", "1d_fd_simple_shear_f.png")
-plot_profile(z, g_ghost[2:end-1], "Fluidity, \$g\$", "1d_fd_simple_shear_g.png")
+end # end let
