diff --git a/examples/dgmulti_2d/elixir_navier_stokes_convergence.jl b/examples/dgmulti_2d/elixir_navierstokes_convergence.jl similarity index 81% rename from examples/dgmulti_2d/elixir_navier_stokes_convergence.jl rename to examples/dgmulti_2d/elixir_navierstokes_convergence.jl index 6d64f662e77..f7039627238 100644 --- a/examples/dgmulti_2d/elixir_navier_stokes_convergence.jl +++ b/examples/dgmulti_2d/elixir_navierstokes_convergence.jl @@ -4,14 +4,14 @@ using Trixi ############################################################################### # semidiscretization of the ideal compressible Navier-Stokes equations -reynolds_number() = 100 prandtl_number() = 0.72 +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, Reynolds=reynolds_number(), Prandtl=prandtl_number(), - Mach_freestream=0.5, gradient_variables=GradientVariablesPrimitive()) +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 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_ = mu() # 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_navierstokes_lid_driven_cavity.jl similarity index 93% rename from examples/dgmulti_2d/elixir_navier_stokes_lid_driven_cavity.jl rename to examples/dgmulti_2d/elixir_navierstokes_lid_driven_cavity.jl index be8ad1a0730..bd986422cc4 100644 --- a/examples/dgmulti_2d/elixir_navier_stokes_lid_driven_cavity.jl +++ b/examples/dgmulti_2d/elixir_navierstokes_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 +mu() = 0.001 equations = CompressibleEulerEquations2D(1.4) -equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, Reynolds=reynolds_number(), - Prandtl=prandtl_number(), - Mach_freestream=mach_number()) +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 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/examples/tree_2d_dgsem/elixir_navier_stokes_convergence.jl b/examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl similarity index 82% rename from examples/tree_2d_dgsem/elixir_navier_stokes_convergence.jl rename to examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl index bd74271dcac..1145c22fcaf 100644 --- a/examples/tree_2d_dgsem/elixir_navier_stokes_convergence.jl +++ b/examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl @@ -4,12 +4,12 @@ using Trixi ############################################################################### # semidiscretization of the ideal compressible Navier-Stokes equations -reynolds_number() = 100 prandtl_number() = 0.72 +mu() = 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, mu=mu(), 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_ = mu() # 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_navierstokes_lid_driven_cavity.jl similarity index 93% 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 index 89882675613..10e6865e719 100644 --- a/examples/tree_2d_dgsem/elixir_navier_stokes_lid_driven_cavity.jl +++ b/examples/tree_2d_dgsem/elixir_navierstokes_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 +mu() = 0.001 equations = CompressibleEulerEquations2D(1.4) -equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, Reynolds=reynolds_number(), - Prandtl=prandtl_number(), - Mach_freestream=mach_number()) +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 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/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..9f7758aa476 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, equations, gradient_variables=GradientVariablesPrimitive()) These equations contain the diffusion (i.e. parabolic) terms applied @@ -7,13 +7,15 @@ to mass, momenta, and total energy together with the advective terms from the [`CompressibleEulerEquations2D`](@ref). - `gamma`: adiabatic constant, -- `Re`: Reynolds number, +- `mu`: dynamic viscosity, - `Pr`: Prandtl number, -- `Ma_inf`: free-stream Mach number - `equations`: instance of 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., +[``\mu``] = kg m⁻¹ s⁻¹. + The particular form of the compressible Navier-Stokes implemented is ```math \frac{\partial}{\partial t} @@ -45,11 +47,17 @@ 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}}. +``` +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 @@ -65,24 +73,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 +82,11 @@ 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 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 +111,19 @@ struct GradientVariablesEntropy end # default to primitive gradient variables function CompressibleNavierStokesDiffusion2D(equations::CompressibleEulerEquations2D; - Reynolds, Prandtl, Mach_freestream, + mu, Prandtl, gradient_variables = GradientVariablesPrimitive()) gamma = equations.gamma inv_gamma_minus_one = equations.inv_gamma_minus_one - Re, Pr, Ma = promote(Reynolds, Prandtl, Mach_freestream) + μ, Pr = promote(mu, Prandtl) # 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 + # constant is kappa = gamma μ / ((gamma-1) Pr). + # Important note! Factor of μ is accounted for later in `flux`. + kappa = gamma * 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, + μ, Pr, kappa, equations, gradient_variables) end @@ -159,15 +140,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) @@ -187,14 +165,17 @@ 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) - # Important note! Due to nondimensional scaling R = 1 / gamma, so the - # temperature T in the gradient computation already contains a factor of gamma + # 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 - # kinematic viscosity is simply 1/Re for this nondimensionalization - mu = 1.0 / equations.Re + # 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 if orientation == 1 # viscous flux components in the x-direction @@ -272,9 +253,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 @@ -291,7 +272,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 @@ -414,8 +395,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 diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index e0552217b3e..d9c3967315a 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -94,19 +94,19 @@ 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.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"), + @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.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 @@ -126,54 +126,54 @@ 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.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 - @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.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] + @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"), - initial_refinement_level=2, tspan=(0,.1), gradient_variables=GradientVariablesEntropy(), + @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.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 - @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.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"), + @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.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