diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl index 03b0f04c7a2..09b8ec05187 100644 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl @@ -131,14 +131,16 @@ lift_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names, l_inf())) drag_coefficient_force = AnalysisSurfaceIntegral(semi, force_boundary_names, - Trixi.DragCoefficientShearStress(aoa(), rho_inf(), - u_inf(equations), - l_inf())) + DragCoefficientShearStress(aoa(), + rho_inf(), + u_inf(equations), + l_inf())) lift_coefficient_force = AnalysisSurfaceIntegral(semi, force_boundary_names, - Trixi.LiftCoefficientShearStress(aoa(), rho_inf(), - u_inf(equations), - l_inf())) + LiftCoefficientShearStress(aoa(), + rho_inf(), + u_inf(equations), + l_inf())) analysis_callback = AnalysisCallback(semi, interval = analysis_interval, output_directory = "out", @@ -156,7 +158,7 @@ save_solution = SaveSolutionCallback(interval = 500, save_final_solution = true, solution_variables = cons2prim) -callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback)#, save_solution) +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution) ############################################################################### # run the simulation diff --git a/src/Trixi.jl b/src/Trixi.jl index 168441513ad..57855320d83 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -262,7 +262,8 @@ export SummaryCallback, SteadyStateCallback, AnalysisCallback, AliveCallback, AMRCallback, StepsizeCallback, GlmSpeedCallback, LBMCollisionCallback, EulerAcousticsCouplingCallback, TrivialCallback, AnalysisCallbackCoupled, - AnalysisSurfaceIntegral, DragCoefficientPressure, LiftCoefficientPressure + AnalysisSurfaceIntegral, DragCoefficientPressure, LiftCoefficientPressure, + DragCoefficientShearStress, LiftCoefficientShearStress export load_mesh, load_time, load_timestep, load_timestep!, load_dt, load_adaptive_time_integrator! diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl index 6d0ee17a0e6..8ac65862e57 100644 --- a/src/callbacks_step/analysis.jl +++ b/src/callbacks_step/analysis.jl @@ -683,9 +683,8 @@ include("analysis_dg3d_parallel.jl") # Note that this needs to be included after `analysis_surface_integral_2d.jl` to # have `VariableViscous` available. function analyze(quantity::Union{typeof(enstrophy), - AnalysisSurfaceIntegral{Semidiscretization, Variable}}, du, - u, - t, + AnalysisSurfaceIntegral{Semidiscretization, Variable}}, + du, u, t, semi::SemidiscretizationHyperbolicParabolic) where {Semidiscretization, Variable <: VariableViscous} diff --git a/src/callbacks_step/analysis_surface_integral_2d.jl b/src/callbacks_step/analysis_surface_integral_2d.jl index e3c6633621e..c40c9bb5741 100644 --- a/src/callbacks_step/analysis_surface_integral_2d.jl +++ b/src/callbacks_step/analysis_surface_integral_2d.jl @@ -62,7 +62,7 @@ struct DragCoefficientPressure{RealT <: Real} end # Abstract base type used for dispatch of `analyze` for quantities -# requiring gradients +# requiring gradients of the velocity field. abstract type VariableViscous end struct LiftCoefficientShearStress{RealT <: Real} <: VariableViscous @@ -189,12 +189,17 @@ function (drag_coefficient::DragCoefficientPressure)(u, normal_direction, equati return p * n / (0.5 * rhoinf * uinf^2 * linf) end -function viscous_stress_vector(u, normal_direction, equations_parabolic, gradients_1, gradients_2) - _, dv1dx, dv2dx, _ = convert_derivative_to_primitive(u, gradients_1, equations_parabolic) - _, dv1dy, dv2dy, _ = convert_derivative_to_primitive(u, gradients_2, equations_parabolic) - - # Normalize normal direction, should point *into* the fluid => *(-1) - n_normal = -normal_direction / norm(normal_direction) +# Compute the three different components of the viscous stress tensor +# (tau_11, tau_12, tau_22) based on the gradients of the velocity field. +# This is required for drag and lift coefficients based on shear stress, +# as well as for the non-integrated quantities such as +# skin friction coefficient (to be added). +function viscous_stress_tensor(u, normal_direction, equations_parabolic, + gradients_1, gradients_2) + _, dv1dx, dv2dx, _ = convert_derivative_to_primitive(u, gradients_1, + equations_parabolic) + _, dv1dy, dv2dy, _ = convert_derivative_to_primitive(u, gradients_2, + equations_parabolic) # Components of viscous stress tensor # (4/3 * (v1)_x - 2/3 * (v2)_y) @@ -205,19 +210,31 @@ function viscous_stress_vector(u, normal_direction, equations_parabolic, gradien # (4/3 * (v2)_y - 2/3 * (v1)_x) tau_22 = 4.0 / 3.0 * dv2dy - 2.0 / 3.0 * dv1dx + mu = dynamic_viscosity(u, equations_parabolic) + + return mu .* (tau_11, tau_12, tau_22) +end + +function viscous_stress_vector(u, normal_direction, equations_parabolic, + gradients_1, gradients_2) + # Normalize normal direction, should point *into* the fluid => *(-1) + n_normal = -normal_direction / norm(normal_direction) + + tau_11, tau_12, tau_22 = viscous_stress_tensor(u, normal_direction, + equations_parabolic, + gradients_1, gradients_2) + # Viscous stress vector: Stress tensor * normal vector visc_stress_vector_1 = tau_11 * n_normal[1] + tau_12 * n_normal[2] visc_stress_vector_2 = tau_12 * n_normal[1] + tau_22 * n_normal[2] - mu = dynamic_viscosity(u, equations_parabolic) - - return mu .* (visc_stress_vector_1, visc_stress_vector_2) + return (visc_stress_vector_1, visc_stress_vector_2) end function (drag_coefficient::LiftCoefficientShearStress)(u, normal_direction, equations_parabolic, gradients_1, gradients_2) - visc_stress_vector = viscous_stress_vector(u, normal_direction, equations_parabolic, + visc_stress_vector = viscous_stress_vector(u, normal_direction, equations_parabolic, gradients_1, gradients_2) @unpack psi, rhoinf, uinf, linf = drag_coefficient.force_state return (visc_stress_vector[1] * psi[1] + visc_stress_vector[2] * psi[2]) / @@ -227,7 +244,7 @@ end function (drag_coefficient::DragCoefficientShearStress)(u, normal_direction, equations_parabolic, gradients_1, gradients_2) - visc_stress_vector = viscous_stress_vector(u, normal_direction, equations_parabolic, + visc_stress_vector = viscous_stress_vector(u, normal_direction, equations_parabolic, gradients_1, gradients_2) @unpack psi, rhoinf, uinf, linf = drag_coefficient.force_state return (visc_stress_vector[1] * psi[1] + visc_stress_vector[2] * psi[2]) / @@ -277,9 +294,9 @@ function analyze(surface_variable::AnalysisSurfaceIntegral, du, u, t, end function analyze(surface_variable::AnalysisSurfaceIntegral{Semidiscretization, - Variable}, du, u, t, - mesh::P4estMesh{2}, - equations, equations_parabolic, + Variable}, + du, u, t, mesh::P4estMesh{2}, + equations, equations_parabolic::AbstractEquationsParabolic, dg::DGSEM, cache, cache_parabolic) where {Semidiscretization, Variable <: VariableViscous}