diff --git a/Project.toml b/Project.toml index 32730ab..8f52c0b 100644 --- a/Project.toml +++ b/Project.toml @@ -2,3 +2,10 @@ name = "MetcalfeGeothermal" uuid = "d178e5ed-affd-4f17-ab72-b61884852ee4" authors = ["Bob Metcalfe 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" diff --git a/src/MetcalfeGeothermal.jl b/src/MetcalfeGeothermal.jl index f294a46..72cc44e 100644 --- a/src/MetcalfeGeothermal.jl +++ b/src/MetcalfeGeothermal.jl @@ -1,5 +1,5 @@ module MetcalfeGeothermal - +using Plots println("Toy Model of Coaxial Geothermal - Bob.Metcalfe@UTexas.edu 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 diff --git a/src/PDESystem/coaxial_geothermal_MOL.jl b/src/PDESystem/coaxial_geothermal_MOL.jl new file mode 100644 index 0000000..66decbe --- /dev/null +++ b/src/PDESystem/coaxial_geothermal_MOL.jl @@ -0,0 +1,116 @@ +using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets, Plots + +println("1D PDE Depth Model of Coaxial Geothermal - Bob.Metcalfe@UTexas.edu 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") diff --git a/src/PDESystem/coaxial_geothermal_MOL_2D.jl b/src/PDESystem/coaxial_geothermal_MOL_2D.jl new file mode 100644 index 0000000..b61d1d2 --- /dev/null +++ b/src/PDESystem/coaxial_geothermal_MOL_2D.jl @@ -0,0 +1,88 @@ +using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets, Plots + +println("1D PDE Model of Coaxial Geothermal - Bob.Metcalfe@UTexas.edu 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] diff --git a/src/PDESystem/metcalfe_geothermal.png b/src/PDESystem/metcalfe_geothermal.png new file mode 100644 index 0000000..419123c Binary files /dev/null and b/src/PDESystem/metcalfe_geothermal.png differ diff --git a/src/PDESystem/metcalfecoaxialgeothermal.pdf b/src/PDESystem/metcalfecoaxialgeothermal.pdf new file mode 100644 index 0000000..7ce8a6c Binary files /dev/null and b/src/PDESystem/metcalfecoaxialgeothermal.pdf differ