From aa9ca7af24f5163976152db7abad0016b2a0ac9a Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Thu, 6 Oct 2022 18:55:47 +0200 Subject: [PATCH 01/12] remove the nondimensionalization and update the docstring for CNS --- src/equations/compressible_euler_2d.jl | 1 - .../compressible_navier_stokes_2d.jl | 81 ++++++------------- 2 files changed, 26 insertions(+), 56 deletions(-) diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index 35852c0d1f9..7bda9999605 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -36,7 +36,6 @@ Here, ``\rho`` is the density, ``v_1``,`v_2` the velocities, ``e`` the specific p = (\gamma - 1) \left( \rho e - \frac{1}{2} \rho (v_1^2+v_2^2) \right) ``` the pressure. - """ struct CompressibleEulerEquations2D{RealT<:Real} <: AbstractCompressibleEulerEquations{2, 4} gamma::RealT # ratio of specific heats diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index f0ba41499f4..4c78add1e19 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -1,5 +1,5 @@ @doc raw""" - CompressibleNavierStokesDiffusion2D(gamma, Re, Pr, Ma_inf, equations, + CompressibleNavierStokesDiffusion2D(gamma, mu, Pr, R, equations, gradient_variables=GradientVariablesPrimitive()) These equations contain the diffusion (i.e. parabolic) terms applied @@ -7,9 +7,9 @@ to mass, momenta, and total energy together with the advective terms from the [`CompressibleEulerEquations2D`](@ref). - `gamma`: adiabatic constant, -- `Re`: Reynolds number, +- `mu`: viscosity, - `Pr`: Prandtl number, -- `Ma_inf`: free-stream Mach number +- `R`: gas constant - `equations`: instance of the [`CompressibleEulerEquations2D`](@ref) - `gradient_variables`: which variables the gradients are taken with respect to. Defaults to `GradientVariablesPrimitive()`. @@ -45,10 +45,10 @@ where ``\underline{I}`` is the ``2\times 2`` identity matrix and the heat flux i \nabla q = -\kappa\nabla\left(T\right),\quad T = \frac{p}{R\rho} ``` where ``T`` is the temperature and ``\kappa`` is the thermal conductivity for Fick's law. -Under the assumption that the gas has a constant Prandtl number +Under the assumption that the gas has a constant Prandtl number, the thermal conductivity is ```math -\kappa = \frac{\gamma \mu R}{(\gamma - 1)\textrm{Pr}} +\kappa = \frac{\gamma \mu R}{(\gamma - 1)\textrm{Pr}}. ``` In two spatial dimensions we require gradients for three quantities, e.g., @@ -65,24 +65,8 @@ where w_2 = \frac{\rho v_1}{p},\, w_3 = \frac{\rho v_2}{p},\, w_4 = -\frac{\rho}{p} ``` -For this particular scaling the vicosity is set internally to be μ = 1/Re. -Further, the nondimensionalization takes the density-temperature-sound speed as -the principle quantities such that -``` -rho_inf = 1.0 -T_ref = 1.0 -c_inf = 1.0 -p_inf = 1.0 / gamma -u_inf = Ma_inf -R = 1.0 / gamma -``` - -Other normalization strategies exist, see the reference below for details. -- Marc Montagnac (2013) - Variable Normalization (nondimensionalization and scaling) for Navier-Stokes - equations: a practical guide - [CERFACS Technical report](https://www.cerfacs.fr/~montagna/TR-CFD-13-77.pdf) -The scaling used herein is Section 4.5 of the reference. +#!!! warning "Experimental code" +# This code is experimental and may be changed or removed in any future release. """ struct CompressibleNavierStokesDiffusion2D{GradientVariables, RealT <: Real, E <: AbstractCompressibleEulerEquations{2}} <: AbstractCompressibleNavierStokesDiffusion{2, 4} # TODO: parabolic @@ -90,15 +74,12 @@ struct CompressibleNavierStokesDiffusion2D{GradientVariables, RealT <: Real, E < # 2) Add NGRADS as a type parameter here and in AbstractEquationsParabolic, add `ngradients(...)` accessor function gamma::RealT # ratio of specific heats inv_gamma_minus_one::RealT # = inv(gamma - 1); can be used to write slow divisions as fast multiplications - Re::RealT # Reynolds number + + mu::RealT # viscosity Pr::RealT # Prandtl number - Ma_inf::RealT # free-stream Mach number + R::RealT # gas constant kappa::RealT # thermal diffusivity for Fick's law - p_inf::RealT # free-stream pressure - u_inf::RealT # free-stream velocity - R::RealT # gas constant (depends on nondimensional scaling!) - equations_hyperbolic::E # CompressibleEulerEquations2D gradient_variables::GradientVariables # GradientVariablesPrimitive or GradientVariablesEntropy end @@ -123,26 +104,19 @@ struct GradientVariablesEntropy end # default to primitive gradient variables function CompressibleNavierStokesDiffusion2D(equations::CompressibleEulerEquations2D; - Reynolds, Prandtl, Mach_freestream, + viscosity, Prandtl, gas_constant = 287.058, gradient_variables = GradientVariablesPrimitive()) gamma = equations.gamma inv_gamma_minus_one = equations.inv_gamma_minus_one - Re, Pr, Ma = promote(Reynolds, Prandtl, Mach_freestream) + mu, Pr, R = promote(viscosity, Prandtl, gas_constant) # Under the assumption of constant Prandtl number the thermal conductivity # constant is kappa = gamma μ R / ((gamma-1) Pr). - # Important note! Due to nondimensional scaling R = 1 / gamma, this constant - # simplifies slightly. Also, the factor of μ is accounted for later. - kappa = inv_gamma_minus_one / Pr - - # From the nondimensionalization discussed above set the remaining free-stream - # quantities - p_inf = 1 / gamma - u_inf = Mach_freestream - R = 1 / gamma + # Important note! Factor of μ is accounted for later in `flux`. + kappa = gamma * R * inv_gamma_minus_one / Pr + CompressibleNavierStokesDiffusion2D{typeof(gradient_variables), typeof(gamma), typeof(equations)}(gamma, inv_gamma_minus_one, - Re, Pr, Ma, kappa, - p_inf, u_inf, R, + mu, Pr, R, kappa, equations, gradient_variables) end @@ -159,15 +133,12 @@ varnames(variable_mapping, equations_parabolic::CompressibleNavierStokesDiffusio gradient_variable_transformation(::CompressibleNavierStokesDiffusion2D{GradientVariablesPrimitive}) = cons2prim gradient_variable_transformation(::CompressibleNavierStokesDiffusion2D{GradientVariablesEntropy}) = cons2entropy -# Explicit formulas for the diffussive Navier-Stokes fluxes are available, e.g. in Section 2 -# of the paper by Svärd, Carpenter and Nordström -# "A stable high-order finite difference scheme for the compressible Navier–Stokes -# equations, far-field boundary conditions" -# Although these authors use a different nondimensionalization so some constants are different -# particularly for Fick's law. -# -# Note, could be generalized to use Sutherland's law to get the molecular and thermal -# diffusivity + +# Explicit formulas for the diffusive Navier-Stokes fluxes are available, e.g., in Section 2 +# of the paper by Rueda-Ramíreza, Hennemann, Hindenlang, Winters, and Gassner +# "An Entropy Stable Nodal Discontinuous Galerkin Method for the resistive +# MHD Equations. Part II: Subcell Finite Volume Shock Capturing" +# where one sets the magnetic field components equal to 0. function flux(u, gradients, orientation::Integer, equations::CompressibleNavierStokesDiffusion2D) # Here, `u` is assumed to be the "transformed" variables specified by `gradient_variable_transformation`. rho, v1, v2, _ = convert_transformed_to_primitive(u, equations) @@ -188,13 +159,13 @@ function flux(u, gradients, orientation::Integer, equations::CompressibleNavierS tau_22 = 4.0 / 3.0 * dv2dy - 2.0 / 3.0 * dv1dx # Fick's law q = -kappa * grad(T); constant is kappa = gamma μ R / ((gamma-1) Pr) - # Important note! Due to nondimensional scaling R = 1 / gamma, so the - # temperature T in the gradient computation already contains a factor of gamma q1 = equations.kappa * dTdx q2 = equations.kappa * dTdy - # kinematic viscosity is simply 1/Re for this nondimensionalization - mu = 1.0 / equations.Re + # Constant viscosity is copied to a variable for readibility. + # Offers flexibility for viscosity via Sutherland's law where it depends + # on temperature and reference values, Ts and Tref such that mu(T) + mu = equations.mu if orientation == 1 # viscous flux components in the x-direction From 0b93d3e8bae67ab8242e1cef1bf059a391dd761f Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Thu, 6 Oct 2022 19:11:20 +0200 Subject: [PATCH 02/12] update TreeMesh elixirs and test values --- .../elixir_navier_stokes_convergence.jl | 61 +++++++++---------- .../elixir_navier_stokes_lid_driven_cavity.jl | 10 ++- test/test_parabolic_2d.jl | 28 ++++----- 3 files changed, 48 insertions(+), 51 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_navier_stokes_convergence.jl b/examples/tree_2d_dgsem/elixir_navier_stokes_convergence.jl index bd74271dcac..9855bf36488 100644 --- a/examples/tree_2d_dgsem/elixir_navier_stokes_convergence.jl +++ b/examples/tree_2d_dgsem/elixir_navier_stokes_convergence.jl @@ -4,12 +4,12 @@ using Trixi ############################################################################### # semidiscretization of the ideal compressible Navier-Stokes equations -reynolds_number() = 100 prandtl_number() = 0.72 +viscosity() = 0.01 equations = CompressibleEulerEquations2D(1.4) -equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, Reynolds=reynolds_number(), Prandtl=prandtl_number(), - Mach_freestream=0.5, gradient_variables=GradientVariablesPrimitive()) +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, viscosity=viscosity(), Prandtl=prandtl_number(), + gradient_variables=GradientVariablesPrimitive()) # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs, @@ -54,7 +54,7 @@ end # see also https://github.com/trixi-framework/Trixi.jl/pull/1160 inv_gamma_minus_one = inv(equations.gamma - 1) Pr = prandtl_number() - Re = reynolds_number() + mu = viscosity() # Same settings as in `initial_condition` # Amplitude and shift @@ -106,7 +106,6 @@ end # Some convenience constants T_const = equations.gamma * inv_gamma_minus_one / Pr - inv_Re = 1.0 / Re inv_rho_cubed = 1.0 / (rho^3) # compute the source terms @@ -120,10 +119,10 @@ end + rho * v1_y * v2 + rho * v1 * v2_y # stress tensor from x-direction - - 4.0 / 3.0 * v1_xx * inv_Re - + 2.0 / 3.0 * v2_xy * inv_Re - - v1_yy * inv_Re - - v2_xy * inv_Re ) + - 4.0 / 3.0 * v1_xx * mu + + 2.0 / 3.0 * v2_xy * mu + - v1_yy * mu + - v2_xy * mu ) # y-momentum equation du3 = ( rho_t * v2 + rho * v2_t + p_y + rho_x * v1 * v2 + rho * v1_x * v2 @@ -131,39 +130,39 @@ end + rho_y * v2^2 + 2.0 * rho * v2 * v2_y # stress tensor from y-direction - - v1_xy * inv_Re - - v2_xx * inv_Re - - 4.0 / 3.0 * v2_yy * inv_Re - + 2.0 / 3.0 * v1_xy * inv_Re ) + - v1_xy * mu + - v2_xx * mu + - 4.0 / 3.0 * v2_yy * mu + + 2.0 / 3.0 * v1_xy * mu ) # total energy equation du4 = ( E_t + v1_x * (E + p) + v1 * (E_x + p_x) + v2_y * (E + p) + v2 * (E_y + p_y) # stress tensor and temperature gradient terms from x-direction - - 4.0 / 3.0 * v1_xx * v1 * inv_Re - + 2.0 / 3.0 * v2_xy * v1 * inv_Re - - 4.0 / 3.0 * v1_x * v1_x * inv_Re - + 2.0 / 3.0 * v2_y * v1_x * inv_Re - - v1_xy * v2 * inv_Re - - v2_xx * v2 * inv_Re - - v1_y * v2_x * inv_Re - - v2_x * v2_x * inv_Re + - 4.0 / 3.0 * v1_xx * v1 * mu + + 2.0 / 3.0 * v2_xy * v1 * mu + - 4.0 / 3.0 * v1_x * v1_x * mu + + 2.0 / 3.0 * v2_y * v1_x * mu + - v1_xy * v2 * mu + - v2_xx * v2 * mu + - v1_y * v2_x * mu + - v2_x * v2_x * mu - T_const * inv_rho_cubed * ( p_xx * rho * rho - 2.0 * p_x * rho * rho_x + 2.0 * p * rho_x * rho_x - - p * rho * rho_xx ) * inv_Re + - p * rho * rho_xx ) * mu # stress tensor and temperature gradient terms from y-direction - - v1_yy * v1 * inv_Re - - v2_xy * v1 * inv_Re - - v1_y * v1_y * inv_Re - - v2_x * v1_y * inv_Re - - 4.0 / 3.0 * v2_yy * v2 * inv_Re - + 2.0 / 3.0 * v1_xy * v2 * inv_Re - - 4.0 / 3.0 * v2_y * v2_y * inv_Re - + 2.0 / 3.0 * v1_x * v2_y * inv_Re + - v1_yy * v1 * mu + - v2_xy * v1 * mu + - v1_y * v1_y * mu + - v2_x * v1_y * mu + - 4.0 / 3.0 * v2_yy * v2 * mu + + 2.0 / 3.0 * v1_xy * v2 * mu + - 4.0 / 3.0 * v2_y * v2_y * mu + + 2.0 / 3.0 * v1_x * v2_y * mu - T_const * inv_rho_cubed * ( p_yy * rho * rho - 2.0 * p_y * rho * rho_y + 2.0 * p * rho_y * rho_y - - p * rho * rho_yy ) * inv_Re ) + - p * rho * rho_yy ) * mu ) return SVector(du1, du2, du3, du4) end diff --git a/examples/tree_2d_dgsem/elixir_navier_stokes_lid_driven_cavity.jl b/examples/tree_2d_dgsem/elixir_navier_stokes_lid_driven_cavity.jl index 89882675613..afc907fbe78 100644 --- a/examples/tree_2d_dgsem/elixir_navier_stokes_lid_driven_cavity.jl +++ b/examples/tree_2d_dgsem/elixir_navier_stokes_lid_driven_cavity.jl @@ -5,14 +5,12 @@ using Trixi # semidiscretization of the ideal compressible Navier-Stokes equations # TODO: parabolic; unify names of these accessor functions -reynolds_number() = 1000.0 prandtl_number() = 0.72 -mach_number() = 0.1 +viscosity() = 0.001 equations = CompressibleEulerEquations2D(1.4) -equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, Reynolds=reynolds_number(), - Prandtl=prandtl_number(), - Mach_freestream=mach_number()) +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, viscosity=viscosity(), + Prandtl=prandtl_number()) # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs) @@ -28,7 +26,7 @@ mesh = TreeMesh(coordinates_min, coordinates_max, function initial_condition_cavity(x, t, equations::CompressibleEulerEquations2D) - Ma = mach_number() + Ma = 0.1 rho = 1.0 u, v = 0, 0 p = 1.0 / (Ma^2 * equations.gamma) diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index e0552217b3e..9ce66147f9c 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -129,8 +129,8 @@ isdir(outdir) && rm(outdir, recursive=true) @trixi_testset "TreeMesh2D: elixir_navier_stokes_convergence.jl" begin @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_navier_stokes_convergence.jl"), initial_refinement_level = 2, tspan=(0.0, 0.1), - l2 = [0.0021116725306635146, 0.0034322351490824465, 0.003874252819611102, 0.012469246082522416], - linf = [0.012006418939297214, 0.03552087120958058, 0.02451274749176294, 0.11191122588577151] + l2 = [0.002111672530658797, 0.0034322351490857846, 0.0038742528195910416, 0.012469246082568561], + linf = [0.012006418939223495, 0.035520871209746126, 0.024512747492231427, 0.11191122588756564] ) end @@ -138,25 +138,25 @@ isdir(outdir) && rm(outdir, recursive=true) @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_navier_stokes_convergence.jl"), initial_refinement_level = 2, tspan=(0.0, 0.1), heat_bc_top_bottom=Isothermal((x, t, equations) -> Trixi.temperature(initial_condition_navier_stokes_convergence_test(x, t, equations), equations)), - l2 = [0.002103629650384378, 0.0034358439333976123, 0.0038673598780978413, 0.012670355349347209], - linf = [0.012006261793021222, 0.035502125190110666, 0.025107947320650532, 0.11647078036915026] + l2 = [0.002103629650383915, 0.003435843933396454, 0.00386735987813341, 0.012670355349235728], + linf = [0.012006261793147788, 0.03550212518982032, 0.025107947319661185, 0.11647078036571124] ) end @trixi_testset "TreeMesh2D: elixir_navier_stokes_convergence.jl (Entropy gradient variables)" begin @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_navier_stokes_convergence.jl"), - initial_refinement_level=2, tspan=(0,.1), gradient_variables=GradientVariablesEntropy(), - l2 = [0.002140374251726679, 0.0034258287094981717, 0.0038915122887464865, 0.012506862342821999], - linf = [0.012244412004805971, 0.035075591861236655, 0.02458089234452718, 0.11425600757951138] + initial_refinement_level=2, tspan=(0.0, 0.1), gradient_variables=GradientVariablesEntropy(), + l2 = [0.0021403742517389513, 0.0034258287094908572, 0.0038915122886898517, 0.012506862343013842], + linf = [0.012244412004628336, 0.03507559186162224, 0.024580892345558894, 0.11425600758350107] ) end @trixi_testset "TreeMesh2D: elixir_navier_stokes_convergence.jl (Entropy gradient variables, isothermal walls)" begin @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_navier_stokes_convergence.jl"), - initial_refinement_level=2, tspan=(0,.1), gradient_variables=GradientVariablesEntropy(), + initial_refinement_level=2, tspan=(0.0, 0.1), gradient_variables=GradientVariablesEntropy(), heat_bc_top_bottom=Isothermal((x, t, equations) -> Trixi.temperature(initial_condition_navier_stokes_convergence_test(x, t, equations), equations)), - l2 = [0.002134973734788134, 0.0034301388278191753, 0.0038928324474145994, 0.012693611436279086], - linf = [0.012244236275815057, 0.035054066314196344, 0.02509959850525358, 0.1179561632485715] + l2 = [0.0021349737347844907, 0.0034301388278203033, 0.0038928324474291572, 0.012693611436230873], + linf = [0.01224423627586213, 0.035054066314102905, 0.025099598504931965, 0.11795616324751634] ) end @@ -164,16 +164,16 @@ isdir(outdir) && rm(outdir, recursive=true) @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_navier_stokes_convergence.jl"), initial_refinement_level = 2, tspan=(0.0, 0.1), volume_integral=VolumeIntegralFluxDifferencing(flux_central), - l2 = [0.0021116725306635146, 0.0034322351490824465, 0.003874252819611102, 0.012469246082522416], - linf = [0.012006418939297214, 0.03552087120958058, 0.02451274749176294, 0.11191122588577151] + l2 = [0.0021116725306633594, 0.0034322351490827557, 0.0038742528196093542, 0.012469246082526909], + linf = [[0.012006418939291663, 0.035520871209594115, 0.024512747491801577, 0.11191122588591007] ) end @trixi_testset "TreeMesh2D: elixir_navier_stokes_lid_driven_cavity.jl" begin @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_navier_stokes_lid_driven_cavity.jl"), initial_refinement_level = 2, tspan=(0.0, 0.5), - l2 = [0.0001514457152968994, 0.018766076072331786, 0.007065070765651992, 0.020839900573430787], - linf = [0.0014523369373645734, 0.12366779944955876, 0.055324509971157544, 0.1609992780534526] + l2 = [0.00015144571529699053, 0.018766076072331623, 0.007065070765652574, 0.0208399005734258], + linf = [0.0014523369373669048, 0.12366779944955864, 0.05532450997115432, 0.16099927805328207] ) end From b23b9a739fa0677aa2ea1235be6f92f58bfe537c Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Thu, 6 Oct 2022 19:18:53 +0200 Subject: [PATCH 03/12] update the elixirs and test values for DGMulti --- .../elixir_navier_stokes_convergence.jl | 61 +++++++++---------- .../elixir_navier_stokes_lid_driven_cavity.jl | 11 ++-- test/test_parabolic_2d.jl | 8 +-- 3 files changed, 39 insertions(+), 41 deletions(-) diff --git a/examples/dgmulti_2d/elixir_navier_stokes_convergence.jl b/examples/dgmulti_2d/elixir_navier_stokes_convergence.jl index 6d64f662e77..8c19082036d 100644 --- a/examples/dgmulti_2d/elixir_navier_stokes_convergence.jl +++ b/examples/dgmulti_2d/elixir_navier_stokes_convergence.jl @@ -4,14 +4,14 @@ using Trixi ############################################################################### # semidiscretization of the ideal compressible Navier-Stokes equations -reynolds_number() = 100 prandtl_number() = 0.72 +viscosity() = 0.01 equations = CompressibleEulerEquations2D(1.4) # Note: If you change the Navier-Stokes parameters here, also change them in the initial condition # I really do not like this structure but it should work for now -equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, Reynolds=reynolds_number(), Prandtl=prandtl_number(), - Mach_freestream=0.5, gradient_variables=GradientVariablesPrimitive()) +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, viscosity=viscosity(), Prandtl=prandtl_number(), + gradient_variables=GradientVariablesPrimitive()) # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux dg = DGMulti(polydeg = 3, element_type = Tri(), approximation_type = Polynomial(), @@ -52,7 +52,7 @@ end # see also https://github.com/trixi-framework/Trixi.jl/pull/1160 inv_gamma_minus_one = inv(equations.gamma - 1) Pr = prandtl_number() - Re = reynolds_number() + mu = viscosity() # Same settings as in `initial_condition` # Amplitude and shift @@ -104,7 +104,6 @@ end # Some convenience constants T_const = equations.gamma * inv_gamma_minus_one / Pr - inv_Re = 1.0 / Re inv_rho_cubed = 1.0 / (rho^3) # compute the source terms @@ -118,10 +117,10 @@ end + rho * v1_y * v2 + rho * v1 * v2_y # stress tensor from x-direction - - 4.0 / 3.0 * v1_xx * inv_Re - + 2.0 / 3.0 * v2_xy * inv_Re - - v1_yy * inv_Re - - v2_xy * inv_Re ) + - 4.0 / 3.0 * v1_xx * mu + + 2.0 / 3.0 * v2_xy * mu + - v1_yy * mu + - v2_xy * mu ) # y-momentum equation du3 = ( rho_t * v2 + rho * v2_t + p_y + rho_x * v1 * v2 + rho * v1_x * v2 @@ -129,39 +128,39 @@ end + rho_y * v2^2 + 2.0 * rho * v2 * v2_y # stress tensor from y-direction - - v1_xy * inv_Re - - v2_xx * inv_Re - - 4.0 / 3.0 * v2_yy * inv_Re - + 2.0 / 3.0 * v1_xy * inv_Re ) + - v1_xy * mu + - v2_xx * mu + - 4.0 / 3.0 * v2_yy * mu + + 2.0 / 3.0 * v1_xy * mu ) # total energy equation du4 = ( E_t + v1_x * (E + p) + v1 * (E_x + p_x) + v2_y * (E + p) + v2 * (E_y + p_y) # stress tensor and temperature gradient terms from x-direction - - 4.0 / 3.0 * v1_xx * v1 * inv_Re - + 2.0 / 3.0 * v2_xy * v1 * inv_Re - - 4.0 / 3.0 * v1_x * v1_x * inv_Re - + 2.0 / 3.0 * v2_y * v1_x * inv_Re - - v1_xy * v2 * inv_Re - - v2_xx * v2 * inv_Re - - v1_y * v2_x * inv_Re - - v2_x * v2_x * inv_Re + - 4.0 / 3.0 * v1_xx * v1 * mu + + 2.0 / 3.0 * v2_xy * v1 * mu + - 4.0 / 3.0 * v1_x * v1_x * mu + + 2.0 / 3.0 * v2_y * v1_x * mu + - v1_xy * v2 * mu + - v2_xx * v2 * mu + - v1_y * v2_x * mu + - v2_x * v2_x * mu - T_const * inv_rho_cubed * ( p_xx * rho * rho - 2.0 * p_x * rho * rho_x + 2.0 * p * rho_x * rho_x - - p * rho * rho_xx ) * inv_Re + - p * rho * rho_xx ) * mu # stress tensor and temperature gradient terms from y-direction - - v1_yy * v1 * inv_Re - - v2_xy * v1 * inv_Re - - v1_y * v1_y * inv_Re - - v2_x * v1_y * inv_Re - - 4.0 / 3.0 * v2_yy * v2 * inv_Re - + 2.0 / 3.0 * v1_xy * v2 * inv_Re - - 4.0 / 3.0 * v2_y * v2_y * inv_Re - + 2.0 / 3.0 * v1_x * v2_y * inv_Re + - v1_yy * v1 * mu + - v2_xy * v1 * mu + - v1_y * v1_y * mu + - v2_x * v1_y * mu + - 4.0 / 3.0 * v2_yy * v2 * mu + + 2.0 / 3.0 * v1_xy * v2 * mu + - 4.0 / 3.0 * v2_y * v2_y * mu + + 2.0 / 3.0 * v1_x * v2_y * mu - T_const * inv_rho_cubed * ( p_yy * rho * rho - 2.0 * p_y * rho * rho_y + 2.0 * p * rho_y * rho_y - - p * rho * rho_yy ) * inv_Re ) + - p * rho * rho_yy ) * mu ) return SVector(du1, du2, du3, du4) end diff --git a/examples/dgmulti_2d/elixir_navier_stokes_lid_driven_cavity.jl b/examples/dgmulti_2d/elixir_navier_stokes_lid_driven_cavity.jl index be8ad1a0730..ebe04b0e799 100644 --- a/examples/dgmulti_2d/elixir_navier_stokes_lid_driven_cavity.jl +++ b/examples/dgmulti_2d/elixir_navier_stokes_lid_driven_cavity.jl @@ -5,14 +5,13 @@ using Trixi # semidiscretization of the ideal compressible Navier-Stokes equations # TODO: parabolic; unify names of these accessor functions -reynolds_number() = 1000.0 prandtl_number() = 0.72 -mach_number() = 0.1 +viscosity() = 0.001 equations = CompressibleEulerEquations2D(1.4) -equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, Reynolds=reynolds_number(), - Prandtl=prandtl_number(), - Mach_freestream=mach_number()) +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, viscosity=viscosity(), + Prandtl=prandtl_number()) + # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux dg = DGMulti(polydeg = 3, element_type = Quad(), approximation_type = GaussSBP(), @@ -25,7 +24,7 @@ is_on_boundary = Dict(:top => top, :rest_of_boundary => rest_of_boundary) mesh = DGMultiMesh(dg, cells_per_dimension=(16, 16); is_on_boundary) function initial_condition_cavity(x, t, equations::CompressibleEulerEquations2D) - Ma = mach_number() + Ma = 0.1 rho = 1.0 u, v = 0, 0 p = 1.0 / (Ma^2 * equations.gamma) diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 9ce66147f9c..24e1bbb5458 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -97,16 +97,16 @@ isdir(outdir) && rm(outdir, recursive=true) @trixi_testset "DGMulti: elixir_navier_stokes_convergence.jl" begin @test_trixi_include(joinpath(examples_dir(), "dgmulti_2d", "elixir_navier_stokes_convergence.jl"), cells_per_dimension = (4, 4), tspan=(0.0, 0.1), - l2 = [0.00153550768125133, 0.0033843168272696357, 0.0036531858107444067, 0.009948436427519428], - linf = [0.005522560467190019, 0.013425258500731063, 0.013962115643483375, 0.027483102120516634] + l2 = [0.0015355076812510957, 0.0033843168272696756, 0.0036531858107443434, 0.009948436427519214], + linf = [0.005522560467190019, 0.013425258500730508, 0.013962115643482154, 0.027483102120502423] ) end @trixi_testset "DGMulti: elixir_navier_stokes_lid_driven_cavity.jl" begin @test_trixi_include(joinpath(examples_dir(), "dgmulti_2d", "elixir_navier_stokes_lid_driven_cavity.jl"), cells_per_dimension = (4, 4), tspan=(0.0, 0.5), - l2 = [0.0002215612522711349, 0.028318325921400257, 0.009509168701069035, 0.028267900513539248], - linf = [0.0015622789413053395, 0.14886653390741342, 0.07163235655334241, 0.19472785105216417] + l2 = [0.00022156125227115747, 0.028318325921401, 0.009509168701070296, 0.028267900513550506], + linf = [0.001562278941298234, 0.14886653390744856, 0.0716323565533752, 0.19472785105241996] ) end From df33b79a198dabcee5f5a859f12e72f7eca784ad Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Thu, 6 Oct 2022 19:20:08 +0200 Subject: [PATCH 04/12] fix unbalanced parentheses --- test/test_parabolic_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 24e1bbb5458..4b31f6c1de8 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -165,7 +165,7 @@ isdir(outdir) && rm(outdir, recursive=true) initial_refinement_level = 2, tspan=(0.0, 0.1), volume_integral=VolumeIntegralFluxDifferencing(flux_central), l2 = [0.0021116725306633594, 0.0034322351490827557, 0.0038742528196093542, 0.012469246082526909], - linf = [[0.012006418939291663, 0.035520871209594115, 0.024512747491801577, 0.11191122588591007] + linf = [0.012006418939291663, 0.035520871209594115, 0.024512747491801577, 0.11191122588591007] ) end From 356630ecb64669db394dcbca8e368585070abd4c Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Thu, 6 Oct 2022 19:25:56 +0200 Subject: [PATCH 05/12] file renaming for consistency --- ...r_stokes_convergence.jl => elixir_navierstokes_convergence.jl} | 0 ..._driven_cavity.jl => elixir_navierstokes_lid_driven_cavity.jl} | 0 ...r_stokes_convergence.jl => elixir_navierstokes_convergence.jl} | 0 ..._driven_cavity.jl => elixir_navierstokes_lid_driven_cavity.jl} | 0 4 files changed, 0 insertions(+), 0 deletions(-) rename examples/dgmulti_2d/{elixir_navier_stokes_convergence.jl => elixir_navierstokes_convergence.jl} (100%) rename examples/dgmulti_2d/{elixir_navier_stokes_lid_driven_cavity.jl => elixir_navierstokes_lid_driven_cavity.jl} (100%) rename examples/tree_2d_dgsem/{elixir_navier_stokes_convergence.jl => elixir_navierstokes_convergence.jl} (100%) rename examples/tree_2d_dgsem/{elixir_navier_stokes_lid_driven_cavity.jl => elixir_navierstokes_lid_driven_cavity.jl} (100%) diff --git a/examples/dgmulti_2d/elixir_navier_stokes_convergence.jl b/examples/dgmulti_2d/elixir_navierstokes_convergence.jl similarity index 100% rename from examples/dgmulti_2d/elixir_navier_stokes_convergence.jl rename to examples/dgmulti_2d/elixir_navierstokes_convergence.jl diff --git a/examples/dgmulti_2d/elixir_navier_stokes_lid_driven_cavity.jl b/examples/dgmulti_2d/elixir_navierstokes_lid_driven_cavity.jl similarity index 100% rename from examples/dgmulti_2d/elixir_navier_stokes_lid_driven_cavity.jl rename to examples/dgmulti_2d/elixir_navierstokes_lid_driven_cavity.jl diff --git a/examples/tree_2d_dgsem/elixir_navier_stokes_convergence.jl b/examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl similarity index 100% rename from examples/tree_2d_dgsem/elixir_navier_stokes_convergence.jl rename to examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl diff --git a/examples/tree_2d_dgsem/elixir_navier_stokes_lid_driven_cavity.jl b/examples/tree_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl similarity index 100% rename from examples/tree_2d_dgsem/elixir_navier_stokes_lid_driven_cavity.jl rename to examples/tree_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl From c3620c6066fffa5c1ff46f964c3d54f7c26c7a00 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Thu, 6 Oct 2022 19:26:21 +0200 Subject: [PATCH 06/12] update file names in the parabolic test script --- test/test_parabolic_2d.jl | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 4b31f6c1de8..d9c3967315a 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -94,16 +94,16 @@ isdir(outdir) && rm(outdir, recursive=true) ) end - @trixi_testset "DGMulti: elixir_navier_stokes_convergence.jl" begin - @test_trixi_include(joinpath(examples_dir(), "dgmulti_2d", "elixir_navier_stokes_convergence.jl"), + @trixi_testset "DGMulti: elixir_navierstokes_convergence.jl" begin + @test_trixi_include(joinpath(examples_dir(), "dgmulti_2d", "elixir_navierstokes_convergence.jl"), cells_per_dimension = (4, 4), tspan=(0.0, 0.1), l2 = [0.0015355076812510957, 0.0033843168272696756, 0.0036531858107443434, 0.009948436427519214], linf = [0.005522560467190019, 0.013425258500730508, 0.013962115643482154, 0.027483102120502423] ) end - @trixi_testset "DGMulti: elixir_navier_stokes_lid_driven_cavity.jl" begin - @test_trixi_include(joinpath(examples_dir(), "dgmulti_2d", "elixir_navier_stokes_lid_driven_cavity.jl"), + @trixi_testset "DGMulti: elixir_navierstokes_lid_driven_cavity.jl" begin + @test_trixi_include(joinpath(examples_dir(), "dgmulti_2d", "elixir_navierstokes_lid_driven_cavity.jl"), cells_per_dimension = (4, 4), tspan=(0.0, 0.5), l2 = [0.00022156125227115747, 0.028318325921401, 0.009509168701070296, 0.028267900513550506], linf = [0.001562278941298234, 0.14886653390744856, 0.0716323565533752, 0.19472785105241996] @@ -126,16 +126,16 @@ isdir(outdir) && rm(outdir, recursive=true) ) end - @trixi_testset "TreeMesh2D: elixir_navier_stokes_convergence.jl" begin - @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_navier_stokes_convergence.jl"), + @trixi_testset "TreeMesh2D: elixir_navierstokes_convergence.jl" begin + @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_navierstokes_convergence.jl"), initial_refinement_level = 2, tspan=(0.0, 0.1), l2 = [0.002111672530658797, 0.0034322351490857846, 0.0038742528195910416, 0.012469246082568561], linf = [0.012006418939223495, 0.035520871209746126, 0.024512747492231427, 0.11191122588756564] ) end - @trixi_testset "TreeMesh2D: elixir_navier_stokes_convergence.jl (isothermal walls)" begin - @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_navier_stokes_convergence.jl"), + @trixi_testset "TreeMesh2D: elixir_navierstokes_convergence.jl (isothermal walls)" begin + @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_navierstokes_convergence.jl"), initial_refinement_level = 2, tspan=(0.0, 0.1), heat_bc_top_bottom=Isothermal((x, t, equations) -> Trixi.temperature(initial_condition_navier_stokes_convergence_test(x, t, equations), equations)), l2 = [0.002103629650383915, 0.003435843933396454, 0.00386735987813341, 0.012670355349235728], @@ -143,16 +143,16 @@ isdir(outdir) && rm(outdir, recursive=true) ) end - @trixi_testset "TreeMesh2D: elixir_navier_stokes_convergence.jl (Entropy gradient variables)" begin - @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_navier_stokes_convergence.jl"), + @trixi_testset "TreeMesh2D: elixir_navierstokes_convergence.jl (Entropy gradient variables)" begin + @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_navierstokes_convergence.jl"), initial_refinement_level=2, tspan=(0.0, 0.1), gradient_variables=GradientVariablesEntropy(), l2 = [0.0021403742517389513, 0.0034258287094908572, 0.0038915122886898517, 0.012506862343013842], linf = [0.012244412004628336, 0.03507559186162224, 0.024580892345558894, 0.11425600758350107] ) end - @trixi_testset "TreeMesh2D: elixir_navier_stokes_convergence.jl (Entropy gradient variables, isothermal walls)" begin - @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_navier_stokes_convergence.jl"), + @trixi_testset "TreeMesh2D: elixir_navierstokes_convergence.jl (Entropy gradient variables, isothermal walls)" begin + @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_navierstokes_convergence.jl"), initial_refinement_level=2, tspan=(0.0, 0.1), gradient_variables=GradientVariablesEntropy(), heat_bc_top_bottom=Isothermal((x, t, equations) -> Trixi.temperature(initial_condition_navier_stokes_convergence_test(x, t, equations), equations)), l2 = [0.0021349737347844907, 0.0034301388278203033, 0.0038928324474291572, 0.012693611436230873], @@ -160,8 +160,8 @@ isdir(outdir) && rm(outdir, recursive=true) ) end - @trixi_testset "TreeMesh2D: elixir_navier_stokes_convergence.jl (flux differencing)" begin - @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_navier_stokes_convergence.jl"), + @trixi_testset "TreeMesh2D: elixir_navierstokes_convergence.jl (flux differencing)" begin + @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_navierstokes_convergence.jl"), initial_refinement_level = 2, tspan=(0.0, 0.1), volume_integral=VolumeIntegralFluxDifferencing(flux_central), l2 = [0.0021116725306633594, 0.0034322351490827557, 0.0038742528196093542, 0.012469246082526909], @@ -169,8 +169,8 @@ isdir(outdir) && rm(outdir, recursive=true) ) end - @trixi_testset "TreeMesh2D: elixir_navier_stokes_lid_driven_cavity.jl" begin - @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_navier_stokes_lid_driven_cavity.jl"), + @trixi_testset "TreeMesh2D: elixir_navierstokes_lid_driven_cavity.jl" begin + @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_navierstokes_lid_driven_cavity.jl"), initial_refinement_level = 2, tspan=(0.0, 0.5), l2 = [0.00015144571529699053, 0.018766076072331623, 0.007065070765652574, 0.0208399005734258], linf = [0.0014523369373669048, 0.12366779944955864, 0.05532450997115432, 0.16099927805328207] From 605b3984f2fd1152098429e613685e3a3dd3b3a1 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 7 Oct 2022 09:35:37 +0200 Subject: [PATCH 07/12] remove the gas constant as a parameter and update comments --- .../compressible_navier_stokes_2d.jl | 43 +++++++++++-------- 1 file changed, 25 insertions(+), 18 deletions(-) diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index 4c78add1e19..f53bd9fd48f 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -1,5 +1,5 @@ @doc raw""" - CompressibleNavierStokesDiffusion2D(gamma, mu, Pr, R, equations, + CompressibleNavierStokesDiffusion2D(gamma, mu, Pr, equations, gradient_variables=GradientVariablesPrimitive()) These equations contain the diffusion (i.e. parabolic) terms applied @@ -7,9 +7,8 @@ to mass, momenta, and total energy together with the advective terms from the [`CompressibleEulerEquations2D`](@ref). - `gamma`: adiabatic constant, -- `mu`: viscosity, +- `mu`: dynamic viscosity, - `Pr`: Prandtl number, -- `R`: gas constant - `equations`: instance of the [`CompressibleEulerEquations2D`](@ref) - `gradient_variables`: which variables the gradients are taken with respect to. Defaults to `GradientVariablesPrimitive()`. @@ -50,6 +49,12 @@ the thermal conductivity is ```math \kappa = \frac{\gamma \mu R}{(\gamma - 1)\textrm{Pr}}. ``` +From this combination of temperature ``T`` and thermal conductivity ``\kappa`` we see +that the gas constant `R` cancels and the heat flux becomes +```math +\nabla q = -\kappa\nabla\left(T\right) = -\frac{\gamma \mu}{(\gamma - 1)\textrm{Pr}}\nabla\left(\frac{p}{\rho}\right) +``` +which is the form implemented below in the [`flux`](@ref) function. In two spatial dimensions we require gradients for three quantities, e.g., primitive quantities @@ -77,7 +82,6 @@ struct CompressibleNavierStokesDiffusion2D{GradientVariables, RealT <: Real, E < mu::RealT # viscosity Pr::RealT # Prandtl number - R::RealT # gas constant kappa::RealT # thermal diffusivity for Fick's law equations_hyperbolic::E # CompressibleEulerEquations2D @@ -104,19 +108,19 @@ struct GradientVariablesEntropy end # default to primitive gradient variables function CompressibleNavierStokesDiffusion2D(equations::CompressibleEulerEquations2D; - viscosity, Prandtl, gas_constant = 287.058, + viscosity, Prandtl, gradient_variables = GradientVariablesPrimitive()) gamma = equations.gamma inv_gamma_minus_one = equations.inv_gamma_minus_one - mu, Pr, R = promote(viscosity, Prandtl, gas_constant) + mu, Pr = promote(viscosity, Prandtl) # Under the assumption of constant Prandtl number the thermal conductivity - # constant is kappa = gamma μ R / ((gamma-1) Pr). + # constant is kappa = gamma μ / ((gamma-1) Pr). # Important note! Factor of μ is accounted for later in `flux`. - kappa = gamma * R * inv_gamma_minus_one / Pr + kappa = gamma * inv_gamma_minus_one / Pr CompressibleNavierStokesDiffusion2D{typeof(gradient_variables), typeof(gamma), typeof(equations)}(gamma, inv_gamma_minus_one, - mu, Pr, R, kappa, + mu, Pr, kappa, equations, gradient_variables) end @@ -158,12 +162,15 @@ function flux(u, gradients, orientation::Integer, equations::CompressibleNavierS # (4/3 * (v2)_y - 2/3 * (v1)_x) tau_22 = 4.0 / 3.0 * dv2dy - 2.0 / 3.0 * dv1dx - # Fick's law q = -kappa * grad(T); constant is kappa = gamma μ R / ((gamma-1) Pr) + # Fick's law q = -kappa * grad(T) = -kappa * grad(p / (R rho)) + # with thermal diffusivity constant kappa = gamma μ R / ((gamma-1) Pr) + # Note, the gas constant cancels under this formulation, so it is not present + # in the implementation q1 = equations.kappa * dTdx q2 = equations.kappa * dTdy - # Constant viscosity is copied to a variable for readibility. - # Offers flexibility for viscosity via Sutherland's law where it depends + # Constant dynamic viscosity is copied to a variable for readibility. + # 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 @@ -243,9 +250,9 @@ end T = temperature(u, equations) return SVector(gradient_entropy_vars[1], - equations.R * T * (gradient_entropy_vars[2] + v1 * gradient_entropy_vars[4]), # grad(u) = R*T*(grad(w_2)+v1*grad(w_4)) - equations.R * T * (gradient_entropy_vars[3] + v2 * gradient_entropy_vars[4]), # grad(v) = R*T*(grad(w_3)+v2*grad(w_4)) - equations.R * T * T * gradient_entropy_vars[4] # grad(T) = R*T^2*grad(w_4)) + T * (gradient_entropy_vars[2] + v1 * gradient_entropy_vars[4]), # grad(u) = T*(grad(w_2)+v1*grad(w_4)) + T * (gradient_entropy_vars[3] + v2 * gradient_entropy_vars[4]), # grad(v) = T*(grad(w_3)+v2*grad(w_4)) + T * T * gradient_entropy_vars[4] # grad(T) = T^2*grad(w_4)) ) end @@ -262,7 +269,7 @@ end rho, rho_v1, rho_v2, rho_e = u p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1^2 + rho_v2^2) / rho) - T = p / (equations.R * rho) + T = p / rho return T end @@ -385,8 +392,8 @@ end v1, v2 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, equations) T = boundary_condition.boundary_condition_heat_flux.boundary_value_function(x, t, equations) - # the entropy variables w2 = rho * v1 / p = v1 / (equations.R * T) = -v1 * w4. Similarly for w3 - w4 = -1 / (equations.R * T) + # the entropy variables w2 = rho * v1 / p = v1 / T = -v1 * w4. Similarly for w3 + w4 = -1 / T return SVector(w_inner[1], -v1 * w4, -v2 * w4, w4) end From 595ef59a04cd1be93af6f76d8e9a5dcf32132e57 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 7 Oct 2022 09:45:04 +0200 Subject: [PATCH 08/12] rename viscosity to dynamic_viscosity --- examples/dgmulti_2d/elixir_navierstokes_convergence.jl | 6 +++--- .../dgmulti_2d/elixir_navierstokes_lid_driven_cavity.jl | 4 ++-- examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl | 6 +++--- .../tree_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl | 4 ++-- 4 files changed, 10 insertions(+), 10 deletions(-) diff --git a/examples/dgmulti_2d/elixir_navierstokes_convergence.jl b/examples/dgmulti_2d/elixir_navierstokes_convergence.jl index 8c19082036d..ecb545cf9e6 100644 --- a/examples/dgmulti_2d/elixir_navierstokes_convergence.jl +++ b/examples/dgmulti_2d/elixir_navierstokes_convergence.jl @@ -5,12 +5,12 @@ using Trixi # semidiscretization of the ideal compressible Navier-Stokes equations prandtl_number() = 0.72 -viscosity() = 0.01 +dynamic_viscosity() = 0.01 equations = CompressibleEulerEquations2D(1.4) # Note: If you change the Navier-Stokes parameters here, also change them in the initial condition # I really do not like this structure but it should work for now -equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, viscosity=viscosity(), Prandtl=prandtl_number(), +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, viscosity=dynamic_viscosity(), Prandtl=prandtl_number(), gradient_variables=GradientVariablesPrimitive()) # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux @@ -52,7 +52,7 @@ end # see also https://github.com/trixi-framework/Trixi.jl/pull/1160 inv_gamma_minus_one = inv(equations.gamma - 1) Pr = prandtl_number() - mu = viscosity() + mu = dynamic_viscosity() # Same settings as in `initial_condition` # Amplitude and shift diff --git a/examples/dgmulti_2d/elixir_navierstokes_lid_driven_cavity.jl b/examples/dgmulti_2d/elixir_navierstokes_lid_driven_cavity.jl index ebe04b0e799..ed9008ec351 100644 --- a/examples/dgmulti_2d/elixir_navierstokes_lid_driven_cavity.jl +++ b/examples/dgmulti_2d/elixir_navierstokes_lid_driven_cavity.jl @@ -6,10 +6,10 @@ using Trixi # TODO: parabolic; unify names of these accessor functions prandtl_number() = 0.72 -viscosity() = 0.001 +dynamic_viscosity() = 0.001 equations = CompressibleEulerEquations2D(1.4) -equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, viscosity=viscosity(), +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, viscosity=dynamic_viscosity(), Prandtl=prandtl_number()) diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl b/examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl index 9855bf36488..1c8960111c7 100644 --- a/examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl +++ b/examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl @@ -5,10 +5,10 @@ using Trixi # semidiscretization of the ideal compressible Navier-Stokes equations prandtl_number() = 0.72 -viscosity() = 0.01 +dynamic_viscosity() = 0.01 equations = CompressibleEulerEquations2D(1.4) -equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, viscosity=viscosity(), Prandtl=prandtl_number(), +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, viscosity=dynamic_viscosity(), Prandtl=prandtl_number(), gradient_variables=GradientVariablesPrimitive()) # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux @@ -54,7 +54,7 @@ end # see also https://github.com/trixi-framework/Trixi.jl/pull/1160 inv_gamma_minus_one = inv(equations.gamma - 1) Pr = prandtl_number() - mu = viscosity() + mu = dynamic_viscosity() # Same settings as in `initial_condition` # Amplitude and shift 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 afc907fbe78..fda93ab87e2 100644 --- a/examples/tree_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl +++ b/examples/tree_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl @@ -6,10 +6,10 @@ using Trixi # TODO: parabolic; unify names of these accessor functions prandtl_number() = 0.72 -viscosity() = 0.001 +dynamic_viscosity() = 0.001 equations = CompressibleEulerEquations2D(1.4) -equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, viscosity=viscosity(), +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, viscosity=dynamic_viscosity(), Prandtl=prandtl_number()) # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux From f459934d0f0a636921ecd621039427f750d14273 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 7 Oct 2022 10:06:30 +0200 Subject: [PATCH 09/12] add comment about units in the docstring --- src/equations/compressible_navier_stokes_2d.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index f53bd9fd48f..a4d0f3df92e 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -13,6 +13,9 @@ the [`CompressibleEulerEquations2D`](@ref). - `gradient_variables`: which variables the gradients are taken with respect to. Defaults to `GradientVariablesPrimitive()`. +Fluid properties such as the dynamic viscosity mu can be provided in any unit system, e.g., +[``\mu``] = kg m⁻¹ s⁻¹. + The particular form of the compressible Navier-Stokes implemented is ```math \frac{\partial}{\partial t} From 507ebe8ae0a02bb572cfeeea01071bc4f002cdfc Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 7 Oct 2022 11:26:13 +0200 Subject: [PATCH 10/12] Update src/equations/compressible_navier_stokes_2d.jl Co-authored-by: Michael Schlottke-Lakemper --- src/equations/compressible_navier_stokes_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index a4d0f3df92e..bc3665ca8fb 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -13,7 +13,7 @@ the [`CompressibleEulerEquations2D`](@ref). - `gradient_variables`: which variables the gradients are taken with respect to. Defaults to `GradientVariablesPrimitive()`. -Fluid properties such as the dynamic viscosity mu can be provided in any unit system, e.g., +Fluid properties such as the dynamic viscosity ``\\mu`` can be provided in any consistent unit system, e.g., [``\mu``] = kg m⁻¹ s⁻¹. The particular form of the compressible Navier-Stokes implemented is From 55aadd87fde4505904967df251f6d640d8f5faf8 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 7 Oct 2022 13:09:18 +0200 Subject: [PATCH 11/12] unify notation to call the viscosity mu --- examples/dgmulti_2d/elixir_navierstokes_convergence.jl | 6 +++--- .../dgmulti_2d/elixir_navierstokes_lid_driven_cavity.jl | 4 ++-- examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl | 6 +++--- .../elixir_navierstokes_lid_driven_cavity.jl | 4 ++-- src/equations/compressible_navier_stokes_2d.jl | 8 ++++---- 5 files changed, 14 insertions(+), 14 deletions(-) diff --git a/examples/dgmulti_2d/elixir_navierstokes_convergence.jl b/examples/dgmulti_2d/elixir_navierstokes_convergence.jl index ecb545cf9e6..a602584663b 100644 --- a/examples/dgmulti_2d/elixir_navierstokes_convergence.jl +++ b/examples/dgmulti_2d/elixir_navierstokes_convergence.jl @@ -5,12 +5,12 @@ using Trixi # semidiscretization of the ideal compressible Navier-Stokes equations prandtl_number() = 0.72 -dynamic_viscosity() = 0.01 +mu() = 0.01 equations = CompressibleEulerEquations2D(1.4) # Note: If you change the Navier-Stokes parameters here, also change them in the initial condition # I really do not like this structure but it should work for now -equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, viscosity=dynamic_viscosity(), Prandtl=prandtl_number(), +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu=mu(), Prandtl=prandtl_number(), gradient_variables=GradientVariablesPrimitive()) # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux @@ -52,7 +52,7 @@ end # see also https://github.com/trixi-framework/Trixi.jl/pull/1160 inv_gamma_minus_one = inv(equations.gamma - 1) Pr = prandtl_number() - mu = dynamic_viscosity() + mu = mu() # Same settings as in `initial_condition` # Amplitude and shift diff --git a/examples/dgmulti_2d/elixir_navierstokes_lid_driven_cavity.jl b/examples/dgmulti_2d/elixir_navierstokes_lid_driven_cavity.jl index ed9008ec351..bd986422cc4 100644 --- a/examples/dgmulti_2d/elixir_navierstokes_lid_driven_cavity.jl +++ b/examples/dgmulti_2d/elixir_navierstokes_lid_driven_cavity.jl @@ -6,10 +6,10 @@ using Trixi # TODO: parabolic; unify names of these accessor functions prandtl_number() = 0.72 -dynamic_viscosity() = 0.001 +mu() = 0.001 equations = CompressibleEulerEquations2D(1.4) -equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, viscosity=dynamic_viscosity(), +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu=mu(), Prandtl=prandtl_number()) diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl b/examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl index 1c8960111c7..53e32bc11b9 100644 --- a/examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl +++ b/examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl @@ -5,10 +5,10 @@ using Trixi # semidiscretization of the ideal compressible Navier-Stokes equations prandtl_number() = 0.72 -dynamic_viscosity() = 0.01 +mu() = 0.01 equations = CompressibleEulerEquations2D(1.4) -equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, viscosity=dynamic_viscosity(), Prandtl=prandtl_number(), +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu=mu(), Prandtl=prandtl_number(), gradient_variables=GradientVariablesPrimitive()) # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux @@ -54,7 +54,7 @@ end # see also https://github.com/trixi-framework/Trixi.jl/pull/1160 inv_gamma_minus_one = inv(equations.gamma - 1) Pr = prandtl_number() - mu = dynamic_viscosity() + mu = mu() # Same settings as in `initial_condition` # Amplitude and shift 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 fda93ab87e2..10e6865e719 100644 --- a/examples/tree_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl +++ b/examples/tree_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl @@ -6,10 +6,10 @@ using Trixi # TODO: parabolic; unify names of these accessor functions prandtl_number() = 0.72 -dynamic_viscosity() = 0.001 +mu() = 0.001 equations = CompressibleEulerEquations2D(1.4) -equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, viscosity=dynamic_viscosity(), +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/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index bc3665ca8fb..9f7758aa476 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -13,7 +13,7 @@ the [`CompressibleEulerEquations2D`](@ref). - `gradient_variables`: which variables the gradients are taken with respect to. Defaults to `GradientVariablesPrimitive()`. -Fluid properties such as the dynamic viscosity ``\\mu`` can be provided in any consistent unit system, e.g., +Fluid properties such as the dynamic viscosity ``\mu`` can be provided in any consistent unit system, e.g., [``\mu``] = kg m⁻¹ s⁻¹. The particular form of the compressible Navier-Stokes implemented is @@ -111,11 +111,11 @@ struct GradientVariablesEntropy end # default to primitive gradient variables function CompressibleNavierStokesDiffusion2D(equations::CompressibleEulerEquations2D; - viscosity, Prandtl, + mu, Prandtl, gradient_variables = GradientVariablesPrimitive()) gamma = equations.gamma inv_gamma_minus_one = equations.inv_gamma_minus_one - mu, Pr = promote(viscosity, Prandtl) + μ, Pr = promote(mu, Prandtl) # Under the assumption of constant Prandtl number the thermal conductivity # constant is kappa = gamma μ / ((gamma-1) Pr). @@ -123,7 +123,7 @@ function CompressibleNavierStokesDiffusion2D(equations::CompressibleEulerEquatio kappa = gamma * inv_gamma_minus_one / Pr CompressibleNavierStokesDiffusion2D{typeof(gradient_variables), typeof(gamma), typeof(equations)}(gamma, inv_gamma_minus_one, - mu, Pr, kappa, + μ, Pr, kappa, equations, gradient_variables) end From 40ef1c34eac015f0f53d52f49f5c4fa668c9f026 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 7 Oct 2022 13:28:31 +0200 Subject: [PATCH 12/12] fix undefined variable --- .../elixir_navierstokes_convergence.jl | 54 +++++++++---------- .../elixir_navierstokes_convergence.jl | 54 +++++++++---------- 2 files changed, 54 insertions(+), 54 deletions(-) diff --git a/examples/dgmulti_2d/elixir_navierstokes_convergence.jl b/examples/dgmulti_2d/elixir_navierstokes_convergence.jl index a602584663b..f7039627238 100644 --- a/examples/dgmulti_2d/elixir_navierstokes_convergence.jl +++ b/examples/dgmulti_2d/elixir_navierstokes_convergence.jl @@ -52,7 +52,7 @@ end # see also https://github.com/trixi-framework/Trixi.jl/pull/1160 inv_gamma_minus_one = inv(equations.gamma - 1) Pr = prandtl_number() - mu = mu() + mu_ = mu() # Same settings as in `initial_condition` # Amplitude and shift @@ -117,10 +117,10 @@ end + rho * v1_y * v2 + rho * v1 * v2_y # stress tensor from x-direction - - 4.0 / 3.0 * v1_xx * mu - + 2.0 / 3.0 * v2_xy * mu - - v1_yy * mu - - v2_xy * mu ) + - 4.0 / 3.0 * v1_xx * mu_ + + 2.0 / 3.0 * v2_xy * mu_ + - v1_yy * mu_ + - v2_xy * mu_ ) # y-momentum equation du3 = ( rho_t * v2 + rho * v2_t + p_y + rho_x * v1 * v2 + rho * v1_x * v2 @@ -128,39 +128,39 @@ end + rho_y * v2^2 + 2.0 * rho * v2 * v2_y # stress tensor from y-direction - - v1_xy * mu - - v2_xx * mu - - 4.0 / 3.0 * v2_yy * mu - + 2.0 / 3.0 * v1_xy * mu ) + - v1_xy * mu_ + - v2_xx * mu_ + - 4.0 / 3.0 * v2_yy * mu_ + + 2.0 / 3.0 * v1_xy * mu_ ) # total energy equation du4 = ( E_t + v1_x * (E + p) + v1 * (E_x + p_x) + v2_y * (E + p) + v2 * (E_y + p_y) # stress tensor and temperature gradient terms from x-direction - - 4.0 / 3.0 * v1_xx * v1 * mu - + 2.0 / 3.0 * v2_xy * v1 * mu - - 4.0 / 3.0 * v1_x * v1_x * mu - + 2.0 / 3.0 * v2_y * v1_x * mu - - v1_xy * v2 * mu - - v2_xx * v2 * mu - - v1_y * v2_x * mu - - v2_x * v2_x * mu + - 4.0 / 3.0 * v1_xx * v1 * mu_ + + 2.0 / 3.0 * v2_xy * v1 * mu_ + - 4.0 / 3.0 * v1_x * v1_x * mu_ + + 2.0 / 3.0 * v2_y * v1_x * mu_ + - v1_xy * v2 * mu_ + - v2_xx * v2 * mu_ + - v1_y * v2_x * mu_ + - v2_x * v2_x * mu_ - T_const * inv_rho_cubed * ( p_xx * rho * rho - 2.0 * p_x * rho * rho_x + 2.0 * p * rho_x * rho_x - - p * rho * rho_xx ) * mu + - p * rho * rho_xx ) * mu_ # stress tensor and temperature gradient terms from y-direction - - v1_yy * v1 * mu - - v2_xy * v1 * mu - - v1_y * v1_y * mu - - v2_x * v1_y * mu - - 4.0 / 3.0 * v2_yy * v2 * mu - + 2.0 / 3.0 * v1_xy * v2 * mu - - 4.0 / 3.0 * v2_y * v2_y * mu - + 2.0 / 3.0 * v1_x * v2_y * mu + - v1_yy * v1 * mu_ + - v2_xy * v1 * mu_ + - v1_y * v1_y * mu_ + - v2_x * v1_y * mu_ + - 4.0 / 3.0 * v2_yy * v2 * mu_ + + 2.0 / 3.0 * v1_xy * v2 * mu_ + - 4.0 / 3.0 * v2_y * v2_y * mu_ + + 2.0 / 3.0 * v1_x * v2_y * mu_ - T_const * inv_rho_cubed * ( p_yy * rho * rho - 2.0 * p_y * rho * rho_y + 2.0 * p * rho_y * rho_y - - p * rho * rho_yy ) * mu ) + - p * rho * rho_yy ) * mu_ ) return SVector(du1, du2, du3, du4) end diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl b/examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl index 53e32bc11b9..1145c22fcaf 100644 --- a/examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl +++ b/examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl @@ -54,7 +54,7 @@ end # see also https://github.com/trixi-framework/Trixi.jl/pull/1160 inv_gamma_minus_one = inv(equations.gamma - 1) Pr = prandtl_number() - mu = mu() + mu_ = mu() # Same settings as in `initial_condition` # Amplitude and shift @@ -119,10 +119,10 @@ end + rho * v1_y * v2 + rho * v1 * v2_y # stress tensor from x-direction - - 4.0 / 3.0 * v1_xx * mu - + 2.0 / 3.0 * v2_xy * mu - - v1_yy * mu - - v2_xy * mu ) + - 4.0 / 3.0 * v1_xx * mu_ + + 2.0 / 3.0 * v2_xy * mu_ + - v1_yy * mu_ + - v2_xy * mu_ ) # y-momentum equation du3 = ( rho_t * v2 + rho * v2_t + p_y + rho_x * v1 * v2 + rho * v1_x * v2 @@ -130,39 +130,39 @@ end + rho_y * v2^2 + 2.0 * rho * v2 * v2_y # stress tensor from y-direction - - v1_xy * mu - - v2_xx * mu - - 4.0 / 3.0 * v2_yy * mu - + 2.0 / 3.0 * v1_xy * mu ) + - v1_xy * mu_ + - v2_xx * mu_ + - 4.0 / 3.0 * v2_yy * mu_ + + 2.0 / 3.0 * v1_xy * mu_ ) # total energy equation du4 = ( E_t + v1_x * (E + p) + v1 * (E_x + p_x) + v2_y * (E + p) + v2 * (E_y + p_y) # stress tensor and temperature gradient terms from x-direction - - 4.0 / 3.0 * v1_xx * v1 * mu - + 2.0 / 3.0 * v2_xy * v1 * mu - - 4.0 / 3.0 * v1_x * v1_x * mu - + 2.0 / 3.0 * v2_y * v1_x * mu - - v1_xy * v2 * mu - - v2_xx * v2 * mu - - v1_y * v2_x * mu - - v2_x * v2_x * mu + - 4.0 / 3.0 * v1_xx * v1 * mu_ + + 2.0 / 3.0 * v2_xy * v1 * mu_ + - 4.0 / 3.0 * v1_x * v1_x * mu_ + + 2.0 / 3.0 * v2_y * v1_x * mu_ + - v1_xy * v2 * mu_ + - v2_xx * v2 * mu_ + - v1_y * v2_x * mu_ + - v2_x * v2_x * mu_ - T_const * inv_rho_cubed * ( p_xx * rho * rho - 2.0 * p_x * rho * rho_x + 2.0 * p * rho_x * rho_x - - p * rho * rho_xx ) * mu + - p * rho * rho_xx ) * mu_ # stress tensor and temperature gradient terms from y-direction - - v1_yy * v1 * mu - - v2_xy * v1 * mu - - v1_y * v1_y * mu - - v2_x * v1_y * mu - - 4.0 / 3.0 * v2_yy * v2 * mu - + 2.0 / 3.0 * v1_xy * v2 * mu - - 4.0 / 3.0 * v2_y * v2_y * mu - + 2.0 / 3.0 * v1_x * v2_y * mu + - v1_yy * v1 * mu_ + - v2_xy * v1 * mu_ + - v1_y * v1_y * mu_ + - v2_x * v1_y * mu_ + - 4.0 / 3.0 * v2_yy * v2 * mu_ + + 2.0 / 3.0 * v1_xy * v2 * mu_ + - 4.0 / 3.0 * v2_y * v2_y * mu_ + + 2.0 / 3.0 * v1_x * v2_y * mu_ - T_const * inv_rho_cubed * ( p_yy * rho * rho - 2.0 * p_y * rho * rho_y + 2.0 * p * rho_y * rho_y - - p * rho * rho_yy ) * mu ) + - p * rho * rho_yy ) * mu_ ) return SVector(du1, du2, du3, du4) end