Add 1D MOL implementation, with draft of 2D

Project.toml

@@ -2,3 +2,10 @@ name = "MetcalfeGeothermal"
 uuid = "d178e5ed-affd-4f17-ab72-b61884852ee4"
 authors = ["Bob Metcalfe <[email protected]> and contributors"]
 version = "0.1.0"
+
+[deps]
+DomainSets = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
+MethodOfLines = "94925ecb-adb7-4558-8ed8-f975c56a0bf4"
+ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"

src/MetcalfeGeothermal.jl

@@ -1,5 +1,5 @@
 module MetcalfeGeothermal
-
+using Plots
 println("Toy Model of Coaxial Geothermal - [email protected] TexasGEO.org")
 # ---------- Earth
 const earth = (EarthSurfaceTemperature = 15.0, EarthTemperatureGradient = 0.025) # [C], [C/m]
@@ -43,7 +43,7 @@ const PumpSpeed = InnerVolume/5.0 # [m3/s]
 const PumpTime = InnerVolume/PumpSpeed # [s]
 # ---------- Operation
 for run in 1:100 # iteration on thermo updates
-    
+
     for ip in 1:NumberOfPipes
         pipe = PipeString[ip]
         # outer pipe to inner pipe Joule flux
@@ -52,16 +52,16 @@ for run in 1:100 # iteration on thermo updates
         O2Ia = InnerArea
         O2Im = InnerVolume * FluidDensity
         O2Ipt = PumpTime
-        
+
         O2Ijf = O2Ic * O2Idt * O2Ia * O2Im * O2Ipt
         # print("O2I ",run, O2Ijf,O2Ic,O2Idt,O2Ia,O2Im,O2Ipt)
-        # ambient to outer Joule flux 
+        # ambient to outer Joule flux
         A2Oc = RockThermalConductivity
         A2Odt = AmbientTemperature(earth, (ip-1)*LengthOfPipe) - pipe.outertemp
         A2Oa = OuterArea
         A2Om = OuterVolume * FluidDensity
         A2Opt = PumpTime
-        
+
         A2Ojf = A2Oc * A2Odt * A2Oa * A2Om * A2Opt
         # print("A2O ",run,A2Ojf,A2Oc,A2Odt,A2Oa,A2Om,A2Opt)
         # Update temperatures
@@ -91,4 +91,6 @@ PipeString[1].outertemp = savedinnerpipe
 PipeString[NumberOfPipes].innertemp = savedouterpipe
 printPipeString("----------- stop -----------")
 
+plot([pipe.innertemp for pipe in PipeString], label="inner")
+plot!([pipe.outertemp for pipe in PipeString], label="outer")
 end # module MetcalfeGeothermal

src/PDESystem/coaxial_geothermal_MOL.jl

@@ -0,0 +1,116 @@
+using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets, Plots
+
+println("1D PDE Depth Model of Coaxial Geothermal - [email protected] TexasGEO.org")
+# ---------- Earth
+const earth = (EarthSurfaceTemperature=15.0, EarthTemperatureGradient=0.025) # [C], [C/m]
+# Ambient temperature in Celsius at depth in meters
+function AmbientTemperature(earth, depth)
+    earth.EarthSurfaceTemperature + depth * earth.EarthTemperatureGradient
+end
+# Rock is Granite
+const RockDensity = GraniteDensity = 2750000 # [g/m3]
+const RockSpecificHeat = GraniteSpecificHeat = 0.790 # [J/gC]
+const RockThermalConductivity = GraniteThermalConductivity = 2.62 # [W/mC] range average
+# Fluid is Water
+const FluidDensity = WaterDensity = 997000.0 # [g/m3]
+const FluidSpecificHeat = WaterSpecificHeat = 4.186 # [J/gC]
+const FluidThermalConductivity = WaterThermalConductivity = 0.6 # [W/mC]
+# ---------- Drilling
+# Volume of inner and outer pipes is equal for now
+const NumberOfPipes = 15 # [n]
+const LengthOfPipe = 10.0 # [m]
+const DepthOfWell = NumberOfPipes * LengthOfPipe # [m]
+const InnerRadius = 0.1 # [m]
+const OuterRadius = InnerRadius * sqrt(2) # [m]
+println("Well depth ", DepthOfWell, " meters, pipe radius ", InnerRadius, " meters.")
+const InnerArea = 2 * InnerRadius * pi * LengthOfPipe # [m2]
+const OuterArea = 2 * OuterRadius * pi * LengthOfPipe # [m2]
+const InnerVolume = pi * InnerRadius^2 * LengthOfPipe # [m3]
+const OuterVolume = (pi * OuterRadius^2 * LengthOfPipe) - InnerVolume # [m3]
+# Check to be sure inner and outer pipe volumes are equal (for now)
+if round(InnerVolume, digits=4) != round(OuterVolume, digits=4)
+    println("Inner and outer volumnes are not equal.")
+end
+
+const PumpSpeed = InnerVolume /5 # [m3/s]
+const PumpTime = InnerVolume / PumpSpeed # [s]
+const WaterSpeed = LengthOfPipe/PumpTime # [m/s]
+
+# Let the drilling begin! - With MethodOfLines.jl
+############################################################################################
+# 1D domain: z
+############################################################################################
+
+println("Water speed ", WaterSpeed, " meters per second.")
+
+
+@parameters t z
+@variables Tinner(..), Touter(..)
+
+Dz = Differential(z)
+Dt = Differential(t)
+
+# Let the drilling begin
+
+@parameters t z
+@variables Tinner(..), Touter(..)
+
+Dz = Differential(z)
+Dt = Differential(t)
+
+# Define constants
+V = WaterSpeed #[m / s]
+c_p = FluidSpecificHeat
+ρ = FluidDensity
+
+u_inner = 2π * InnerRadius
+u_outer = 2π * OuterRadius
+
+# ! These are far too small, we need to find another way to get the correct values
+α_inner = FluidThermalConductivity/InnerRadius
+α_outer = RockThermalConductivity/OuterRadius
+
+A_inner = π * InnerRadius^2 # cross sectional area of inner pipe
+A_outer = π * (OuterRadius^2 - InnerRadius^2) # cross sectional area of outer pipe
+
+tmax = 3000.0
+
+# Structure defined in the PDF
+eqs = [A_inner * ρ * c_p * Dt(Tinner(t, z)) - V * Dz(Tinner(t, z)) ~ -u_inner * α_inner * (Tinner(t, z) - Touter(t, z)),
+    A_outer * ρ * c_p * Dt(Touter(t, z)) + V * Dz(Touter(t, z)) ~ u_inner * α_inner * (Tinner(t, z) - Touter(t, z)) +
+                                                              u_outer * α_outer * (AmbientTemperature(earth, z) - Touter(t, z)),
+]
+
+bcs = [Tinner(0, z) ~ AmbientTemperature(earth, z),
+    Touter(0, z) ~ AmbientTemperature(earth, z),
+    Touter(t, 0) ~ AmbientTemperature(earth, 0),
+    Tinner(t, DepthOfWell) ~ Touter(t, DepthOfWell),
+]
+
+domains = [z in Interval(0.0, DepthOfWell),
+    t in Interval(0.0, tmax),
+]
+
+@named pde = PDESystem(eqs, bcs, domains, [t, z], [Tinner(t, z), Touter(t, z)])
+
+# Create the discretization
+
+disc = MOLFiniteDifference([z => 100], t)
+
+# Generate the ODE problem
+prob = discretize(pde, disc)
+
+# Solve the ODE problem
+sol = solve(prob, QBDF(), saveat=10.0)
+
+sol_Tinner = sol[Tinner(t, z)]
+sol_Touter = sol[Touter(t, z)]
+
+solt = sol[t]
+solz = sol[z]
+
+println("At t = $tmax, the inner pipe temperature at ground level is ", sol_Tinner[end, 1], " and the outer pipe temperature is ", sol_Touter[end, 1], ".")
+
+# Plot the solution
+p = plot(solz, [sol_Tinner[end, 1:end], sol_Touter[end, 1:end]], label=["Tinner" "Touter"], xlabel="z", ylabel="T", title="Temperature vs. Depth")
+savefig(p, "src/PDESystem/metcalfe_geothermal.png")

src/PDESystem/coaxial_geothermal_MOL_2D.jl

@@ -0,0 +1,88 @@
+using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets, Plots
+
+println("1D PDE Model of Coaxial Geothermal - [email protected] TexasGEO.org")
+# ---------- Earth
+const earth = (EarthSurfaceTemperature=15.0, EarthTemperatureGradient=0.025) # [C], [C/m]
+# Ambient temperature in Celsius at depth in meters
+function AmbientTemperature(earth, depth)
+    earth.EarthSurfaceTemperature + depth * earth.EarthTemperatureGradient
+end
+# Rock is Granite
+const RockDensity = GraniteDensity = 2750000 # [g/m3]
+const RockSpecificHeat = GraniteSpecificHeat = 0.790 # [J/gC]
+const RockThermalConductivity = GraniteThermalConductivity = 2.62 # [W/mC] range average
+# Fluid is Water
+const FluidDensity = WaterDensity = 997000.0 # [g/m3]
+const FluidSpecificHeat = WaterSpecificHeat = 4.186 # [J/gC]
+const FluidThermalConductivity = WaterThermalConductivity = 0.6 # [W/mC]
+# ---------- Drilling
+# Volume of inner and outer pipes is equal for now
+const NumberOfPipes = 15 # [n]
+const LengthOfPipe = 10.0 # [m]
+const DepthOfWell = NumberOfPipes * LengthOfPipe # [m]
+const InnerRadius = 0.1 # [m]
+const OuterRadius = InnerRadius * sqrt(2) # [m]
+println("Well depth ", DepthOfWell, " meters, pipe radius ", InnerRadius, " meters.")
+const InnerArea = 2 * InnerRadius * pi * LengthOfPipe # [m2]
+const OuterArea = 2 * OuterRadius * pi * LengthOfPipe # [m2]
+const InnerVolume = pi * InnerRadius^2 * LengthOfPipe # [m3]
+const OuterVolume = (pi * OuterRadius^2 * LengthOfPipe) - InnerVolume # [m3]
+# Check to be sure inner and outer pipe volumes are equal (for now)
+if round(InnerVolume, digits=4) != round(OuterVolume, digits=4)
+    println("Inner and outer volumnes are not equal.")
+end
+
+const PumpSpeed = InnerVolume / 5.0 # [m3/s]
+const PumpTime = InnerVolume / PumpSpeed # [s]
+
+############################################################################################
+# cylindrical domain **WIP**
+############################################################################################
+
+## Now try to do the same thing with radial gradients and a cylindrical domain
+@parameters z r t
+@variables T(..)
+
+Dt = Differential(t)
+Dz = Differential(z)
+Dr = Differential(r)
+Drr = Differential(r)^2
+
+# r dependent velocity
+v(x) = ifelse(x < InnerRadius, PumpSpeed, -PumpSpeed)
+
+# ? TODO: Define constants
+
+# Single equation this time
+eqs = [K * Dt(T(t, z, r)) + v(r) * Dz(T(t, z, r)) ~ u_inner * α_inner * 1 / r^2 * Dr(r^2 * Dr(T(t, z, r)))]
+
+bcs = [T(0, z, r) ~ AmbientTemperature(earth, z),
+    T(0, z, r) ~ AmbientTemperature(earth, z),
+    T(t, z, 0) ~ 6 * Drr(T(u, z, r)), # from the paper https://web.mit.edu/braatzgroup/analysis_of_finite_difference_discretization_schemes_for_diffusion_in_spheres_with_variable_diffusivity.pdf
+    T(t, z, OuterRadius) ~ AmbientTemperature(earth, z),
+    T(t, 0, r) ~ AmbientTemperature(earth, 0),
+    Touter(t, DepthOfWell, r) ~ Tinner(t, DepthOfWell, r),
+]
+
+domains = [z in Interval(0.0, DepthOfWell),
+    r in Interval(0.0, OuterRadius),
+    t in Interval(0.0, 100.0),
+]
+
+@named pde = PDESystem(eqs, bcs, domains, [t, z, r], [T(t, z, r)])
+
+# Create the discretization
+
+disc = MOLFiniteDifference([z => 32, r => 16], t)
+
+# Generate the ODE problem
+prob = discretize(pde, disc)
+
+# Solve the ODE problem
+sol = solve(prob, QBDF(), saveat=10.0)
+
+sol_Tinner = sol[Tinner(t, z, r)]
+sol_Touter = sol[Touter(t, z, r)]
+
+solt = sol[t]
+solz = sol[z]