cngf-pf-exp1

commit 63709c7bdabc7953803e85e4b1ac902388f619de
parent b944650f811707194c52ed1621b4672d3caaa472
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date:   Thu, 27 Jun 2019 15:04:14 +0200

Add analytical solution to pressure variation

Diffstat:
A
p_f_analytical.jl
 | 
27
+++++++++++++++++++++++++++

1 file changed, 27 insertions(+), 0 deletions(-)
diff --git a/p_f_analytical.jl b/p_f_analytical.jl
@@ -0,0 +1,27 @@
+#!/usr/bin/env julia
+import PyPlot
+
+n = 50                          # resolution
+z = range(0, stop=2, length=n)  # spatial grid
+k = 2e-17                       # permeability
+ϕ = 0.25                        # porosity
+μ_f = 1e-3                      # water viscosity
+β_f = 4.5e-10                   # water compressibility
+κ = k/(ϕ*μ_f*β_f)               # diffusivity
+ω = 2*π*1/(3600*24)             # diurnal circular frequency
+
+# T is water pressure in kPa
+T0 = 50                         # initial pressure
+ΔT = 50                         # pressure perturbation amplitude
+
+# Turcotte and Schubert, eq. 4.89
+T_(z,t) = T0 .+ ΔT.*exp.(-z.*sqrt(ω/(2κ))).*cos.(ω*t .- z.*sqrt(ω/(2κ)))
+
+# plot hourly curves for a full day
+for t=0:3600:24*3600
+	PyPlot.plot(reverse(T_(z,t)), z)
+end
+
+PyPlot.xlabel("Water pressure [kPa]")
+PyPlot.ylabel("Vertical position [m]")
+PyPlot.savefig("p_f_analytical.pdf")