Merge pull request #80 from UM-PEPL/plume_divergence

Plume divergence and some fixes for mass conservation

Manifest.toml

@@ -2,7 +2,7 @@
 
 julia_version = "1.9.0-beta2"
 manifest_format = "2.0"
-project_hash = "6dc93f3c34307df9551e9e0376ceb5119b058357"
+project_hash = "16c3410dc94142d4a6d58ae6d032899994167fed"
 
 [[deps.ANSIColoredPrinters]]
 git-tree-sha1 = "574baf8110975760d391c710b6341da1afa48d8c"
@@ -133,12 +133,6 @@ git-tree-sha1 = "ec8a9c9f0ecb1c687e34c1fda2699de4d054672a"
 uuid = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 version = "0.4.29"
 
-[[deps.JLLWrappers]]
-deps = ["Preferences"]
-git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1"
-uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
-version = "1.4.1"
-
 [[deps.JSON]]
 deps = ["Dates", "Mmap", "Parsers", "Unicode"]
 git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e"
@@ -234,17 +228,6 @@ deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
 uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
 version = "0.3.21+0"
 
-[[deps.OpenLibm_jll]]
-deps = ["Artifacts", "Libdl"]
-uuid = "05823500-19ac-5b8b-9628-191a04bc5112"
-version = "0.8.1+0"
-
-[[deps.OpenSpecFun_jll]]
-deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
-uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
-version = "0.5.5+0"
-
 [[deps.OrderedCollections]]
 git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
 uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
@@ -333,12 +316,6 @@ version = "1.1.0"
 deps = ["Libdl", "LinearAlgebra", "Random", "Serialization", "SuiteSparse_jll"]
 uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
-[[deps.SpecialFunctions]]
-deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
-git-tree-sha1 = "d75bda01f8c31ebb72df80a46c88b25d1c79c56d"
-uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
-version = "2.1.7"
-
 [[deps.Statistics]]
 deps = ["LinearAlgebra", "SparseArrays"]
 uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

Project.toml

@@ -1,7 +1,7 @@
 name = "HallThruster"
 uuid = "2311f341-5e6d-4941-9e3e-3ce0ae0d9ed6"
 authors = ["Thomas Marks <[email protected]>"]
-version = "0.5.2"
+version = "0.6.0"
 
 [deps]
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
@@ -19,7 +19,6 @@ QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 SnoopPrecompile = "66db9d55-30c0-4569-8b51-7e840670fc0c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
-SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"

src/HallThruster.jl

@@ -1,6 +1,5 @@
 module HallThruster
 
-using SpecialFunctions
 using LinearAlgebra
 using DocStringExtensions
 
@@ -61,6 +60,7 @@ include("simulation/potential.jl")
 include("simulation/heavy_species.jl")
 include("simulation/electronenergy.jl")
 include("simulation/sourceterms.jl")
+include("simulation/plume.jl")
 include("simulation/update_electrons.jl")
 include("simulation/configuration.jl")
 include("simulation/restart.jl")

src/numerics/flux.jl

@@ -107,11 +107,11 @@ for NUM_CONSERVATIVE in 1:3
 
         charge_factor = Z * e * coupled
 
-         aL = sqrt((charge_factor * TeL + γ * kB * TL) / mi)
-         aR = sqrt((charge_factor * TeR + γ * kB * TR) / mi)
+        aL = sqrt((charge_factor * TeL + γ * kB * TL) / mi)
+        aR = sqrt((charge_factor * TeR + γ * kB * TR) / mi)
 
-        sL_max = max(abs(uL - aL), abs(uL + aL), #=abs(uL)=#)
-        sR_max = max(abs(uR - aR), abs(uR + aR), #=abs(uR)=#)
+        sL_max = max(abs(uL - aL), abs(uL + aL))
+        sR_max = max(abs(uR - aR), abs(uR + aR))
 
         smax = max(sL_max, sR_max)
 
@@ -192,7 +192,7 @@ function reconstruct(uⱼ₋₁, uⱼ, uⱼ₊₁, limiter)
     return uⱼ₋½ᴿ, uⱼ₊½ᴸ
 end
 
-function compute_edge_states!(UL, UR, U, params)
+function compute_edge_states!(UL, UR, U, params; apply_boundary_conditions = false)
     (nvars,  ncells) = size(U)
     (;config, index) = params
     (;scheme) = config
@@ -234,11 +234,16 @@ function compute_edge_states!(UL, UR, U, params)
         end
     end
 
-    @. @views UL[:, 1]   = U[:, 1]
-    @. @views UR[:, end] = U[:, end]
+    if apply_boundary_conditions
+        @views left_boundary_state!(UL[:, 1], U, params)
+        @views right_boundary_state!(UR[:, end], U, params)
+    else
+        @. @views UL[:, 1] = U[:, 1]
+        @. @views UR[:, end] = U[:, end]
+    end
 end
 
-function compute_fluxes!(F, UL, UR, U, params)
+function compute_fluxes!(F, UL, UR, U, params; apply_boundary_conditions = false)
     (;config, index, fluids, Δz_edge, z_cell, num_neutral_fluids) = params
     λ_global = params.cache.λ_global
     (;propellant, electron_pressure_coupled, scheme, ncharge) = config
@@ -252,7 +257,7 @@ function compute_fluxes!(F, UL, UR, U, params)
     params.max_timestep[1] = Inf
 
     # Reconstruct the states at the left and right edges using MUSCL scheme
-    compute_edge_states!(UL, UR, U, params)
+    compute_edge_states!(UL, UR, U, params; apply_boundary_conditions)
 
     # The contribution to the electron-pressure-coupled method will be 3/2 Te if we're using LANDMARK, since pe = nϵ in that benchmark
     if params.config.LANDMARK

src/simulation/boundaryconditions.jl

@@ -1,8 +1,8 @@
 function left_boundary_state!(bc_state, U, params)
-    (;index, A_ch, config,) = params
+    (;index, config,) = params
+    (;Tev, channel_area) = params.cache
     mi = config.propellant.m
     Ti = config.ion_temperature
-    (;Tev) = params.cache
 
     un = config.neutral_velocity
     mdot_a = config.anode_mass_flow_rate
@@ -25,13 +25,8 @@ function left_boundary_state!(bc_state, U, params)
         bohm_factor = 1.0
     end
 
-    bohm_factor = 1.0
-
-    # Precompute bohm velocity
-    bohm_velocity = bohm_factor * sqrt(e * Tev[1] / mi)
-
     # Add inlet neutral density
-    bc_state[index.ρn[1]] = mdot_a / A_ch / un
+    bc_state[index.ρn[1]] = mdot_a / channel_area[1] / un
 
     if config.solve_background_neutrals
         # Background neutrals are accomodated and re-emitted as anode neutrals
@@ -41,44 +36,46 @@ function left_boundary_state!(bc_state, U, params)
         bc_state[index.ρn[1]] += background_neutral_flux / un
         bc_state[index.ρn[2]] = U[index.ρn[2], begin+1]
     end
-
-    sound_speed = bohm_velocity #sqrt(config.propellant.γ * kB * Ti / mi)
-    boundary_velocity = -sound_speed
+
 
     @inbounds for Z in 1:params.config.ncharge
         interior_density = U[index.ρi[Z],   begin+1]
         interior_flux    = U[index.ρiui[Z], begin+1]
         interior_velocity = interior_flux / interior_density
 
+        sound_speed = sqrt((kB * Ti + Z * e * Tev[1]) / mi)  # Sound speed considering electron pressure-coupled terms
+        boundary_velocity = -bohm_factor * sound_speed # Want to drive flow to (negative) bohm velocity
+
         if interior_velocity <= -sound_speed
+            # Supersonic outflow, pure Neumann boundary condition
             boundary_density = interior_density
             boundary_flux = interior_flux
         else
+            # Subsonic outflow, need to drive the flow toward sonic
+            # For the isothermal Euler equations, the Riemann invariants are
+            # J⁺ = u + c ln ρ
+            # J⁻ = u - c ln ρ
+            # For the boundary condition, we take c = u_bohm
+
+            # 1. Compute outgoing characteristic using interior state
             J⁻ = interior_velocity - sound_speed * log(interior_density)
+
+            # 2. Compute incoming characteristic using J⁻ invariant
             J⁺ = 2 * boundary_velocity - J⁻
+
+            # 3. Compute boundary density using J⁺ and J⁻ invariants
             boundary_density = exp(0.5 * (J⁺ - J⁻) / sound_speed)
+
+            # Compute boundary flux
             boundary_flux = boundary_velocity * boundary_density
 
-            #@show interior_velocity, -sound_speed
-            #@show J⁺, J⁻
-            #@show interior_density, boundary_density
+            #boundary_flux = interior_flux
+            #boundary_density = boundary_flux / boundary_velocity
         end
 
         bc_state[index.ρn[1]] -= boundary_flux / un
         bc_state[index.ρi[Z]] = boundary_density
         bc_state[index.ρiui[Z]] = boundary_flux
-
-
-        #= Enforce Bohm condition at left boundary
-        boundary_velocity = min(-sqrt(Z) * bohm_velocity, boundary_velocity)
-
-        recombination_density = -(boundary_density * boundary_velocity) / un
-
-        bc_state[index.ρn[1]] += recombination_density
-        bc_state[index.ρi[Z]] = boundary_density # Neumann BC for ion density at left boundary
-        bc_state[index.ρiui[Z]] = boundary_velocity * boundary_density=#
-
-
     end
 end
 

src/simulation/configuration.jl

@@ -45,6 +45,7 @@ struct Config{A<:AnomalousTransportModel, W<:WallLossModel, IZ<:IonizationModel,
     background_neutral_temperature::Float64
     anode_boundary_condition::Symbol
     anom_smoothing_iters::Int
+    solve_plume::Bool
 end
 
 function Config(;
@@ -87,6 +88,7 @@ function Config(;
         background_neutral_temperature      = 100.0u"K",
         anode_boundary_condition            = :sheath,
         anom_smoothing_iters                = 0,
+        solve_plume = false
     ) where {IC, S_N, S_IC, S_IM, S_ϕ, S_E}
 
     # check that number of ion source terms matches number of charges for both
@@ -148,7 +150,8 @@ function Config(;
         background_pressure,
         background_neutral_temperature,
         anode_boundary_condition,
-        anom_smoothing_iters
+        anom_smoothing_iters,
+        solve_plume
     )
 end
 

src/simulation/electronenergy.jl

@@ -1,13 +1,14 @@
 const LOOKUP_ZS = 1.0:1.0:5.0
 const LOOKUP_CONDUCTIVITY_COEFFS = [4.66, 4.0, 3.7, 3.6, 3.2]
-const CONDUCTIVITY_LOOKUP_TABLE = LinearInterpolation(LOOKUP_ZS, LOOKUP_CONDUCTIVITY_COEFFS)
+const ELECTRON_CONDUCTIVITY_LOOKUP = LinearInterpolation(LOOKUP_ZS, LOOKUP_CONDUCTIVITY_COEFFS)
 
 function update_electron_energy!(U, params)
     (;Δz_cell, Δz_edge, dt, index, config, cache, Te_L, Te_R) = params
-    (;Aϵ, bϵ, μ, ue, ne, Tev) = cache
+    (;Aϵ, bϵ, μ, ue, ne, Tev, channel_area, dA_dz) = cache
     implicit = params.config.implicit_energy
     explicit = 1 - implicit
     ncells = size(U, 2)
+    mi = params.config.propellant.m
 
     nϵ = @views U[index.nϵ, :]
     Aϵ.d[1] = 1.0
@@ -18,10 +19,10 @@ function update_electron_energy!(U, params)
     if config.anode_boundary_condition == :dirichlet || ue[1] > 0
         bϵ[1] = 1.5 * Te_L * ne[1]
     else
-        # Neumann BC for internal energy
+        # Neumann BC for electron temperature
         bϵ[1] = 0
-        Aϵ.d[1] = 1.0
-        Aϵ.du[1] = -1.0
+        Aϵ.d[1] = 1.0 / ne[1]
+        Aϵ.du[1] = -1.0 / ne[2]
     end
 
     bϵ[end] = 1.5 * Te_R * ne[end]
@@ -61,7 +62,7 @@ function update_electron_energy!(U, params)
 
         else
             #get adjusted coeffient for higher charge states
-            κ_charge = CONDUCTIVITY_LOOKUP_TABLE(params.cache.Z_eff[i])
+            κ_charge = ELECTRON_CONDUCTIVITY_LOOKUP(params.cache.Z_eff[i])
             correction_factor = κ_charge/4.7
             # Adjust thermal conductivity to be slightly more accurate
             κL = 10/9 * 24/25 * (1 / (1 + params.cache.νei[i-1] / √(2) / params.cache.νc[i-1])) * μnϵL * correction_factor
@@ -90,13 +91,31 @@ function update_electron_energy!(U, params)
                 # left flux is sheath heat flux
                 Te0 = 2/3 * nϵ0 / ne0
 
-                uth = -0.25 * sqrt(8 * e * Te0 / π / me) * exp(-params.cache.Vs[] / Te0)
+                # Sheath heat flux = je_sheath_edge * (2 * Te_sheath + Vs) / e
+                # Assume Te_sheath is the same as that in first interior cell
+                # and that ne_sheath is that computed by using the bohm condition bc
+
+                # je_sheath_edg * (2 * Te_sheath + Vs) / e = - 2 ne ue Te - 2 ne Vs
+                # = - 4/3 * nϵ[1] * ue_sheath_edge - 2 ne ue Vs
+
+                # discharge current density
+                jd = params.cache.Id[] / channel_area[1]
+
+                # current densities at sheath edge
+                ji_sheath_edge = e * sum(Z * U[index.ρiui[Z], 1] for Z in 1:params.config.ncharge) / mi
+
+                je_sheath_edge = jd - ji_sheath_edge
+
+                #uth = -0.25 * sqrt(8 * e * Te0 / π / me) * exp(-params.cache.Vs[] / Te0)
+
+                ne_sheath_edge = sum(Z * U[index.ρi[Z], 1] for Z in 1:params.config.ncharge) / mi
+                ue_sheath_edge = -je_sheath_edge / ne_sheath_edge / e
 
                 FL_factor_L = 0.0
-                FL_factor_C = 4/3 * uth
+                FL_factor_C = 4/3 * ue_sheath_edge * (1 + params.cache.Vs[] / Te0)
                 FL_factor_R = 0.0
 
-                Q += ne0 * uth * params.cache.Vs[] / Δz
+                #Q = ne0 * ue_sheath_edge * params.cache.Vs[] / Δz
             elseif i == 2
                 # central differences at left boundary for compatibility with dirichlet BC
                 FL_factor_L = 5/3 * ueL + κL / ΔzL / neL
@@ -131,8 +150,13 @@ function update_electron_energy!(U, params)
         # Explicit flux
         F_explicit = (FR - FL) / Δz
 
+        # Term to allow for changing area
+        dlnA_dz = dA_dz[i] / channel_area[i]
+        flux = 5/3 * nϵ0 * ue0
+
         # Explicit right-hand-side
-        bϵ[i] = nϵ[i] + dt * (Q - explicit * F_explicit)
+        bϵ[i] = nϵ[i] + dt * (Q - explicit * F_explicit) 
+        bϵ[i] -= dt * flux * dlnA_dz
     end
 
     # Solve equation system using Thomas algorithm

src/simulation/heavy_species.jl

@@ -1,26 +1,26 @@
 
-function update_heavy_species!(dU, U, params, t)
-    ####################################################################
-    #extract some useful stuff from params
+function update_heavy_species!(dU, U, params, t; apply_boundary_conditions = true)
 
     (;index, Δz_cell, config, cache) = params
     (;
         source_neutrals, source_ion_continuity, source_ion_momentum,
         ncharge, ion_wall_losses
     ) = config
 
-    (;F, UL, UR) = cache
+    (;F, UL, UR, channel_area, dA_dz) = cache
 
     ncells = size(U, 2)
 
-    compute_fluxes!(F, UL, UR, U, params)
+    compute_fluxes!(F, UL, UR, U, params; apply_boundary_conditions)
 
     @inbounds for i in 2:ncells-1
         left = left_edge(i)
         right = right_edge(i)
 
         Δz = Δz_cell[i]
 
+        dlnA_dz = dA_dz[i] / channel_area[i]
+
         # Handle neutrals
         for j in 1:params.num_neutral_fluids
             # Neutral fluxes
@@ -33,8 +33,11 @@ function update_heavy_species!(dU, U, params, t)
         # Handle ions
         for Z in 1:ncharge
             # Ion fluxes
-            dU[index.ρi[Z]  , i] = (F[index.ρi[Z],   left] - F[index.ρi[Z],   right]) / Δz
-            dU[index.ρiui[Z], i] = (F[index.ρiui[Z], left] - F[index.ρiui[Z], right]) / Δz
+            # ∂ρ/∂t + ∂/∂z(ρu) = Q - ρu * ∂/∂z(lnA)
+            ρi = U[index.ρi[Z], i]
+            ρiui = U[index.ρiui[Z], i]
+            dU[index.ρi[Z]  , i] = (F[index.ρi[Z],   left] - F[index.ρi[Z],   right]) / Δz - ρiui * dlnA_dz
+            dU[index.ρiui[Z], i] = (F[index.ρiui[Z], left] - F[index.ρiui[Z], right]) / Δz - ρiui^2 / ρi * dlnA_dz
 
             # User-provided ion source terms
             dU[index.ρi[Z],   i] += source_ion_continuity[Z](U, params, i)