From 9dbce9f58fe0aac7ba05069b3489bd55710d4b53 Mon Sep 17 00:00:00 2001 From: Walter Dal'Maz Silva Date: Sat, 30 Mar 2024 17:12:35 +0100 Subject: [PATCH] Adding DryGranular --- Manifest.toml | 2 +- Project.toml | 1 + docs/Manifest.toml | 2 +- docs/Project.toml | 1 + docs/make.jl | 3 + docs/src/DryGranular/index.md | 220 +++++++++++++++++++ docs/src/DryGranular/kramers.md | 376 ++++++++++++++++++++++++++++++++ src/DryGranular.jl | 374 +++++++++++++++++++++++++++++++ 8 files changed, 977 insertions(+), 2 deletions(-) create mode 100644 docs/src/DryGranular/index.md create mode 100644 docs/src/DryGranular/kramers.md create mode 100644 src/DryGranular.jl diff --git a/Manifest.toml b/Manifest.toml index 4216e3937..59e6ad798 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -2,7 +2,7 @@ julia_version = "1.10.1" manifest_format = "2.0" -project_hash = "63593c2b584fbc32a0e5777ac3429c279bcc4c17" +project_hash = "aa40df391343477c67256f2b12be80ad030257dd" [[deps.ADTypes]] git-tree-sha1 = "016833eb52ba2d6bea9fcb50ca295980e728ee24" diff --git a/Project.toml b/Project.toml index c30ea6b15..5e994b11b 100644 --- a/Project.toml +++ b/Project.toml @@ -4,6 +4,7 @@ authors = ["Walter Dal'Maz Silva "] version = "0.1.0" [deps] +BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" diff --git a/docs/Manifest.toml b/docs/Manifest.toml index 6965551d6..d571724eb 100644 --- a/docs/Manifest.toml +++ b/docs/Manifest.toml @@ -2,7 +2,7 @@ julia_version = "1.10.1" manifest_format = "2.0" -project_hash = "ae56c51e2accfbdfa2d99de62643dc58ae6de2bf" +project_hash = "86498443abd203753402feebd36f4237b78f6ae6" [[deps.ADTypes]] git-tree-sha1 = "016833eb52ba2d6bea9fcb50ca295980e728ee24" diff --git a/docs/Project.toml b/docs/Project.toml index 18b3fe50f..0379da6fd 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,4 +1,5 @@ [deps] +BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" diff --git a/docs/make.jl b/docs/make.jl index a757f1c1b..685181c6c 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -20,6 +20,7 @@ using WallyToolbox using DryAbstract using DryConstants using DryFlowsheet +using DryGranular using DryMaterials using DryUtilities using OpenFOAM @@ -40,6 +41,7 @@ modules = [ DryAbstract, DryConstants, DryFlowsheet, + DryGranular, DryMaterials, DryUtilities, OpenFOAM, @@ -73,6 +75,7 @@ pages = [ "Dry Packages" => [ "DryFlowsheet" => "DryFlowsheet/index.md", + "DryGranular" => "DryGranular/index.md", "DryMaterials" => "DryMaterials/index.md", "Helpers" => "helpers.md" ], diff --git a/docs/src/DryGranular/index.md b/docs/src/DryGranular/index.md new file mode 100644 index 000000000..be6a9f8b1 --- /dev/null +++ b/docs/src/DryGranular/index.md @@ -0,0 +1,220 @@ +# DryGranular + +```@meta +CurrentModule = DryGranular +DocTestSetup = quote + using Statistics + using DryGranular +end +``` + +## General porous media + +Modeling of geometrical characteristics of porous beds is required for including +both their thermal effect or role over chemistry in chemical reactors. A +classical approach used in several commercial and open source tools is that of +Gunn [Gunn1978](@cite). In what follows we develop the ideas that lead to an +analogous model which is implemented by this structure. + +To build the model we will assume a reactor of constant rectangular +cross-section ``{A}_{r}={b}{w}`` and volume ``{V}_{R}={b}{w}{h}``. Its +cross-section perimeter is then ``{P}_{R}=2({b}+{w})``. Inside this reactor we +randomly pack cubic particles ``\beta`` of surface area +``{A}_{\beta}=6{l}_{\beta}^2`` and volume ``{V}_{\beta}={l}_{\beta}^3`` at a +porosity level ``{\phi}``. Thus the total volume of solids inside the reactor is +``{V}_{S}=(1-{\phi}){V}_{R}`` and the approximate number of particles +``{N}=\frac{{V}_{S}}{{V}_{\beta}}``. Following a similar reasoning the total +surface area of particles is ``{A}_{S}={N}{A}_{\beta}``. Performing all the +substitutions so far one finds the following expression + +```math +{A}_{S}=\frac{6(1-{\phi}){b}{w}{h}}{{l}_{\beta}} +``` + +Since the differential ``d{A}={P}d{l}`` holds for the surface of a body over its +length ``{l}``, one can divide the above expression by the reactor length to get +the perimeter of particles in a cross-section. We can further divide by the +cross-section area itself and find the *perimeter density* which is a more +general result, and find the expression proposed by Gunn [Gunn1978](@cite). This +result is summarized in the next equation where the subscript of particle size +was dropped for generality. + +```math +{P} = \frac{6(1-{\phi})}{{l}} +``` + +An estimator of the number of channels per unit cross-section of reactor ``{N}`` +can be related to the porosity through ``{N}\pi{R}^2={\phi}``. Because the above +perimeter is shared between the fluid volume and solids, it holds that +``{N}2\pi{R}=P``. Using these expressions one can solve for the porosity +channels characteristic *radius* ``{R}`` as given below, which is also a result +reported by Gunn [Gunn1978](@cite). + +```math +{R}=\frac{{\phi}{l}}{3(1-{\phi})} +``` + +This model is probided in [`PackedBedPorosityDescriptor`](@ref). + +```@docs +DryGranular.PackedBedPorosityDescriptor +``` + +[`PackedBedPorosityDescriptor`](@ref) can be used to describe the geometry of +exchange section of a packed bed for a single set of arguments. + +```jldoctest +julia> PackedBedPorosityDescriptor(; ϕ = 0.65, l = 0.10, area = 1.0) +PackedBedPorosityDescriptor(P = 21.000000 m, D = 0.123810 m) +``` + +It can also be used to describe randomly varying reactors, what is a more +realistic thing to do when using this structure to simulate real world systems. + +```jldoctest +julia> PackedBedPorosityDescriptor(; + ϕ = 0.65, l = 0.10, + σϕ = 0.03, σl = 0.01, + N = 2, + ϕlims = (0.4, 0.8), + llims = (0.0, 0.3), + seed = 42, + area = 1.0 + ) +PackedBedPorosityDescriptor( + P from 21.455749 m to 24.370742 m + D from 0.125589 m to 0.102353 m +) +``` + +## Rotary kiln models + +In a rotary kiln as proposed by Kramers and Croockewite (1952) +[Kramers1952](@cite). Its goal is to be used as a process support tool or to +integrate more complex models requiring integration of the bed profile. In its +classical statement, the bed height profile ``h(z)`` can be evaluated from +*volume* of flowing material conservation through the following equations. +Coordinate ``z=0`` represents the discharge position where initial condition +must be applied. This is given by the dam height, if any, or particle size. + +```math +\begin{aligned} +\dfrac{dh}{dz} &= C₁ \left[\frac{h}{R}\left(2 - \frac{h}{R}\right)\right]^{-\frac{3}{2}} - C₂\\[6pt] +C₁ &= \frac{3}{4}\dfrac{Φ\tan{γ}}{π R^3 ω}\\[6pt] +C₂ &= \dfrac{\tan{β}}{\cos{γ}} +\end{aligned} +``` + +The structure [`SymbolicLinearKramersModel`](@ref) implements the Kramers' +ordinary differential equation for prediction of bed height profile in a rotary +kiln. This equation is implemented under the formalism of `ModelingToolkit`. + +```@docs +DryGranular.SymbolicLinearKramersModel +``` + + +For integration of this model we implement [`RotaryKilnBedSolution`](@ref). It +provides the solved description of a rotary kiln bed geometry computed from the +solution of bed height along the kiln length. The main goal of the quantities +computed here is their use with heat and mass transfer models for the simulation +of rotary kiln process. A simple post-processing utilitiy +[`plotlinearkramersmodel`](@ref) is also provided. + +```@docs +DryGranular.RotaryKilnBedSolution +DryGranular.plotlinearkramersmodel +``` + +Data in next example is an SI conversion of an example from Kramers and +Croockewite (1952) [Kramers1952](@cite). + +```jldoctest +julia> L = 13.715999999999998; # Kiln length [m] + +julia> D = 1.8897599999999999; # Kiln diameter [m] + +julia> β = 2.3859440303888126; # Kiln slope [°] + +julia> γ = 45.0; # Repose angle [°] + +julia> d = 1.0; # Particle/dam size [mm] + +julia> Φ = 10.363965852671996; # Feed rate [m³/h] + +julia> ω = 3.0300000000000002; # Rotation rate [rev/min] + +julia> bed = RotaryKilnBedSolution(; + model = SymbolicLinearKramersModel(), + L = L, + R = D / 2.0, + Φ = Φ / 3600.0, + ω = ω / 60.0, + β = deg2rad(β), + γ = deg2rad(γ), + d = d / 1000.0 + ); + +julia> bed +RotaryKilnBedSolution(τ = 13.169938 min, ηₘ = 5.913271 %) + +julia> bed.τ +790.1963002204403 +``` + +In the following dummy example we force a very thick *analytical* bed solution, +filling the radius of the rotary drum. + +```jldoctest dummy-1 +julia> R = 1.0e+00; + +julia> Φ = 1.0e-02; + +julia> z = collect(0.0:0.1:10.0); + +julia> h = R * ones(size(z)); + +julia> Aₐ = π * R^2 / 2; + +julia> Vₐ = Aₐ * z[end]; + +julia> bed = RotaryKilnBedSolution(z, h, 0, R, Φ) +RotaryKilnBedSolution(τ = 26.179939 min, ηₘ = 50.000000 %) +``` + +Next we confirm the *internal* evaluations of the model match the expected *analytical* values. + +```jldoctest dummy-1 +julia> mean(bed.θ) ≈ π +true + +julia> mean(bed.l) ≈ 2R +true + +julia> mean(bed.A) ≈ Aₐ +true + +julia> mean(bed.η) ≈ 0.5 +true + +julia> bed.ηₘ ≈ 50.0 +true + +julia> bed.V ≈ Vₐ +true + +julia> bed.τ ≈ Vₐ / Φ +true +``` + +Validation of Kramers' model is provided [here](kramers.md). + +Finally a set of basic equations provided for process analysis. + +```@docs +DryGranular.sullivansηₘ +DryGranular.dimlessNΦ +DryGranular.dimlessNₖ +DryGranular.perrayresidence +DryGranular.kramersnlapprox +``` diff --git a/docs/src/DryGranular/kramers.md b/docs/src/DryGranular/kramers.md new file mode 100644 index 000000000..044f730c4 --- /dev/null +++ b/docs/src/DryGranular/kramers.md @@ -0,0 +1,376 @@ +# Kramers' model + +```@setup kramers +using CairoMakie +using CSV +using DataFrames +using Latexify +using Printf +using DryGranular + +"Partial data from Kramers (1952) Table 3" +const DATA_TABLE3 = """\ +ρ,γ,tan(β),n,ṁ,prod_dimless,η̄ᵣ,hold_real +1480.0,36.0,0.0094,0.059,5.15e-03,18.3,0.111,8.10 +1480.0,36.0,0.0094,0.090,2.68e-03,6.25,0.054,5.00 +1480.0,36.0,0.0094,0.195,1.32e-02,14.2,0.088,7.75 +1480.0,36.0,0.0094,0.232,7.24e-03,6.55,0.043,3.85 +1480.0,36.0,0.0100,0.040,6.38e-03,29.7,0.169,13.3 +1480.0,36.0,0.0100,0.040,5.00e-03,23.2,0.144,11.2 +1480.0,36.0,0.0100,0.069,9.20e-03,24.8,0.150,10.6 +1480.0,36.0,0.0100,0.069,6.53e-03,17.6,0.113,8.50 +1480.0,36.0,0.0100,0.106,1.50e-02,27.8,0.162,12.2 +1480.0,36.0,0.0100,0.159,1.20e-02,14.0,0.092,7.49 +1480.0,36.0,0.0100,0.238,1.55e-02,12.1,0.083,7.48 +1480.0,36.0,0.0100,0.238,1.19e-02,9.22,0.068,6.13 +""" + +"Compares approximate analytical to numerical solution." +function solvekiln(; L, D, Φ, ω, β, γ, d, show = true) + model = RotaryKilnBedSolution(; + model = SymbolicLinearKramersModel(), + L = L, + R = D / 2.0, + Φ = Φ / 3600.0, + ω = ω / 60.0, + β = deg2rad(β), + γ = deg2rad(γ), + d = d / 1000.0 + ) + + optim = kramersnlapprox(; + z = model.z, + R = D / 2.0, + Φ = Φ / 3600.0, + ω = ω / 60.0, + β = deg2rad(β), + γ = deg2rad(γ), + d = d / 1000.0 + ) + + f = nothing + ax = nothing + + if show + f = Figure() + ax = Axis(f[1, 1]) + + lines!(ax, 100model.z/L, 100model.h, linewidth = 3, label = "Numerical") + lines!(ax, 100optim.z/L, 100optim.h, linewidth = 3, label = "Analytical") + + a = @sprintf("%.1f", model.ηₘ) + b = @sprintf("%.1f", optim.ηₘ) + title = "Loading: $(a)% (numerical) | $(b)% (analytical)" + + ax.title = title + ax.xlabel = "Coordinate [%]" + ax.ylabel = "Bed height [cm]" + ax.xticks = 0.0:20.0:100.0 + xlims!(ax, extrema(ax.xticks.val)) + end + + return model, optim, f, ax +end + +"Reference case for alumina kiln testing." +function aluminakiln(ṁ, ω; show = false) + # Density of bed [kg/m³] + ρ = 800.0 + L = 34.0 + D = 1.5 + β = atan(0.025) + + model, optim, f, ax = solvekiln( + L = L, + D = D, + Φ = (1000// 24) * ṁ / ρ, + ω = ω, + β = rad2deg(β), + γ = 33.0, + d = 0.050, + show = show + ) + + τₚ = perrayresidence(L, ω, D, β) + + return model, optim, f, ax, τₚ +end + +"Run `aluminakiln` against some known conditions." +function scanaluminakiln() + ṁlist = [33.6, 43.2] + ωlist = [0.85, 1.20] + + df = DataFrame( + ṁ = Float64[], + ω = Float64[], + η̄ = Float64[], + τᵢ = Float64[], + τₚ = Float64[] + ) + + for ṁ ∈ ṁlist, ω ∈ ωlist + model, _, _, _, τ = aluminakiln(ṁ, ω, show = false) + η̄ = round(model.ηₘ, digits = 0) + τᵢ = round(model.τ / 60.0, digits = 0) + τₚ = round(τ, digits = 0) + push!(df, [ṁ ω η̄ τᵢ τₚ]) + end + + return df +end + +let + # @info("Solution of reference case") + + in1_to_m1(v) = 0.0254 * v + ft1_to_m1(v) = in1_to_m1(12.0) * v + ft3_to_m3(v) = ft1_to_m1(1.0)^3 * v + + # Kiln length [m] + L = ft1_to_m1(45.0) + + # Kiln diameter [m] + D = 2 * ft1_to_m1(3.1) + + # Volume flow rate [m³/h] + Φ = ft3_to_m3(6.1) * 60 + + # Rotation rate (+0.0005) [rev/min] + ω = 0.0505 * 60.0 + + # Kiln slope (0.5in/ft) [°] + β = rad2deg(atan(0.5 / 12)) + + # Repose angle [°] + γ = 45.0 + + # Particle size [mm] + d = 0.050 + + # Conversions to match model inputs. + R = D / 2.0 + Φ = Φ / 3600.0 + ω = ω / 60.0 + β = deg2rad(β) + γ = deg2rad(γ) + d = d / 1000.0 + + # Create problem container. + kramers = RotaryKilnBedSolution(; + model = SymbolicLinearKramersModel(), + L = L, + R = R, + Φ = Φ, + ω = ω, + β = β, + γ = γ, + d = d + ) + + optim = kramersnlapprox(; + z = kramers.z, + R = R, + Φ = Φ, + ω = ω, + β = β, + γ = γ, + d = d + ) + + global kramers_NΦ = dimlessNΦ(R, β, ω, Φ, γ) + global kramers_Nₖ = dimlessNₖ(L, R, β, γ) + global kramers_η̄ₛ = sullivansηₘ(R, β, ω, Φ, γ) + global kramers_ref = kramers + global optim_ref = optim + + global RESULTS_TABLE = DataFrame( + Quantity = [ + "NΦ", + "Nₖ", + "η̄ᵣ", + "η̄ᵢ", + ], + Reference = [ + "1.15", + "1.17", + "5.65", + @sprintf("%.2f", optim_ref.ηₘ) + ], + Computed = [ + @sprintf("%.2f", kramers_NΦ), + @sprintf("%.2f", kramers_Nₖ), + @sprintf("%.2f", kramers_η̄ₛ), + @sprintf("%.2f", kramers_ref.ηₘ) + ] + ) +end + +const TABLE3 = let + # @info("Verification of *Table 3*") + + Dₖ = 0.197 + Lₖ = 1.780 + dₖ = 0.0012 + + table3 = DataFrame(CSV.File(IOBuffer(DATA_TABLE3))) + table3[!, "η̄ᵢ"] = zeros(length(table3[!, "η̄ᵣ"])) + table3[!, "η̄ᵣ"] *= 100 + + model = SymbolicLinearKramersModel() + + for (i, row) in enumerate(eachrow(table3)) + Φ = 3600.0 * row["ṁ"] / row["ρ"] + ω = row["n"] * 60.0 + β = rad2deg(atan(row["tan(β)"])) + γ = row["γ"] + + kramers = RotaryKilnBedSolution(; + model = model, + L = Lₖ, + R = Dₖ / 2.0, + Φ = Φ / 3600.0, + ω = ω / 60.0, + β = deg2rad(β), + γ = deg2rad(γ), + d = dₖ / 1000.0 + ) + + table3[i, "η̄ᵢ"] = round(kramers.ηₘ, digits = 1) + end + + exclude = ["ρ", "γ", "prod_dimless", "hold_real"] + select(table3, Not(exclude)) +end + +const DIMLESSPLOT = let + @info("Dimensionless profiles solution") + + ρ = 1480.0 + L = 20.0 + D = 0.197 + Φ = 5.15e-03 / ρ * 3600 + ω = 0.059 * 60 + β = rad2deg(atan(0.0094)) + γ = 36.0 + + # Conversions to match model inputs. + R = D / 2.0 + Φ = Φ / 3600.0 + ω = ω / 60.0 + β = deg2rad(β) + γ = deg2rad(γ) + + # Things held constant in loop. + NΦ = dimlessNΦ(R, β, ω, Φ, γ) + Nₖ = dimlessNₖ(L, R, β, γ) + model = SymbolicLinearKramersModel() + + f = Figure() + ax = Axis(f[1, 1]) + + for d in [0.05, 0.10, 0.15, 0.193, 0.25] + kramers = RotaryKilnBedSolution(; + model = model, + L = L, + R = R, + Φ = Φ, + ω = ω, + β = β, + γ = γ, + d = d * R * NΦ + ) + + # Dimensionless axes. + z = kramers.z + h = kramers.h / (R * NΦ) + z = @. (L - z) / L * 1 / (NΦ * Nₖ) + z = @. z[1] - z + + label = @sprintf("%.3f", d) + lines!(ax, z, h; linewidth = 2, label = label) + end + + ax.title = "Dimensionless loading curves" + ax.xlabel = "Coordinate" + ax.ylabel = "Bed height" + ax.xticks.val = 0.0:0.1:0.5 + ax.yticks.val = 0.05:0.05:0.25 + xlims!(ax, extrema(ax.xticks.val)) + ylims!(ax, extrema(ax.yticks.val)) + axislegend(ax; position = :rb) + + f +end +``` + +## Sample reference case + +Here we make use of the current implementation to check if it correctly approximates the last example provided in reference paper from [Kramers1952](@cite). To minimize rounding errors causes by unit conversions, we provide the required functions to convert from imperial to international system in the solution process. + +The next table summarizes the results. It is seen that the dimensionless numbers are well approximated. It must be emphasized that the reference estimates η̄ᵣ by a graphical method -- it was 1952 -- and the current value is considered a good enough approximation. Additionally, the equation was not integrated numerically as done here, but engineering relationships were used in the approximation. That said, the proper loading to be considered in our days is η̄ᵢ. + +```@example kramers +mdtable(RESULTS_TABLE, latex=false) # hide +``` + +**Note:** the last value in column `Reference` above is not provided in Kramers' paper but computed from the approximate analytical solution provided by the authors. As we see here, it may get >20% error under some circumstances. + +## Verification of *Table 3* + +In the next cell we provide the kiln dimensions used by Kramers (1952) to experimentally validate the model. Some data from their Tab. 3 is then loaded and all rows are simulated with current model. Fractional hold-up seems to be well correlated at least to a few percent of the reference value. + +```@example kramers +mdtable(TABLE3, latex=false) # hide +``` + +## Dimensionless profiles + +Next step in validation is to check profiles in dimensionless format, as done by Kramers in their Fig. 3. Notice that here we used the numerical integration curves instead of the analytical approximation of profiles, so reproduction and consequences of results are not exactly the same. + +```@example kramers +DIMLESSPLOT # hide +``` + +## Comparison with analytical + +The final step in model validation is to compare the approximate analytical solution proposed by Kramers and the results of numerical integration. It is worth mentioning that numerical integration remains the recommended method because one does not need to verify the ranges of validity of analytical approximation for every use case. + +```@example kramers +let # hide + _, _, f, ax = solvekiln( # hide + L = 10.0, # hide + D = 1.0, # hide + Φ = 1.0, # hide + ω = 1.0, # hide + β = 3.0, # hide + γ = 45.0, # hide + d = 0.001 # hide + ) # hide + # hide + ax.yticks = 0.0:4.0:20.0 # hide + ylims!(ax, extrema(ax.yticks.val)) # hide + f # hide +end # hide +``` + +## Industrial cases + +The following illustrates a practical use case of the model. Next we scan a parameter space to confirm once again the model suitability as an alternative to analytical engineering estimations as per Peray's notebook. + +```@example kramers +let # hide + ṁ = 33.6 # hide + ω = 0.85 # hide + _, _, f, ax, _ = aluminakiln(ṁ, ω, show = true) # hide + ax.yticks = 0.0:6.0:30.0 # hide + ylims!(ax, extrema(ax.yticks.val)) # hide + f # hide +end # hide +``` + +The following table confirms the expected values as per Peray. + +```@example kramers +mdtable(scanaluminakiln(), latex=false) # hide +``` diff --git a/src/DryGranular.jl b/src/DryGranular.jl new file mode 100644 index 000000000..8e68d6d7f --- /dev/null +++ b/src/DryGranular.jl @@ -0,0 +1,374 @@ +module DryGranular + +using CairoMakie +using DifferentialEquations: ODEProblem, Tsit5 +using DifferentialEquations: solve +using Distributions +using DocStringExtensions: TYPEDFIELDS +using Ipopt +using ModelingToolkit +using Printf +using Random +using Trapz: trapz +import JuMP + +export PackedBedPorosityDescriptor +export SymbolicLinearKramersModel +export RotaryKilnBedSolution +export solvelinearkramersmodel +export plotlinearkramersmodel +export dimlessNΦ +export dimlessNₖ +export sullivansηₘ +export perrayresidence +export kramersnlapprox + +""" +Provides description of porosity parameters with stochastic behavior. + +$(TYPEDFIELDS) +""" +struct PackedBedPorosityDescriptor + "Porosity volume fraction in medium [-]." + ϕ::Union{Float64, Vector{Float64}} + + "Characteristic particle size in medium [m]." + l::Union{Float64, Vector{Float64}} + + "Optional standard deviation of porosity volume fraction [-]." + σϕ::Union{Float64, Nothing} + + "Optional standard deviation of characteristic particle size [m]." + σl::Union{Float64, Nothing} + + "Perimeter in reactor cross-section [m]." + P::Union{Float64, Vector{Float64}} + + "Characteristic diameter of porosity channels [m]." + D::Union{Float64, Vector{Float64}} + + "Reactor area used for scaling perimeter [m²]." + A::Float64 + + function PackedBedPorosityDescriptor(; + ϕ::Float64, + l::Float64, + σϕ::Union{Float64, Nothing} = nothing, + σl::Union{Float64, Nothing} = nothing, + N::Union{Int64, Nothing} = nothing, + ϕlims::Tuple{Float64, Float64} = (0.4, 0.8), + llims::Tuple{Float64, Float64} = (0.0, 0.3), + seed::Int64 = 42, + area::Float64 = 1.0 + ) + if !any(isnothing, [σϕ, σl, N]) + Random.seed!(seed) + ϕᵤ = rand(truncated(Normal(ϕ, σϕ), ϕlims...), N) + lᵤ = rand(truncated(Normal(l, σl), llims...), N) + else + ϕᵤ = ϕ + lᵤ = l + end + + P = @. 6.0 * (1.0 - ϕᵤ) / lᵤ + D = @. 4.0 * ϕᵤ / P + + return new(ϕᵤ, lᵤ, σϕ, σl, area * P, D, area) + end +end + +function Base.length(p::PackedBedPorosityDescriptor) + return length(p.ϕ) +end + +function Base.show(io::IO, obj::PackedBedPorosityDescriptor) + if any(isnothing, [obj.σϕ, obj.σl]) + P = @sprintf("%.6f m", obj.P) + D = @sprintf("%.6f m", obj.D) + print(io, "PackedBedPorosityDescriptor(P = $(P), D = $(D))") + else + P = map(x->@sprintf("%10.6f m", x), obj.P) + D = map(x->@sprintf("%10.6f m", x), obj.D) + print(io, """\ + PackedBedPorosityDescriptor( + P from $(P[1]) to $(P[2]) + D from $(D[1]) to $(D[2]) + )\ + """) + end +end + +""" +Creates a reusable linear Kramers model for rotary kiln simulation. + +$(TYPEDFIELDS) +""" +struct SymbolicLinearKramersModel + + "Symbolic kiln internal radius" + R::Num + + "Symbolic kiln feed rate" + Φ::Num + + "Symbolic kiln rotation rate" + ω::Num + + "Symbolic kiln slope" + β::Num + + "Symbolic solids repose angle" + γ::Num + + "Symbolic kiln axial coordinates" + z::Num + + "Symbolic bed height profile" + h::Num + + "Problem ordinary differential equation" + sys::ODESystem + + """ Symbolic model constructor. """ + function SymbolicLinearKramersModel() + # Declare symbols and unknowns. + ps = @parameters R Φ ω β γ + @variables z + @variables h(z) + + # Declare a derivative. + D = Differential(z) + + # Compose problem right-hand side. + C = (3//4) * tan(γ) * Φ / (π * R^3 * ω) + f = C * ((h / R) * (2 - h / R))^(-3//2) + + # *Stack* equation. + eqs = D(h) ~ f - tan(β) / cos(γ) + + # Assembly system for solution. + @named sys = ODESystem(eqs, z, [h], ps) + sys = structural_simplify(sys) + + return new(R, Φ, ω, β, γ, z, h, sys) + end +end + +""" +General geometric description of a bed from Kramers equation solution. + +$(TYPEDFIELDS) + +# Arguments + +Internal elements are initialized through the following constructor: + +```julia +RotaryKilnBedSolution(z, h, β, R, Φ) +``` + +Where parameters are given as: + +- `z`: solution coordinates over length, [m]. +- `h`: bed profile solution over length, [m]. +- `R`: kiln internal radius, [m]. +- `Φ`: kiln feed rate, [m³/s]. + +An outer constructor is also provided for managing the integration of an +instance of `SymbolicLinearKramersModel`. This is the recommended usage +that is illustrated below. + +**Important:** inputs must be provided in international system (SI) units +as a better physical practice. The only exception is the rotation rate `ω` +provided in revolution multiples. If the discharge end is held by a dam, +its height must be provided instead of the particle size, as it is used +as the ODE initial condition. + +- `model`: a symbolic kiln model. +- `L`: kiln length, [m]. +- `R`: kiln internal radius, [m]. +- `Φ`: kiln feed rate, [m³/s]. +- `ω`: kiln rotation rate, [rev/s]. +- `β`: kiln slope, [rad]. +- `γ`: solids repose angle, [rad]. +- `d`: particle size or dam height, [m]. +- `solver`: Solver for `DifferentialEquations`. Defaults to `Tsit5`. +- `rtol`: Relative integration tolerance. Defaults to 1.0e-08. +- `atol`: Absolute integration tolerance. Defaults to 1.0e-08. +""" +struct RotaryKilnBedSolution + "Solution coordinates [m]" + z::Vector{Float64} + + "Solution bed height [m]" + h::Vector{Float64} + + "View angle from kiln center [rad]" + θ::Vector{Float64} + + "Bed-freeboard cord length [m]" + l::Vector{Float64} + + "Local bed cross section area [m²]" + A::Vector{Float64} + + "Local loading based on height [-]" + η::Vector{Float64} + + "Mean loading of kiln [%]" + ηₘ::Float64 + + "Bed integral volume [m³]" + V::Float64 + + "Residence time of particles" + τ::Float64 + + "Kiln slope [rad]" + β::Float64 + + function RotaryKilnBedSolution(z, h, β, R, Φ) + L = z[end] + θ = @. 2acos(1 - h / R) + l = @. 2R * sin(θ / 2) + A = @. (θ * R^2 - l * (R - h)) / 2 + η = @. (θ - sin(θ)) / 2π + ηₘ = 100trapz(z, η) / L + + # Integrate mid-point volume approximation. + Aₘ = (1//2) * (A[1:end-1] + A[2:end]) + δz = z[2:end] - z[1:end-1] + V = sum(@. Aₘ * δz) + + # Residence time is bed volume divided by flow rate. + τ = V / Φ + + # Construct solution object. + return new(z, h, θ, l, A, η, ηₘ, V, τ, β) + end +end + +function RotaryKilnBedSolution(; + model::SymbolicLinearKramersModel, + L::Float64, + R::Float64, + Φ::Float64, + ω::Float64, + β::Float64, + γ::Float64, + d::Float64, + solver::Any = Tsit5(), + rtol::Float64 = 1.0e-08, + atol::Float64 = 1.0e-08 + ) + h = [model.h => d] + p = [model.R => R, + model.Φ => Φ, + model.ω => ω, + model.β => β, + model.γ => γ] + + prob = ODEProblem(model.sys, h, (0.0, L), p, jac = true) + sol = solve(prob, solver, reltol = rtol, abstol = atol) + return RotaryKilnBedSolution(sol.t, sol[1, :], β, R, Φ) +end + +function Base.show(io::IO, obj::RotaryKilnBedSolution) + τ = @sprintf("%.6f min", obj.τ/60) + ηₘ = @sprintf("%.6f", obj.ηₘ) + print(io, "RotaryKilnBedSolution(τ = $(τ), ηₘ = $(ηₘ) %)") +end + +""" + plotlinearkramersmodel( + model::RotaryKilnBedSolution; + normz::Bool = false, + normh::Bool = false + )::Figure + +Standardized plotting of `RotaryKilnBedSolution` bed profile. It +supports normalization of axes throught keywords `normz` for axial +coordinate and `normh` for bed depth. +""" +function plotlinearkramersmodel( + model::RotaryKilnBedSolution; + normz::Bool = false, + normh::Bool = false + )::Figure + z = model.z + h = tan(model.β) * z + model.h + + z = normz ? (100z / maximum(z[end])) : z + h = normh ? (100h / maximum(h[end])) : 100h + + unitz = normz ? "%" : "m" + unith = normh ? "%" : "cm" + + η = @sprintf("%.1f", model.ηₘ) + τ = @sprintf("%.0f", model.τ / 60) + + title = "Loading $(η)% | Residence $(τ) min" + xlab = "Coordinate [$(unitz)]" + ylab = "Bed height [$(unith)]" + + xlims = (normz) ? (0.0, 100.0) : (0.0, model.z[end]) + ylims = (normh) ? (0.0, 100.0) : (0.0, round(maximum(h)+1)) + + fig = Figure() + ax = Axis(fig[1, 1], title = title, xlabel = xlab, ylabel = ylab, + xticks = range(xlims..., 6), yticks = range(ylims..., 6)) + lines!(ax, z, h, color = :red, label = "Profile") + limits!(ax, xlims, ylims) + axislegend(position = :lt) + + return fig +end + +"Kramers (1952) dimensionless group NΦ." +function dimlessNΦ(R, β, ω, Φ, γ) + return Φ * sin(γ) / (ω * R^3 * tan(β)) +end + +"Kramers (1952) dimensionless group Nₖ." +function dimlessNₖ(L, R, β, γ) + return R * cos(γ) / (L * tan(β)) +end + +"Sullivans approximation to kiln filling." +function sullivansηₘ(R, β, ω, Φ, γ) + return 3.8 * dimlessNΦ(R, β, ω, Φ, γ) * sqrt(γ) / sin(γ) +end + +"Compute residence time from Peray's equation." +function perrayresidence(L, ω, D, β) + return 0.19 * L / (ω * D * tan(β)) +end + +"Nonlinear formulation of Kramers model approximate solution." +function kramersnlapprox(; z, R, Φ, ω, β, γ, d) + L = z[end] + N = length(z) + + NΦ = dimlessNΦ(R, β, ω, Φ, γ) + Nₖ = dimlessNₖ(L, R, β, γ) + + C₁ = R * NΦ + C₂ = 3C₁ / (4π * 1.24) + C₃ = C₁ / (L * NΦ * Nₖ) + + optim = JuMP.Model(Ipopt.Optimizer) + JuMP.set_silent(optim) + + @JuMP.variable(optim, h[1:N]) + @JuMP.NLconstraint( + optim, + [i = 1:N], + C₂ * log((d - C₂) / (h[i] - C₂)) - C₃ * z[i] - h[i] + d == 0, + ) + @JuMP.NLconstraint(optim, [i = 1:N], h[i] >= 0.0) + JuMP.optimize!(optim) + + return RotaryKilnBedSolution(z, JuMP.value.(h), β, R, Φ) +end + +end # module DryGranular