From 73da92c446c6f60fc09ad1cf5e579bbf5eccb34d Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Wed, 27 Mar 2024 15:55:37 +0100 Subject: [PATCH] Sutherlands Law for temperature dependent viscosity (#1808) * sample 2D * bug fixes * typo * sutherlands 1d * sutherlands 3d * fmt * main merge * test var mu * variable mu * fmt * example using sutherland * comment * example * reset toml * Shorten * Unit test * fmt * Update examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex_sutherland.jl * remove todo * Apply suggestions from code review Co-authored-by: Andrew Winters * test vals * non-functionalized mu * comments & dosctrings * docstrings * fix tests * avoid if * docstring and comments * Update src/equations/compressible_navier_stokes.jl Co-authored-by: Michael Schlottke-Lakemper * fmt --------- Co-authored-by: Andrew Winters Co-authored-by: Michael Schlottke-Lakemper Co-authored-by: Jesse Chan <1156048+jlchan@users.noreply.github.com> --- .../elixir_navierstokes_lid_driven_cavity.jl | 5 +- ...elixir_navierstokes_taylor_green_vortex.jl | 5 +- .../elixir_navierstokes_lid_driven_cavity.jl | 5 +- ...ixir_navierstokes_lid_driven_cavity_amr.jl | 5 +- .../elixir_navierstokes_blast_wave_amr.jl | 5 +- ...elixir_navierstokes_taylor_green_vortex.jl | 5 +- ...ir_navierstokes_taylor_green_vortex_amr.jl | 5 +- ...lixir_navierstokes_convergence_periodic.jl | 1 - .../elixir_navierstokes_lid_driven_cavity.jl | 5 +- .../elixir_navierstokes_shearlayer_amr.jl | 5 +- ...elixir_navierstokes_taylor_green_vortex.jl | 5 +- ...erstokes_taylor_green_vortex_sutherland.jl | 94 +++++++++++++++++++ ...elixir_navierstokes_taylor_green_vortex.jl | 5 +- src/equations/compressible_navier_stokes.jl | 14 +++ .../compressible_navier_stokes_1d.jl | 26 +++-- .../compressible_navier_stokes_2d.jl | 26 +++-- .../compressible_navier_stokes_3d.jl | 26 +++-- test/test_parabolic_2d.jl | 34 +++++-- test/test_unit.jl | 39 ++++++++ 19 files changed, 243 insertions(+), 72 deletions(-) create mode 100644 examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex_sutherland.jl diff --git a/examples/dgmulti_2d/elixir_navierstokes_lid_driven_cavity.jl b/examples/dgmulti_2d/elixir_navierstokes_lid_driven_cavity.jl index 7c55cbf0cc..a612dd0e0d 100644 --- a/examples/dgmulti_2d/elixir_navierstokes_lid_driven_cavity.jl +++ b/examples/dgmulti_2d/elixir_navierstokes_lid_driven_cavity.jl @@ -4,12 +4,11 @@ using Trixi ############################################################################### # semidiscretization of the ideal compressible Navier-Stokes equations -# TODO: parabolic; unify names of these accessor functions prandtl_number() = 0.72 -mu() = 0.001 +mu = 0.001 equations = CompressibleEulerEquations2D(1.4) -equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu(), +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu, Prandtl = prandtl_number()) # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux diff --git a/examples/dgmulti_3d/elixir_navierstokes_taylor_green_vortex.jl b/examples/dgmulti_3d/elixir_navierstokes_taylor_green_vortex.jl index dedd8267a3..9ae90ac47b 100644 --- a/examples/dgmulti_3d/elixir_navierstokes_taylor_green_vortex.jl +++ b/examples/dgmulti_3d/elixir_navierstokes_taylor_green_vortex.jl @@ -5,12 +5,11 @@ using Trixi ############################################################################### # semidiscretization of the compressible Navier-Stokes equations -# TODO: parabolic; unify names of these accessor functions prandtl_number() = 0.72 -mu() = 6.25e-4 # equivalent to Re = 1600 +mu = 6.25e-4 # equivalent to Re = 1600 equations = CompressibleEulerEquations3D(1.4) -equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu(), +equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu, Prandtl = prandtl_number()) """ diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl index bc28ae6ffb..728736fe49 100644 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl @@ -4,12 +4,11 @@ using Trixi ############################################################################### # semidiscretization of the ideal compressible Navier-Stokes equations -# TODO: parabolic; unify names of these accessor functions prandtl_number() = 0.72 -mu() = 0.001 +mu = 0.001 equations = CompressibleEulerEquations2D(1.4) -equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu(), +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu, Prandtl = prandtl_number()) # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity_amr.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity_amr.jl index 898366969a..d9eaf4fdb7 100644 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity_amr.jl +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity_amr.jl @@ -4,12 +4,11 @@ using Trixi ############################################################################### # semidiscretization of the ideal compressible Navier-Stokes equations -# TODO: parabolic; unify names of these accessor functions prandtl_number() = 0.72 -mu() = 0.001 +mu = 0.001 equations = CompressibleEulerEquations2D(1.4) -equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu(), +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu, Prandtl = prandtl_number()) # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux diff --git a/examples/p4est_3d_dgsem/elixir_navierstokes_blast_wave_amr.jl b/examples/p4est_3d_dgsem/elixir_navierstokes_blast_wave_amr.jl index 5df89fbcdf..d556d0ab70 100644 --- a/examples/p4est_3d_dgsem/elixir_navierstokes_blast_wave_amr.jl +++ b/examples/p4est_3d_dgsem/elixir_navierstokes_blast_wave_amr.jl @@ -5,12 +5,11 @@ using Trixi ############################################################################### # semidiscretization of the compressible Navier-Stokes equations -# TODO: parabolic; unify names of these accessor functions prandtl_number() = 0.72 -mu() = 6.25e-4 # equivalent to Re = 1600 +mu = 6.25e-4 # equivalent to Re = 1600 equations = CompressibleEulerEquations3D(1.4) -equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu(), +equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu, Prandtl = prandtl_number()) function initial_condition_3d_blast_wave(x, t, equations::CompressibleEulerEquations3D) diff --git a/examples/p4est_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl b/examples/p4est_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl index d785013f5a..9c90e4d321 100644 --- a/examples/p4est_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl +++ b/examples/p4est_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl @@ -5,12 +5,11 @@ using Trixi ############################################################################### # semidiscretization of the compressible Navier-Stokes equations -# TODO: parabolic; unify names of these accessor functions prandtl_number() = 0.72 -mu() = 6.25e-4 # equivalent to Re = 1600 +mu = 6.25e-4 # equivalent to Re = 1600 equations = CompressibleEulerEquations3D(1.4) -equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu(), +equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu, Prandtl = prandtl_number()) """ diff --git a/examples/p4est_3d_dgsem/elixir_navierstokes_taylor_green_vortex_amr.jl b/examples/p4est_3d_dgsem/elixir_navierstokes_taylor_green_vortex_amr.jl index c15227a1c2..2741f0df17 100644 --- a/examples/p4est_3d_dgsem/elixir_navierstokes_taylor_green_vortex_amr.jl +++ b/examples/p4est_3d_dgsem/elixir_navierstokes_taylor_green_vortex_amr.jl @@ -5,12 +5,11 @@ using Trixi ############################################################################### # semidiscretization of the compressible Navier-Stokes equations -# TODO: parabolic; unify names of these accessor functions prandtl_number() = 0.72 -mu() = 6.25e-4 # equivalent to Re = 1600 +mu = 6.25e-4 # equivalent to Re = 1600 equations = CompressibleEulerEquations3D(1.4) -equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu(), +equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu, Prandtl = prandtl_number()) """ diff --git a/examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic.jl b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic.jl index 33ad0d5271..eab0840f38 100644 --- a/examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic.jl +++ b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic.jl @@ -5,7 +5,6 @@ using Trixi ############################################################################### # semidiscretization of the compressible Navier-Stokes equations -# TODO: parabolic; unify names of these accessor functions prandtl_number() = 0.72 mu() = 6.25e-4 # equivalent to Re = 1600 diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl b/examples/tree_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl index 70d76fc907..b8e20e27f6 100644 --- a/examples/tree_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl +++ b/examples/tree_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl @@ -4,12 +4,11 @@ using Trixi ############################################################################### # semidiscretization of the ideal compressible Navier-Stokes equations -# TODO: parabolic; unify names of these accessor functions prandtl_number() = 0.72 -mu() = 0.001 +mu = 0.001 equations = CompressibleEulerEquations2D(1.4) -equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu(), +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu, Prandtl = prandtl_number()) # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl b/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl index 4d92ea261e..a7492bafb4 100644 --- a/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl +++ b/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl @@ -5,12 +5,11 @@ using Trixi ############################################################################### # semidiscretization of the compressible Navier-Stokes equations -# TODO: parabolic; unify names of these accessor functions prandtl_number() = 0.72 -mu() = 1.0 / 3.0 * 10^(-5) # equivalent to Re = 30,000 +mu = 1.0 / 3.0 * 10^(-4) # equivalent to Re = 30,000 equations = CompressibleEulerEquations2D(1.4) -equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu(), +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu, Prandtl = prandtl_number()) """ diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex.jl b/examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex.jl index a3e38bf6d9..c6e5f0bc40 100644 --- a/examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex.jl +++ b/examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex.jl @@ -5,12 +5,11 @@ using Trixi ############################################################################### # semidiscretization of the compressible Navier-Stokes equations -# TODO: parabolic; unify names of these accessor functions prandtl_number() = 0.72 -mu() = 6.25e-4 # equivalent to Re = 1600 +mu = 6.25e-4 # equivalent to Re = 1600 equations = CompressibleEulerEquations2D(1.4) -equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu(), +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu, Prandtl = prandtl_number()) """ diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex_sutherland.jl b/examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex_sutherland.jl new file mode 100644 index 0000000000..9598ae184c --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex_sutherland.jl @@ -0,0 +1,94 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Navier-Stokes equations + +prandtl_number() = 0.72 + +# Use Sutherland's law for a temperature-dependent viscosity. +# For details, see e.g. +# Frank M. White: Viscous Fluid Flow, 2nd Edition. +# 1991, McGraw-Hill, ISBN, 0-07-069712-4 +# Pages 28 and 29. +@inline function mu(u, equations) + T_ref = 291.15 + + R_specific_air = 287.052874 + T = R_specific_air * Trixi.temperature(u, equations) + + C_air = 120.0 + mu_ref_air = 1.827e-5 + + return mu_ref_air * (T_ref + C_air) / (T + C_air) * (T / T_ref)^1.5 +end + +equations = CompressibleEulerEquations2D(1.4) +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu, + Prandtl = prandtl_number()) + +""" + initial_condition_taylor_green_vortex(x, t, equations::CompressibleEulerEquations2D) + +The classical viscous Taylor-Green vortex in 2D. +This forms the basis behind the 3D case found for instance in + - Jonathan R. Bull and Antony Jameson + Simulation of the Compressible Taylor Green Vortex using High-Order Flux Reconstruction Schemes + [DOI: 10.2514/6.2014-3210](https://doi.org/10.2514/6.2014-3210) +""" +function initial_condition_taylor_green_vortex(x, t, + equations::CompressibleEulerEquations2D) + A = 1.0 # magnitude of speed + Ms = 0.1 # maximum Mach number + + rho = 1.0 + v1 = A * sin(x[1]) * cos(x[2]) + v2 = -A * cos(x[1]) * sin(x[2]) + p = (A / Ms)^2 * rho / equations.gamma # scaling to get Ms + p = p + 1.0 / 4.0 * A^2 * rho * (cos(2 * x[1]) + cos(2 * x[2])) + + return prim2cons(SVector(rho, v1, v2, p), equations) +end +initial_condition = initial_condition_taylor_green_vortex + +volume_flux = flux_ranocha +solver = DGSEM(polydeg = 3, surface_flux = flux_hllc, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +coordinates_min = (-1.0, -1.0) .* pi +coordinates_max = (1.0, 1.0) .* pi +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 4, + n_cells_max = 100_000) + +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), + initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 20.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + save_analysis = true, + extra_analysis_integrals = (energy_kinetic, + energy_internal)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback) + +############################################################################### +# run the simulation + +time_int_tol = 1e-9 +sol = solve(ode, RDPK3SpFSAL49(); abstol = time_int_tol, reltol = time_int_tol, + ode_default_options()..., callback = callbacks) +summary_callback() # print the timer summary diff --git a/examples/tree_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl b/examples/tree_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl index 65bd9aa133..3e54c791ec 100644 --- a/examples/tree_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl +++ b/examples/tree_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl @@ -5,12 +5,11 @@ using Trixi ############################################################################### # semidiscretization of the compressible Navier-Stokes equations -# TODO: parabolic; unify names of these accessor functions prandtl_number() = 0.72 -mu() = 6.25e-4 # equivalent to Re = 1600 +mu = 6.25e-4 # equivalent to Re = 1600 equations = CompressibleEulerEquations3D(1.4) -equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu(), +equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu, Prandtl = prandtl_number()) """ diff --git a/src/equations/compressible_navier_stokes.jl b/src/equations/compressible_navier_stokes.jl index 3059771197..c530625f22 100644 --- a/src/equations/compressible_navier_stokes.jl +++ b/src/equations/compressible_navier_stokes.jl @@ -62,3 +62,17 @@ Under `GradientVariablesEntropy`, the Navier-Stokes discretization is provably e """ struct GradientVariablesPrimitive end struct GradientVariablesEntropy end + +""" + dynamic_viscosity(u, equations) + +Wrapper for the dynamic viscosity that calls +`dynamic_viscosity(u, equations.mu, equations)`, which dispatches on the type of +`equations.mu`. +For constant `equations.mu`, i.e., `equations.mu` is of `Real`-type it is returned directly. +In all other cases, `equations.mu` is assumed to be a function with arguments +`u` and `equations` and is called with these arguments. +""" +dynamic_viscosity(u, equations) = dynamic_viscosity(u, equations.mu, equations) +dynamic_viscosity(u, mu::Real, equations) = mu +dynamic_viscosity(u, mu::T, equations) where {T} = mu(u, equations) diff --git a/src/equations/compressible_navier_stokes_1d.jl b/src/equations/compressible_navier_stokes_1d.jl index d2c46ecc7d..3dbdf11c56 100644 --- a/src/equations/compressible_navier_stokes_1d.jl +++ b/src/equations/compressible_navier_stokes_1d.jl @@ -21,6 +21,9 @@ the [`CompressibleEulerEquations1D`](@ref). Fluid properties such as the dynamic viscosity ``\mu`` can be provided in any consistent unit system, e.g., [``\mu``] = kg m⁻¹ s⁻¹. +The viscosity ``\mu`` may be a constant or a function of the current state, e.g., +depending on temperature (Sutherland's law): ``\mu = \mu(T)``. +In the latter case, the function `mu` needs to have the signature `mu(u, equations)`. The particular form of the compressible Navier-Stokes implemented is ```math @@ -80,7 +83,7 @@ where w_2 = \frac{\rho v1}{p},\, w_3 = -\frac{\rho}{p} ``` """ -struct CompressibleNavierStokesDiffusion1D{GradientVariables, RealT <: Real, +struct CompressibleNavierStokesDiffusion1D{GradientVariables, RealT <: Real, Mu, E <: AbstractCompressibleEulerEquations{1}} <: AbstractCompressibleNavierStokesDiffusion{1, 3, GradientVariables} # TODO: parabolic @@ -89,7 +92,7 @@ struct CompressibleNavierStokesDiffusion1D{GradientVariables, RealT <: Real, gamma::RealT # ratio of specific heats inv_gamma_minus_one::RealT # = inv(gamma - 1); can be used to write slow divisions as fast multiplications - mu::RealT # viscosity + mu::Mu # viscosity Pr::RealT # Prandtl number kappa::RealT # thermal diffusivity for Fick's law @@ -103,16 +106,17 @@ function CompressibleNavierStokesDiffusion1D(equations::CompressibleEulerEquatio gradient_variables = GradientVariablesPrimitive()) gamma = equations.gamma inv_gamma_minus_one = equations.inv_gamma_minus_one - μ, Pr = promote(mu, Prandtl) # Under the assumption of constant Prandtl number the thermal conductivity - # constant is kappa = gamma μ / ((gamma-1) Pr). + # constant is kappa = gamma μ / ((gamma-1) Prandtl). # Important note! Factor of μ is accounted for later in `flux`. - kappa = gamma * inv_gamma_minus_one / Pr + # This avoids recomputation of kappa for non-constant μ. + kappa = gamma * inv_gamma_minus_one / Prandtl CompressibleNavierStokesDiffusion1D{typeof(gradient_variables), typeof(gamma), + typeof(mu), typeof(equations)}(gamma, inv_gamma_minus_one, - μ, Pr, kappa, + mu, Prandtl, kappa, equations, gradient_variables) end @@ -159,10 +163,12 @@ function flux(u, gradients, orientation::Integer, # in the implementation q1 = equations.kappa * dTdx - # Constant dynamic viscosity is copied to a variable for readability. - # Offers flexibility for dynamic viscosity via Sutherland's law where it depends - # on temperature and reference values, Ts and Tref such that mu(T) - mu = equations.mu + # In the simplest cases, the user passed in `mu` or `mu()` + # (which returns just a constant) but + # more complex functions like Sutherland's law are possible. + # `dynamic_viscosity` is a helper function that handles both cases + # by dispatching on the type of `equations.mu`. + mu = dynamic_viscosity(u, equations) # viscous flux components in the x-direction f1 = zero(rho) diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index 5df7c01ca5..3256343703 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -21,6 +21,9 @@ the [`CompressibleEulerEquations2D`](@ref). Fluid properties such as the dynamic viscosity ``\mu`` can be provided in any consistent unit system, e.g., [``\mu``] = kg m⁻¹ s⁻¹. +The viscosity ``\mu`` may be a constant or a function of the current state, e.g., +depending on temperature (Sutherland's law): ``\mu = \mu(T)``. +In the latter case, the function `mu` needs to have the signature `mu(u, equations)`. The particular form of the compressible Navier-Stokes implemented is ```math @@ -80,7 +83,7 @@ where w_2 = \frac{\rho v_1}{p},\, w_3 = \frac{\rho v_2}{p},\, w_4 = -\frac{\rho}{p} ``` """ -struct CompressibleNavierStokesDiffusion2D{GradientVariables, RealT <: Real, +struct CompressibleNavierStokesDiffusion2D{GradientVariables, RealT <: Real, Mu, E <: AbstractCompressibleEulerEquations{2}} <: AbstractCompressibleNavierStokesDiffusion{2, 4, GradientVariables} # TODO: parabolic @@ -89,7 +92,7 @@ struct CompressibleNavierStokesDiffusion2D{GradientVariables, RealT <: Real, gamma::RealT # ratio of specific heats inv_gamma_minus_one::RealT # = inv(gamma - 1); can be used to write slow divisions as fast multiplications - mu::RealT # viscosity + mu::Mu # viscosity Pr::RealT # Prandtl number kappa::RealT # thermal diffusivity for Fick's law @@ -103,16 +106,17 @@ function CompressibleNavierStokesDiffusion2D(equations::CompressibleEulerEquatio gradient_variables = GradientVariablesPrimitive()) gamma = equations.gamma inv_gamma_minus_one = equations.inv_gamma_minus_one - μ, Pr = promote(mu, Prandtl) # Under the assumption of constant Prandtl number the thermal conductivity - # constant is kappa = gamma μ / ((gamma-1) Pr). + # constant is kappa = gamma μ / ((gamma-1) Prandtl). # Important note! Factor of μ is accounted for later in `flux`. - kappa = gamma * inv_gamma_minus_one / Pr + # This avoids recomputation of kappa for non-constant μ. + kappa = gamma * inv_gamma_minus_one / Prandtl CompressibleNavierStokesDiffusion2D{typeof(gradient_variables), typeof(gamma), + typeof(mu), typeof(equations)}(gamma, inv_gamma_minus_one, - μ, Pr, kappa, + mu, Prandtl, kappa, equations, gradient_variables) end @@ -168,10 +172,12 @@ function flux(u, gradients, orientation::Integer, q1 = equations.kappa * dTdx q2 = equations.kappa * dTdy - # Constant dynamic viscosity is copied to a variable for readability. - # Offers flexibility for dynamic viscosity via Sutherland's law where it depends - # on temperature and reference values, Ts and Tref such that mu(T) - mu = equations.mu + # In the simplest cases, the user passed in `mu` or `mu()` + # (which returns just a constant) but + # more complex functions like Sutherland's law are possible. + # `dynamic_viscosity` is a helper function that handles both cases + # by dispatching on the type of `equations.mu`. + mu = dynamic_viscosity(u, equations) if orientation == 1 # viscous flux components in the x-direction diff --git a/src/equations/compressible_navier_stokes_3d.jl b/src/equations/compressible_navier_stokes_3d.jl index e5567ae578..9833122eb3 100644 --- a/src/equations/compressible_navier_stokes_3d.jl +++ b/src/equations/compressible_navier_stokes_3d.jl @@ -21,6 +21,9 @@ the [`CompressibleEulerEquations3D`](@ref). Fluid properties such as the dynamic viscosity ``\mu`` can be provided in any consistent unit system, e.g., [``\mu``] = kg m⁻¹ s⁻¹. +The viscosity ``\mu`` may be a constant or a function of the current state, e.g., +depending on temperature (Sutherland's law): ``\mu = \mu(T)``. +In the latter case, the function `mu` needs to have the signature `mu(u, equations)`. The particular form of the compressible Navier-Stokes implemented is ```math @@ -80,7 +83,7 @@ where w_2 = \frac{\rho v_1}{p},\, w_3 = \frac{\rho v_2}{p},\, w_4 = \frac{\rho v_3}{p},\, w_5 = -\frac{\rho}{p} ``` """ -struct CompressibleNavierStokesDiffusion3D{GradientVariables, RealT <: Real, +struct CompressibleNavierStokesDiffusion3D{GradientVariables, RealT <: Real, Mu, E <: AbstractCompressibleEulerEquations{3}} <: AbstractCompressibleNavierStokesDiffusion{3, 5, GradientVariables} # TODO: parabolic @@ -89,7 +92,7 @@ struct CompressibleNavierStokesDiffusion3D{GradientVariables, RealT <: Real, gamma::RealT # ratio of specific heats inv_gamma_minus_one::RealT # = inv(gamma - 1); can be used to write slow divisions as fast multiplications - mu::RealT # viscosity + mu::Mu # viscosity Pr::RealT # Prandtl number kappa::RealT # thermal diffusivity for Fick's law @@ -103,16 +106,17 @@ function CompressibleNavierStokesDiffusion3D(equations::CompressibleEulerEquatio gradient_variables = GradientVariablesPrimitive()) gamma = equations.gamma inv_gamma_minus_one = equations.inv_gamma_minus_one - μ, Pr = promote(mu, Prandtl) # Under the assumption of constant Prandtl number the thermal conductivity - # constant is kappa = gamma μ / ((gamma-1) Pr). + # constant is kappa = gamma μ / ((gamma-1) Prandtl). # Important note! Factor of μ is accounted for later in `flux`. - kappa = gamma * inv_gamma_minus_one / Pr + # This avoids recomputation of kappa for non-constant μ. + kappa = gamma * inv_gamma_minus_one / Prandtl CompressibleNavierStokesDiffusion3D{typeof(gradient_variables), typeof(gamma), + typeof(mu), typeof(equations)}(gamma, inv_gamma_minus_one, - μ, Pr, kappa, + mu, Prandtl, kappa, equations, gradient_variables) end @@ -181,10 +185,12 @@ function flux(u, gradients, orientation::Integer, q2 = equations.kappa * dTdy q3 = equations.kappa * dTdz - # Constant dynamic viscosity is copied to a variable for readability. - # Offers flexibility for dynamic viscosity via Sutherland's law where it depends - # on temperature and reference values, Ts and Tref such that mu(T) - mu = equations.mu + # In the simplest cases, the user passed in `mu` or `mu()` + # (which returns just a constant) but + # more complex functions like Sutherland's law are possible. + # `dynamic_viscosity` is a helper function that handles both cases + # by dispatching on the type of `equations.mu`. + mu = dynamic_viscosity(u, equations) if orientation == 1 # viscous flux components in the x-direction diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 9f1382caa6..d47c34f9e7 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -476,20 +476,38 @@ end @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_navierstokes_shearlayer_amr.jl"), l2=[ - 0.00526017743452336, - 0.4130430692895672, - 0.4310996183791349, - 1.1544344171604635, + 0.005155557460409018, + 0.4048446934219344, + 0.43040068852937047, + 1.1255130552079322, ], linf=[ - 0.03492185879198495, - 1.392635891671335, - 1.357551616406459, - 8.713760873018146, + 0.03287305649809613, + 1.1656793717431393, + 1.3917196016246969, + 8.146587380114653, ], tspan=(0.0, 0.7)) end +@trixi_testset "TreeMesh2D: elixir_navierstokes_taylor_green_vortex_sutherland.jl" begin + @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", + "elixir_navierstokes_taylor_green_vortex_sutherland.jl"), + l2=[ + 0.001452856280034929, + 0.0007538775539989481, + 0.0007538775539988681, + 0.011035506549989587, + ], + linf=[ + 0.003291912841311362, + 0.002986462478096974, + 0.0029864624780958637, + 0.0231954665514138, + ], + tspan=(0.0, 1.0)) +end + @trixi_testset "P4estMesh2D: elixir_advection_diffusion_periodic.jl" begin @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", "elixir_advection_diffusion_periodic.jl"), diff --git a/test/test_unit.jl b/test/test_unit.jl index 03a78f6918..79950f58d5 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -1576,6 +1576,45 @@ end @test mesh.boundary_faces[:entire_boundary] == [1, 2] end + +@testset "Sutherlands Law" begin + function mu(u, equations) + T_ref = 291.15 + + R_specific_air = 287.052874 + T = R_specific_air * Trixi.temperature(u, equations) + + C_air = 120.0 + mu_ref_air = 1.827e-5 + + return mu_ref_air * (T_ref + C_air) / (T + C_air) * (T / T_ref)^1.5 + end + + function mu_control(u, equations, T_ref, R_specific, C, mu_ref) + T = R_specific * Trixi.temperature(u, equations) + + return mu_ref * (T_ref + C) / (T + C) * (T / T_ref)^1.5 + end + + # Dry air (values from Wikipedia: https://de.wikipedia.org/wiki/Sutherland-Modell) + T_ref = 291.15 + C = 120.0 # Sutherland's constant + R_specific = 287.052874 + mu_ref = 1.827e-5 + prandtl_number() = 0.72 + gamma = 1.4 + + equations = CompressibleEulerEquations2D(gamma) + equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu, + Prandtl = prandtl_number()) + + # Flow at rest + u = prim2cons(SVector(1.0, 0.0, 0.0, 1.0), equations_parabolic) + + # Comparison value from https://www.engineeringtoolbox.com/air-absolute-kinematic-viscosity-d_601.html at 18°C + @test isapprox(mu_control(u, equations_parabolic, T_ref, R_specific, C, mu_ref), + 1.803e-5, atol = 5e-8) +end end end #module