From d52a0419f0fe5dbee8dd58e3b12a23bc9fa67fc1 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Fri, 11 Aug 2023 09:43:31 +0200 Subject: [PATCH] Navier-Stokes 1D (#1597) * Remove doubled implementations * kepp main updated with true main * Avoid allocations in parabolic boundary fluxes * Correct shear layer IC * Whitespaces * Restore main * restore main * 1D Navier Stokes * Conventional notation for heat flux * remove multi-dim artefacts * Move general part into own file * Slip Wall BC for 1D Compressible Euler * Correct arguments for 1D BCs * format * Add convergence test with walls * Test gradient with entropy variables * Test isothermal BC, test gradient in entropy vars * Correct test data --------- Co-authored-by: Hendrik Ranocha --- ...lixir_navierstokes_convergence_periodic.jl | 136 ++++++ .../elixir_navierstokes_convergence_walls.jl | 160 +++++++ src/Trixi.jl | 3 +- src/equations/compressible_euler_1d.jl | 51 +++ src/equations/compressible_navier_stokes.jl | 70 +++ .../compressible_navier_stokes_1d.jl | 403 ++++++++++++++++++ .../compressible_navier_stokes_2d.jl | 77 +--- .../compressible_navier_stokes_3d.jl | 6 +- src/equations/equations_parabolic.jl | 2 + test/test_parabolic_1d.jl | 35 +- 10 files changed, 864 insertions(+), 79 deletions(-) create mode 100644 examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic.jl create mode 100644 examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls.jl create mode 100644 src/equations/compressible_navier_stokes.jl create mode 100644 src/equations/compressible_navier_stokes_1d.jl diff --git a/examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic.jl b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic.jl new file mode 100644 index 00000000000..3f72d319b0b --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic.jl @@ -0,0 +1,136 @@ + +using OrdinaryDiffEq +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 + +equations = CompressibleEulerEquations1D(1.4) +equations_parabolic = CompressibleNavierStokesDiffusion1D(equations, mu=mu(), + Prandtl=prandtl_number()) + +# This convergence test setup was originally derived by Andrew Winters (@andrewwinters5000) +# (Simplified version of the 2D) +function initial_condition_navier_stokes_convergence_test(x, t, equations) + # Amplitude and shift + A = 0.5 + c = 2.0 + + # convenience values for trig. functions + pi_x = pi * x[1] + pi_t = pi * t + + rho = c + A * sin(pi_x) * cos(pi_t) + v1 = sin(pi_x) * cos(pi_t) + p = rho^2 + + return prim2cons(SVector(rho, v1, p), equations) +end +initial_condition = initial_condition_navier_stokes_convergence_test + +@inline function source_terms_navier_stokes_convergence_test(u, x, t, equations) + # we currently need to hardcode these parameters until we fix the "combined equation" issue + # see also https://github.com/trixi-framework/Trixi.jl/pull/1160 + inv_gamma_minus_one = inv(equations.gamma - 1) + Pr = prandtl_number() + mu_ = mu() + + # Same settings as in `initial_condition` + # Amplitude and shift + A = 0.5 + c = 2.0 + + # convenience values for trig. functions + pi_x = pi * x[1] + pi_t = pi * t + + # compute the manufactured solution and all necessary derivatives + rho = c + A * sin(pi_x) * cos(pi_t) + rho_t = -pi * A * sin(pi_x) * sin(pi_t) + rho_x = pi * A * cos(pi_x) * cos(pi_t) + rho_xx = -pi * pi * A * sin(pi_x) * cos(pi_t) + + v1 = sin(pi_x) * cos(pi_t) + v1_t = -pi * sin(pi_x) * sin(pi_t) + v1_x = pi * cos(pi_x) * cos(pi_t) + v1_xx = -pi * pi * sin(pi_x) * cos(pi_t) + + p = rho * rho + p_t = 2.0 * rho * rho_t + p_x = 2.0 * rho * rho_x + p_xx = 2.0 * rho * rho_xx + 2.0 * rho_x * rho_x + + E = p * inv_gamma_minus_one + 0.5 * rho * v1^2 + E_t = p_t * inv_gamma_minus_one + 0.5 * rho_t * v1^2 + rho * v1 * v1_t + E_x = p_x * inv_gamma_minus_one + 0.5 * rho_x * v1^2 + rho * v1 * v1_x + + # Some convenience constants + T_const = equations.gamma * inv_gamma_minus_one / Pr + inv_rho_cubed = 1.0 / (rho^3) + + # compute the source terms + # density equation + du1 = rho_t + rho_x * v1 + rho * v1_x + + # x-momentum equation + du2 = ( rho_t * v1 + rho * v1_t + + p_x + rho_x * v1^2 + 2.0 * rho * v1 * v1_x + # stress tensor from x-direction + - v1_xx * mu_) + + # total energy equation + du3 = ( E_t + v1_x * (E + p) + v1 * (E_x + p_x) + # stress tensor and temperature gradient terms from x-direction + - v1_xx * v1 * mu_ + - v1_x * v1_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_) + + return SVector(du1, du2, du3) +end + +volume_flux = flux_ranocha +solver = DGSEM(polydeg=3, surface_flux=flux_hllc, + volume_integral=VolumeIntegralFluxDifferencing(volume_flux)) + +coordinates_min = -1.0 +coordinates_max = 1.0 +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level=4, + n_cells_max=100_000) + + +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), + initial_condition, solver, + source_terms = source_terms_navier_stokes_convergence_test) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 10.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval) + +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 \ No newline at end of file diff --git a/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls.jl b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls.jl new file mode 100644 index 00000000000..181a2cb209f --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls.jl @@ -0,0 +1,160 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the ideal compressible Navier-Stokes equations + +prandtl_number() = 0.72 +mu() = 0.01 + +equations = CompressibleEulerEquations1D(1.4) +equations_parabolic = CompressibleNavierStokesDiffusion1D(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, + volume_integral=VolumeIntegralWeakForm()) + +coordinates_min = -1.0 +coordinates_max = 1.0 + +# Create a uniformly refined mesh with periodic boundaries +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level=3, + periodicity=false, + n_cells_max=30_000) # set maximum capacity of tree data structure + +# Note: the initial condition cannot be specialized to `CompressibleNavierStokesDiffusion1D` +# since it is called by both the parabolic solver (which passes in `CompressibleNavierStokesDiffusion1D`) +# and by the initial condition (which passes in `CompressibleEulerEquations1D`). +# This convergence test setup was originally derived by Andrew Winters (@andrewwinters5000) +function initial_condition_navier_stokes_convergence_test(x, t, equations) + # Amplitude and shift + A = 0.5 + c = 2.0 + + # convenience values for trig. functions + pi_x = pi * x[1] + pi_t = pi * t + + rho = c + A * cos(pi_x) * cos(pi_t) + v1 = log(x[1] + 2.0) * (1.0 - exp(-A * (x[1] - 1.0)) ) * cos(pi_t) + p = rho^2 + + return prim2cons(SVector(rho, v1, p), equations) +end + +@inline function source_terms_navier_stokes_convergence_test(u, x, t, equations) + x = x[1] + + # TODO: parabolic + # we currently need to hardcode these parameters until we fix the "combined equation" issue + # see also https://github.com/trixi-framework/Trixi.jl/pull/1160 + inv_gamma_minus_one = inv(equations.gamma - 1) + Pr = prandtl_number() + mu_ = mu() + + # Same settings as in `initial_condition` + # Amplitude and shift + A = 0.5 + c = 2.0 + + # convenience values for trig. functions + pi_x = pi * x + pi_t = pi * t + + # compute the manufactured solution and all necessary derivatives + rho = c + A * cos(pi_x) * cos(pi_t) + rho_t = -pi * A * cos(pi_x) * sin(pi_t) + rho_x = -pi * A * sin(pi_x) * cos(pi_t) + rho_xx = -pi * pi * A * cos(pi_x) * cos(pi_t) + + v1 = log(x + 2.0) * (1.0 - exp(-A * (x - 1.0))) * cos(pi_t) + v1_t = -pi * log(x + 2.0) * (1.0 - exp(-A * (x - 1.0))) * sin(pi_t) + v1_x = (A * log(x + 2.0) * exp(-A * (x - 1.0)) + (1.0 - exp(-A * (x - 1.0))) / (x + 2.0)) * cos(pi_t) + v1_xx = (( 2.0 * A * exp(-A * (x - 1.0)) / (x + 2.0) + - A * A * log(x + 2.0) * exp(-A * (x - 1.0)) + - (1.0 - exp(-A * (x - 1.0))) / ((x + 2.0) * (x + 2.0))) * cos(pi_t)) + + p = rho * rho + p_t = 2.0 * rho * rho_t + p_x = 2.0 * rho * rho_x + p_xx = 2.0 * rho * rho_xx + 2.0 * rho_x * rho_x + + # Note this simplifies slightly because the ansatz assumes that v1 = v2 + E = p * inv_gamma_minus_one + 0.5 * rho * v1^2 + E_t = p_t * inv_gamma_minus_one + 0.5 * rho_t * v1^2 + rho * v1 * v1_t + E_x = p_x * inv_gamma_minus_one + 0.5 * rho_x * v1^2 + rho * v1 * v1_x + + # Some convenience constants + T_const = equations.gamma * inv_gamma_minus_one / Pr + inv_rho_cubed = 1.0 / (rho^3) + + # compute the source terms + # density equation + du1 = rho_t + rho_x * v1 + rho * v1_x + + # y-momentum equation + du2 = ( rho_t * v1 + rho * v1_t + + p_x + rho_x * v1^2 + 2.0 * rho * v1 * v1_x + # stress tensor from y-direction + - v1_xx * mu_) + + # total energy equation + du3 = ( E_t + v1_x * (E + p) + v1 * (E_x + p_x) + # stress tensor and temperature gradient terms from x-direction + - v1_xx * v1 * mu_ + - v1_x * v1_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_ ) + + return SVector(du1, du2, du3) +end + +initial_condition = initial_condition_navier_stokes_convergence_test + +# BC types +velocity_bc_left_right = NoSlip((x, t, equations) -> initial_condition_navier_stokes_convergence_test(x, t, equations)[2]) + +heat_bc_left = Isothermal((x, t, equations) -> + Trixi.temperature(initial_condition_navier_stokes_convergence_test(x, t, equations), + equations_parabolic)) +heat_bc_right = Adiabatic((x, t, equations) -> 0.0) + +boundary_condition_left = BoundaryConditionNavierStokesWall(velocity_bc_left_right, heat_bc_left) +boundary_condition_right = BoundaryConditionNavierStokesWall(velocity_bc_left_right, heat_bc_right) + +# define inviscid boundary conditions +boundary_conditions = (; x_neg = boundary_condition_slip_wall, + x_pos = boundary_condition_slip_wall) + +# define viscous boundary conditions +boundary_conditions_parabolic = (; x_neg = boundary_condition_left, + x_pos = boundary_condition_right) + +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver; + boundary_conditions=(boundary_conditions, boundary_conditions_parabolic), + source_terms=source_terms_navier_stokes_convergence_test) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span `tspan` +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() +alive_callback = AliveCallback(alive_interval=10) +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval) +callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback) + +############################################################################### +# run the simulation + +time_int_tol = 1e-8 +sol = solve(ode, RDPK3SpFSAL49(); abstol=time_int_tol, reltol=time_int_tol, dt = 1e-5, + ode_default_options()..., callback=callbacks) +summary_callback() # print the timer summary \ No newline at end of file diff --git a/src/Trixi.jl b/src/Trixi.jl index 990c33f3c94..78ddaa3ca7f 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -152,7 +152,8 @@ export AcousticPerturbationEquations2D, LinearizedEulerEquations2D export LaplaceDiffusion1D, LaplaceDiffusion2D, - CompressibleNavierStokesDiffusion2D, CompressibleNavierStokesDiffusion3D + CompressibleNavierStokesDiffusion1D, CompressibleNavierStokesDiffusion2D, + CompressibleNavierStokesDiffusion3D export GradientVariablesPrimitive, GradientVariablesEntropy diff --git a/src/equations/compressible_euler_1d.jl b/src/equations/compressible_euler_1d.jl index e4fd0997eae..9204989e8be 100644 --- a/src/equations/compressible_euler_1d.jl +++ b/src/equations/compressible_euler_1d.jl @@ -198,6 +198,57 @@ function initial_condition_eoc_test_coupled_euler_gravity(x, t, return prim2cons(SVector(rho, v1, p), equations) end +""" + boundary_condition_slip_wall(u_inner, orientation, direction, x, t, + surface_flux_function, equations::CompressibleEulerEquations1D) +Determine the boundary numerical surface flux for a slip wall condition. +Imposes a zero normal velocity at the wall. +Density is taken from the internal solution state and pressure is computed as an +exact solution of a 1D Riemann problem. Further details about this boundary state +are available in the paper: +- J. J. W. van der Vegt and H. van der Ven (2002) + Slip flow boundary conditions in discontinuous Galerkin discretizations of + the Euler equations of gas dynamics + [PDF](https://reports.nlr.nl/bitstream/handle/10921/692/TP-2002-300.pdf?sequence=1) + + Should be used together with [`TreeMesh`](@ref). +""" +@inline function boundary_condition_slip_wall(u_inner, orientation, + direction, x, t, + surface_flux_function, + equations::CompressibleEulerEquations1D) + # compute the primitive variables + rho_local, v_normal, p_local = cons2prim(u_inner, equations) + + if isodd(direction) # flip sign of normal to make it outward pointing + v_normal *= -1 + end + + # Get the solution of the pressure Riemann problem + # See Section 6.3.3 of + # Eleuterio F. Toro (2009) + # Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction + # [DOI: 10.1007/b79761](https://doi.org/10.1007/b79761) + if v_normal <= 0.0 + sound_speed = sqrt(equations.gamma * p_local / rho_local) # local sound speed + p_star = p_local * + (1 + 0.5 * (equations.gamma - 1) * v_normal / sound_speed)^(2 * + equations.gamma * + equations.inv_gamma_minus_one) + else # v_normal > 0.0 + A = 2 / ((equations.gamma + 1) * rho_local) + B = p_local * (equations.gamma - 1) / (equations.gamma + 1) + p_star = p_local + + 0.5 * v_normal / A * + (v_normal + sqrt(v_normal^2 + 4 * A * (p_local + B))) + end + + # For the slip wall we directly set the flux as the normal velocity is zero + return SVector(zero(eltype(u_inner)), + p_star, + zero(eltype(u_inner))) +end + # Calculate 1D flux for a single point @inline function flux(u, orientation::Integer, equations::CompressibleEulerEquations1D) rho, rho_v1, rho_e = u diff --git a/src/equations/compressible_navier_stokes.jl b/src/equations/compressible_navier_stokes.jl new file mode 100644 index 00000000000..af7897d4586 --- /dev/null +++ b/src/equations/compressible_navier_stokes.jl @@ -0,0 +1,70 @@ +# TODO: can we generalize this to MHD? +""" + struct BoundaryConditionNavierStokesWall + +Creates a wall-type boundary conditions for the compressible Navier-Stokes equations. +The fields `boundary_condition_velocity` and `boundary_condition_heat_flux` are intended +to be boundary condition types such as the `NoSlip` velocity boundary condition and the +`Adiabatic` or `Isothermal` heat boundary condition. + +!!! warning "Experimental feature" + This is an experimental feature and may change in future releases. +""" +struct BoundaryConditionNavierStokesWall{V, H} + boundary_condition_velocity::V + boundary_condition_heat_flux::H +end + +""" + struct NoSlip + +Use to create a no-slip boundary condition with `BoundaryConditionNavierStokesWall`. The field `boundary_value_function` +should be a function with signature `boundary_value_function(x, t, equations)` +and should return a `SVector{NDIMS}` whose entries are the velocity vector at a +point `x` and time `t`. +""" +struct NoSlip{F} + boundary_value_function::F # value of the velocity vector on the boundary +end + +""" + struct Isothermal + +Used to create a no-slip boundary condition with [`BoundaryConditionNavierStokesWall`](@ref). +The field `boundary_value_function` should be a function with signature +`boundary_value_function(x, t, equations)` and return a scalar value for the +temperature at point `x` and time `t`. +""" +struct Isothermal{F} + boundary_value_function::F # value of the temperature on the boundary +end + +""" + struct Adiabatic + +Used to create a no-slip boundary condition with [`BoundaryConditionNavierStokesWall`](@ref). +The field `boundary_value_normal_flux_function` should be a function with signature +`boundary_value_normal_flux_function(x, t, equations)` and return a scalar value for the +normal heat flux at point `x` and time `t`. +""" +struct Adiabatic{F} + boundary_value_normal_flux_function::F # scaled heat flux 1/T * kappa * dT/dn +end + +""" +!!! warning "Experimental code" + This code is experimental and may be changed or removed in any future release. + +`GradientVariablesPrimitive` and `GradientVariablesEntropy` are gradient variable type parameters +for `CompressibleNavierStokesDiffusion1D`. By default, the gradient variables are set to be +`GradientVariablesPrimitive`. Specifying `GradientVariablesEntropy` instead uses the entropy variable +formulation from +- Hughes, Mallet, Franca (1986) + A new finite element formulation for computational fluid dynamics: I. Symmetric forms of the + compressible Euler and Navier-Stokes equations and the second law of thermodynamics. + [https://doi.org/10.1016/0045-7825(86)90127-1](https://doi.org/10.1016/0045-7825(86)90127-1) + +Under `GradientVariablesEntropy`, the Navier-Stokes discretization is provably entropy stable. +""" +struct GradientVariablesPrimitive end +struct GradientVariablesEntropy end diff --git a/src/equations/compressible_navier_stokes_1d.jl b/src/equations/compressible_navier_stokes_1d.jl new file mode 100644 index 00000000000..dca846cac1e --- /dev/null +++ b/src/equations/compressible_navier_stokes_1d.jl @@ -0,0 +1,403 @@ +@doc raw""" + CompressibleNavierStokesDiffusion1D(equations; mu, Pr, + gradient_variables=GradientVariablesPrimitive()) + +Contains the diffusion (i.e. parabolic) terms applied +to mass, momenta, and total energy together with the advective terms from +the [`CompressibleEulerEquations1D`](@ref). + +- `equations`: instance of the [`CompressibleEulerEquations1D`](@ref) +- `mu`: dynamic viscosity, +- `Pr`: Prandtl number, +- `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} +\begin{pmatrix} +\rho \\ \rho v \\ \rho e +\end{pmatrix} ++ +\frac{\partial}{\partial x} +\begin{pmatrix} + \rho v \\ \rho v^2 + p \\ (\rho e + p) v +\end{pmatrix} += +\frac{\partial}{\partial x} +\begin{pmatrix} +0 \\ \tau \\ \tau v - q +\end{pmatrix} +``` +where the system is closed with the ideal gas assumption giving +```math +p = (\gamma - 1) \left( \rho e - \frac{1}{2} \rho v^2 \right) +``` +as the pressure. The value of the adiabatic constant `gamma` is taken from the [`CompressibleEulerEquations1D`](@ref). +The terms on the right hand side of the system above +are built from the viscous stress +```math +\tau = \mu \frac{\partial}{\partial x} v +``` +where the heat flux is +```math +q = -\kappa \frac{\partial}{\partial x} \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, +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 +q = -\kappa \frac{\partial}{\partial x} \left(T\right) = -\frac{\gamma \mu}{(\gamma - 1)\textrm{Pr}} \frac{\partial}{\partial x} \left(\frac{p}{\rho}\right) +``` +which is the form implemented below in the [`flux`](@ref) function. + +In one spatial dimensions we require gradients for two quantities, e.g., +primitive quantities +```math +\frac{\partial}{\partial x} v,\, \frac{\partial}{\partial x} T +``` +or the entropy variables +```math +\frac{\partial}{\partial x} w_2,\, \frac{\partial}{\partial x} w_3 +``` +where +```math +w_2 = \frac{\rho v1}{p},\, w_3 = -\frac{\rho}{p} +``` + +!!! warning "Experimental code" + This code is experimental and may be changed or removed in any future release. +""" +struct CompressibleNavierStokesDiffusion1D{GradientVariables, RealT <: Real, + E <: AbstractCompressibleEulerEquations{1}} <: + AbstractCompressibleNavierStokesDiffusion{1, 3} + # TODO: parabolic + # 1) For now save gamma and inv(gamma-1) again, but could potentially reuse them from the Euler equations + # 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 + + mu::RealT # viscosity + Pr::RealT # Prandtl number + kappa::RealT # thermal diffusivity for Fick's law + + equations_hyperbolic::E # CompressibleEulerEquations1D + gradient_variables::GradientVariables # GradientVariablesPrimitive or GradientVariablesEntropy +end + +# default to primitive gradient variables +function CompressibleNavierStokesDiffusion1D(equations::CompressibleEulerEquations1D; + mu, Prandtl, + 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). + # Important note! Factor of μ is accounted for later in `flux`. + kappa = gamma * inv_gamma_minus_one / Pr + + CompressibleNavierStokesDiffusion1D{typeof(gradient_variables), typeof(gamma), + typeof(equations)}(gamma, inv_gamma_minus_one, + μ, Pr, kappa, + equations, gradient_variables) +end + +# TODO: parabolic +# This is the flexibility a user should have to select the different gradient variable types +# varnames(::typeof(cons2prim) , ::CompressibleNavierStokesDiffusion1D) = ("v1", "v2", "T") +# varnames(::typeof(cons2entropy), ::CompressibleNavierStokesDiffusion1D) = ("w2", "w3", "w4") + +function varnames(variable_mapping, + equations_parabolic::CompressibleNavierStokesDiffusion1D) + varnames(variable_mapping, equations_parabolic.equations_hyperbolic) +end + +# we specialize this function to compute gradients of primitive variables instead of +# conservative variables. +function gradient_variable_transformation(::CompressibleNavierStokesDiffusion1D{ + GradientVariablesPrimitive + }) + cons2prim +end +function gradient_variable_transformation(::CompressibleNavierStokesDiffusion1D{ + GradientVariablesEntropy + }) + cons2entropy +end + +# Explicit formulas for the diffusive Navier-Stokes fluxes are available, e.g., in Section 2 +# of the paper by Rueda-Ramírez, 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::CompressibleNavierStokesDiffusion1D) + # Here, `u` is assumed to be the "transformed" variables specified by `gradient_variable_transformation`. + rho, v1, _ = convert_transformed_to_primitive(u, equations) + # Here `gradients` is assumed to contain the gradients of the primitive variables (rho, v1, v2, T) + # either computed directly or reverse engineered from the gradient of the entropy variables + # by way of the `convert_gradient_variables` function. + _, dv1dx, dTdx = convert_derivative_to_primitive(u, gradients, equations) + + # Viscous stress (tensor) + tau_11 = dv1dx + + # 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 + + # 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 + + # viscous flux components in the x-direction + f1 = zero(rho) + f2 = tau_11 * mu + f3 = (v1 * tau_11 + q1) * mu + + return SVector(f1, f2, f3) +end + +# Convert conservative variables to primitive +@inline function cons2prim(u, equations::CompressibleNavierStokesDiffusion1D) + rho, rho_v1, _ = u + + v1 = rho_v1 / rho + T = temperature(u, equations) + + return SVector(rho, v1, T) +end + +# Convert conservative variables to entropy +# TODO: parabolic. We can improve efficiency by not computing w_1, which involves logarithms +# This can be done by specializing `cons2entropy` and `entropy2cons` to `CompressibleNavierStokesDiffusion1D`, +# but this may be confusing to new users. +function cons2entropy(u, equations::CompressibleNavierStokesDiffusion1D) + cons2entropy(u, equations.equations_hyperbolic) +end +function entropy2cons(w, equations::CompressibleNavierStokesDiffusion1D) + entropy2cons(w, equations.equations_hyperbolic) +end + +# the `flux` function takes in transformed variables `u` which depend on the type of the gradient variables. +# For CNS, it is simplest to formulate the viscous terms in primitive variables, so we transform the transformed +# variables into primitive variables. +@inline function convert_transformed_to_primitive(u_transformed, + equations::CompressibleNavierStokesDiffusion1D{ + GradientVariablesPrimitive + }) + return u_transformed +end + +# TODO: parabolic. Make this more efficient! +@inline function convert_transformed_to_primitive(u_transformed, + equations::CompressibleNavierStokesDiffusion1D{ + GradientVariablesEntropy + }) + # note: this uses CompressibleNavierStokesDiffusion1D versions of cons2prim and entropy2cons + return cons2prim(entropy2cons(u_transformed, equations), equations) +end + +# Takes the solution values `u` and gradient of the entropy variables (w_2, w_3, w_4) and +# reverse engineers the gradients to be terms of the primitive variables (v1, v2, T). +# Helpful because then the diffusive fluxes have the same form as on paper. +# Note, the first component of `gradient_entropy_vars` contains gradient(rho) which is unused. +# TODO: parabolic; entropy stable viscous terms +@inline function convert_derivative_to_primitive(u, gradient, + ::CompressibleNavierStokesDiffusion1D{ + GradientVariablesPrimitive + }) + return gradient +end + +# the first argument is always the "transformed" variables. +@inline function convert_derivative_to_primitive(w, gradient_entropy_vars, + equations::CompressibleNavierStokesDiffusion1D{ + GradientVariablesEntropy + }) + + # TODO: parabolic. This is inefficient to pass in transformed variables but then transform them back. + # We can fix this if we directly compute v1, v2, T from the entropy variables + u = entropy2cons(w, equations) # calls a "modified" entropy2cons defined for CompressibleNavierStokesDiffusion1D + rho, rho_v1, _ = u + + v1 = rho_v1 / rho + T = temperature(u, equations) + + return SVector(gradient_entropy_vars[1], + T * (gradient_entropy_vars[2] + v1 * gradient_entropy_vars[3]), # grad(u) = T*(grad(w_2)+v1*grad(w_3)) + T * T * gradient_entropy_vars[3]) +end + +# This routine is required because `prim2cons` is called in `initial_condition`, which +# is called with `equations::CompressibleEulerEquations1D`. This means it is inconsistent +# with `cons2prim(..., ::CompressibleNavierStokesDiffusion1D)` as defined above. +# TODO: parabolic. Is there a way to clean this up? +@inline function prim2cons(u, equations::CompressibleNavierStokesDiffusion1D) + prim2cons(u, equations.equations_hyperbolic) +end + +@inline function temperature(u, equations::CompressibleNavierStokesDiffusion1D) + rho, rho_v1, rho_e = u + + p = (equations.gamma - 1) * (rho_e - 0.5 * rho_v1^2 / rho) + T = p / rho + return T +end + +@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, + <:Adiabatic})(flux_inner, + u_inner, + orientation::Integer, + direction, + x, t, + operator_type::Gradient, + equations::CompressibleNavierStokesDiffusion1D{ + GradientVariablesPrimitive + }) + v1 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, + equations) + return SVector(u_inner[1], v1, u_inner[3]) +end + +@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, + <:Adiabatic})(flux_inner, + u_inner, + orientation::Integer, + direction, + x, t, + operator_type::Divergence, + equations::CompressibleNavierStokesDiffusion1D{ + GradientVariablesPrimitive + }) + # rho, v1, v2, _ = u_inner + normal_heat_flux = boundary_condition.boundary_condition_heat_flux.boundary_value_normal_flux_function(x, + t, + equations) + v1 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, + equations) + _, tau_1n, _ = flux_inner # extract fluxes for 2nd equation + normal_energy_flux = v1 * tau_1n + normal_heat_flux + return SVector(flux_inner[1], flux_inner[2], normal_energy_flux) +end + +@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, + <:Isothermal})(flux_inner, + u_inner, + orientation::Integer, + direction, + x, t, + operator_type::Gradient, + equations::CompressibleNavierStokesDiffusion1D{ + GradientVariablesPrimitive + }) + v1 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, + equations) + T = boundary_condition.boundary_condition_heat_flux.boundary_value_function(x, t, + equations) + return SVector(u_inner[1], v1, T) +end + +@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, + <:Isothermal})(flux_inner, + u_inner, + orientation::Integer, + direction, + x, t, + operator_type::Divergence, + equations::CompressibleNavierStokesDiffusion1D{ + GradientVariablesPrimitive + }) + return flux_inner +end + +# specialized BC impositions for GradientVariablesEntropy. + +# This should return a SVector containing the boundary values of entropy variables. +# Here, `w_inner` are the transformed variables (e.g., entropy variables). +# +# Taken from "Entropy stable modal discontinuous Galerkin schemes and wall boundary conditions +# for the compressible Navier-Stokes equations" by Chan, Lin, Warburton 2022. +# DOI: 10.1016/j.jcp.2021.110723 +@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, + <:Adiabatic})(flux_inner, + w_inner, + orientation::Integer, + direction, + x, t, + operator_type::Gradient, + equations::CompressibleNavierStokesDiffusion1D{ + GradientVariablesEntropy + }) + v1 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, + equations) + negative_rho_inv_p = w_inner[3] # w_3 = -rho / p + return SVector(w_inner[1], -v1 * negative_rho_inv_p, negative_rho_inv_p) +end + +# this is actually identical to the specialization for GradientVariablesPrimitive, but included for completeness. +@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, + <:Adiabatic})(flux_inner, + w_inner, + orientation::Integer, + direction, + x, t, + operator_type::Divergence, + equations::CompressibleNavierStokesDiffusion1D{ + GradientVariablesEntropy + }) + normal_heat_flux = boundary_condition.boundary_condition_heat_flux.boundary_value_normal_flux_function(x, + t, + equations) + v1 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, + equations) + _, tau_1n, _ = flux_inner # extract fluxes for 2nd equation + normal_energy_flux = v1 * tau_1n + normal_heat_flux + return SVector(flux_inner[1], flux_inner[2], normal_energy_flux) +end + +@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, + <:Isothermal})(flux_inner, + w_inner, + orientation::Integer, + direction, + x, t, + operator_type::Gradient, + equations::CompressibleNavierStokesDiffusion1D{ + GradientVariablesEntropy + }) + v1 = 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 / T = -v1 * w3. + w3 = -1 / T + return SVector(w_inner[1], -v1 * w3, w3) +end + +@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, + <:Isothermal})(flux_inner, + w_inner, + orientation::Integer, + direction, + x, t, + operator_type::Divergence, + equations::CompressibleNavierStokesDiffusion1D{ + GradientVariablesEntropy + }) + return SVector(flux_inner[1], flux_inner[2], flux_inner[3]) +end diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index a1f11717e69..f762fe5d5ee 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -29,7 +29,7 @@ The particular form of the compressible Navier-Stokes implemented is = \nabla \cdot \begin{pmatrix} -0 \\ \underline{\tau} \\ \underline{\tau}\mathbf{v} - \nabla q +0 \\ \underline{\tau} \\ \underline{\tau}\mathbf{v} - \mathbf{q} \end{pmatrix} ``` where the system is closed with the ideal gas assumption giving @@ -44,7 +44,7 @@ are built from the viscous stress tensor ``` where ``\underline{I}`` is the ``2\times 2`` identity matrix and the heat flux is ```math -\nabla q = -\kappa\nabla\left(T\right),\quad T = \frac{p}{R\rho} +\mathbf{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, @@ -55,7 +55,7 @@ the thermal conductivity is 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) +\mathbf{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. @@ -93,24 +93,6 @@ struct CompressibleNavierStokesDiffusion2D{GradientVariables, RealT <: Real, gradient_variables::GradientVariables # GradientVariablesPrimitive or GradientVariablesEntropy end -""" -!!! warning "Experimental code" - This code is experimental and may be changed or removed in any future release. - -`GradientVariablesPrimitive` and `GradientVariablesEntropy` are gradient variable type parameters -for `CompressibleNavierStokesDiffusion2D`. By default, the gradient variables are set to be -`GradientVariablesPrimitive`. Specifying `GradientVariablesEntropy` instead uses the entropy variable -formulation from -- Hughes, Mallet, Franca (1986) - A new finite element formulation for computational fluid dynamics: I. Symmetric forms of the - compressible Euler and Navier-Stokes equations and the second law of thermodynamics. - [https://doi.org/10.1016/0045-7825(86)90127-1](https://doi.org/10.1016/0045-7825(86)90127-1) - -Under `GradientVariablesEntropy`, the Navier-Stokes discretization is provably entropy stable. -""" -struct GradientVariablesPrimitive end -struct GradientVariablesEntropy end - # default to primitive gradient variables function CompressibleNavierStokesDiffusion2D(equations::CompressibleEulerEquations2D; mu, Prandtl, @@ -315,59 +297,6 @@ end return dv2dx - dv1dy end -# TODO: can we generalize this to MHD? -""" - struct BoundaryConditionNavierStokesWall - -Creates a wall-type boundary conditions for the compressible Navier-Stokes equations. -The fields `boundary_condition_velocity` and `boundary_condition_heat_flux` are intended -to be boundary condition types such as the `NoSlip` velocity boundary condition and the -`Adiabatic` or `Isothermal` heat boundary condition. - -!!! warning "Experimental feature" - This is an experimental feature and may change in future releases. -""" -struct BoundaryConditionNavierStokesWall{V, H} - boundary_condition_velocity::V - boundary_condition_heat_flux::H -end - -""" - struct NoSlip - -Use to create a no-slip boundary condition with `BoundaryConditionNavierStokesWall`. The field `boundary_value_function` -should be a function with signature `boundary_value_function(x, t, equations)` -and should return a `SVector{NDIMS}` whose entries are the velocity vector at a -point `x` and time `t`. -""" -struct NoSlip{F} - boundary_value_function::F # value of the velocity vector on the boundary -end - -""" - struct Isothermal - -Used to create a no-slip boundary condition with [`BoundaryConditionNavierStokesWall`](@ref). -The field `boundary_value_function` should be a function with signature -`boundary_value_function(x, t, equations)` and return a scalar value for the -temperature at point `x` and time `t`. -""" -struct Isothermal{F} - boundary_value_function::F # value of the temperature on the boundary -end - -""" - struct Adiabatic - -Used to create a no-slip boundary condition with [`BoundaryConditionNavierStokesWall`](@ref). -The field `boundary_value_normal_flux_function` should be a function with signature -`boundary_value_normal_flux_function(x, t, equations)` and return a scalar value for the -normal heat flux at point `x` and time `t`. -""" -struct Adiabatic{F} - boundary_value_normal_flux_function::F # scaled heat flux 1/T * kappa * dT/dn -end - @inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, <:Adiabatic})(flux_inner, u_inner, diff --git a/src/equations/compressible_navier_stokes_3d.jl b/src/equations/compressible_navier_stokes_3d.jl index 0b770dff1ca..166b53bf615 100644 --- a/src/equations/compressible_navier_stokes_3d.jl +++ b/src/equations/compressible_navier_stokes_3d.jl @@ -29,7 +29,7 @@ The particular form of the compressible Navier-Stokes implemented is = \nabla \cdot \begin{pmatrix} -0 \\ \underline{\tau} \\ \underline{\tau}\mathbf{v} - \nabla q +0 \\ \underline{\tau} \\ \underline{\tau}\mathbf{v} - \mathbf{q} \end{pmatrix} ``` where the system is closed with the ideal gas assumption giving @@ -44,7 +44,7 @@ are built from the viscous stress tensor ``` where ``\underline{I}`` is the ``3\times 3`` identity matrix and the heat flux is ```math -\nabla q = -\kappa\nabla\left(T\right),\quad T = \frac{p}{R\rho} +\mathbf{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, @@ -55,7 +55,7 @@ the thermal conductivity is 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) +\mathbf{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. diff --git a/src/equations/equations_parabolic.jl b/src/equations/equations_parabolic.jl index 6c0be43798a..66214025044 100644 --- a/src/equations/equations_parabolic.jl +++ b/src/equations/equations_parabolic.jl @@ -11,5 +11,7 @@ include("laplace_diffusion_2d.jl") # Compressible Navier-Stokes equations abstract type AbstractCompressibleNavierStokesDiffusion{NDIMS, NVARS} <: AbstractEquationsParabolic{NDIMS, NVARS} end +include("compressible_navier_stokes.jl") +include("compressible_navier_stokes_1d.jl") include("compressible_navier_stokes_2d.jl") include("compressible_navier_stokes_3d.jl") diff --git a/test/test_parabolic_1d.jl b/test/test_parabolic_1d.jl index 1aaf23d576a..06a55100d62 100644 --- a/test/test_parabolic_1d.jl +++ b/test/test_parabolic_1d.jl @@ -19,7 +19,40 @@ isdir(outdir) && rm(outdir, recursive=true) linf = [2.847421658558336e-05] ) end - + + @trixi_testset "TreeMesh1D: elixir_navierstokes_convergence_periodic.jl" begin + @test_trixi_include(joinpath(examples_dir(), "tree_1d_dgsem", "elixir_navierstokes_convergence_periodic.jl"), + l2 = [0.0001133835907077494, 6.226282245610444e-5, 0.0002820171699999139], + linf = [0.0006255102377159538, 0.00036195501456059986, 0.0016147729485886941] + ) + end + + @trixi_testset "TreeMesh1D: elixir_navierstokes_convergence_periodic.jl: GradientVariablesEntropy" begin + @test_trixi_include(joinpath(examples_dir(), "tree_1d_dgsem", "elixir_navierstokes_convergence_periodic.jl"), + equations_parabolic = CompressibleNavierStokesDiffusion1D(equations, mu=mu(), + Prandtl=prandtl_number(), + gradient_variables = GradientVariablesEntropy()), + l2 = [0.00011310615871043463, 6.216495207074201e-5, 0.00028195843110817814], + linf = [0.0006240837363233886, 0.0003616694320713876, 0.0016147339542413874] + ) + end + + @trixi_testset "TreeMesh1D: elixir_navierstokes_convergence_walls.jl" begin + @test_trixi_include(joinpath(examples_dir(), "tree_1d_dgsem", "elixir_navierstokes_convergence_walls.jl"), + l2 = [0.00047023310868269237, 0.00032181736027057234, 0.0014966266486095025], + linf = [0.002996375101363302, 0.002863904256059634, 0.012691132946258676] + ) + end + + @trixi_testset "TreeMesh1D: elixir_navierstokes_convergence_walls.jl: GradientVariablesEntropy" begin + @test_trixi_include(joinpath(examples_dir(), "tree_1d_dgsem", "elixir_navierstokes_convergence_walls.jl"), + equations_parabolic = CompressibleNavierStokesDiffusion1D(equations, mu=mu(), + Prandtl=prandtl_number(), + gradient_variables = GradientVariablesEntropy()), + l2 = [0.0004608500483647771, 0.00032431091222851285, 0.0015159733360626845], + linf = [0.002754803146635787, 0.0028567714697580906, 0.012941794048176192] + ) + end end # Clean up afterwards: delete Trixi output directory