Merge pull request #566 from CliMA/he/ventilation

Add ventilation_factor function to Common

Author: haakon-e

docs/src/API.md

@@ -244,6 +244,7 @@ Common.Chen2022_vel_coeffs_B2
 Common.Chen2022_vel_coeffs_B4
 Common.volume_sphere_D
 Common.volume_sphere_R
+Common.ventilation_factor
 ```
 
 # Parameters

src/Common.jl

@@ -21,6 +21,7 @@ export Chen2022_exponential_pdf
 export Chen2022_vel_coeffs_B1
 export Chen2022_vel_coeffs_B2
 export Chen2022_vel_coeffs_B4
+export ventilation_factor
 
 """
     G_func(air_props, tps, T, Liquid())
@@ -361,4 +362,30 @@ See also [`volume_sphere_D`](@ref).
 """
 volume_sphere_R(R) = volume_sphere_D(2R)
 
+"""
+    ventilation_factor(vent, aps, v_term)
+
+Returns a function that computes the ventilation factor for a particle as a function of its diameter, `D`.
+
+The ventilation factor parameterizes the increase in the mass and heat exchange for falling particles.
+
+# Arguments
+- `vent`: Ventilation parameterization constants, [`CMP.VentilationSB2005`](@ref)
+- `aps`: Parameters with air properties, [`CMP.AirProperties`](@ref)
+- `v_term`: A function `v_term(D)` that returns the terminal velocity of a particle with diameter `D`
+
+# Returns
+- `F_v(D)`: The ventilation factor as a function of diameter, `D`
+
+See e.g. [SeifertBeheng2006](@cite) Eq. (24) for the definition of the ventilation factor.
+"""
+function ventilation_factor(vent, aps, v_term)
+    (; vent_a, vent_b) = vent
+    (; ν_air, D_vapor) = aps
+    N_sc = ν_air / D_vapor           # Schmidt number
+    N_Re(D) = D * v_term(D) / ν_air  # Reynolds number
+    F_v(D) = vent_a + vent_b * ∛(N_sc) * √(N_Re(D))  # Ventilation factor
+    return F_v
+end
+
 end # module end

test/runtests.jl

@@ -17,3 +17,4 @@ include("aqua.jl")
 include("performance_tests.jl")
 include("aerosol_activation_calibration.jl")
 include("ice_nucleation_calibration.jl")
+include("ventilation_tests.jl")

test/ventilation_tests.jl

@@ -0,0 +1,34 @@
+using Test: @testset, @test, @test_throws, @test_broken
+import CloudMicrophysics.P3Scheme as P3
+import CloudMicrophysics.Common as CO
+import CloudMicrophysics.Parameters as CMP
+import CloudMicrophysics.Microphysics2M as CM2
+import Thermodynamics as TD
+import ClimaParams as CP
+
+@info("Ventilation factor tests")
+
+function test_ventilation_factor(FT)
+    @testset "Ventilation factor / P3 terminal velocity smoke test ($FT)" begin
+        params = CMP.ParametersP3(FT)
+        vel = CMP.Chen2022VelType(FT)
+        aps = CMP.AirProperties(FT)
+        F_rim = FT(0.5)  # Riming fraction [-]
+        ρ_r = FT(500)    # Riming density [kg/m³]
+        ρₐ = FT(1.2)     # Air density [kg/m³]
+        state = P3.get_state(params; F_rim, ρ_r)
+        vent = state.params.vent
+        v_term = D -> P3.ice_particle_terminal_velocity(state, D, vel, ρₐ)
+        vent_factor = CO.ventilation_factor(vent, aps, v_term)
+        Ds = range(FT(0.5e-4), stop = FT(4.5e-4), length = 5)
+        calc_vents = vent_factor.(Ds)
+        smoke_vents = [0.91818553, 1.3191912, 1.7451854, 2.1598392, 2.5553002]
+        for (calc_vent, smoke_vent) in zip(calc_vents, smoke_vents)
+            @test calc_vent ≈ smoke_vent rtol = 1e-6
+        end
+    end
+end
+
+for FT in [Float32, Float64]
+    test_ventilation_factor(FT)
+end