Skip to content

Commit

Permalink
Modularize stress tensor comp for skin friction
Browse files Browse the repository at this point in the history
  • Loading branch information
@DanielDoehring
DanielDoehring committed Apr 1, 2024
1 parent 4fd1593 commit 422e99a
Show file tree
Hide file tree
Showing 4 changed files with 45 additions and 26 deletions.
16 changes: 9 additions & 7 deletions examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand All @@ -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
Expand Down
3 changes: 2 additions & 1 deletion src/Trixi.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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!
Expand Down
5 changes: 2 additions & 3 deletions src/callbacks_step/analysis.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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}
Expand Down
47 changes: 32 additions & 15 deletions src/callbacks_step/analysis_surface_integral_2d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand All @@ -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]) /
Expand All @@ -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]) /
Expand Down Expand Up @@ -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}
Expand Down

0 comments on commit 422e99a

Please sign in to comment.