From c025873b9bacf26da4f8155f425b3a47bab5f510 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Mon, 8 Apr 2024 07:55:49 +0200 Subject: [PATCH] Boundary integration of viscous forces (#1893) * Viscous force computation * Modularize stress tensor comp for skin friction * typo * try aqua * remove semi from AnalysisSurfaceIntegral * bugfix * two methods for passing aqua tests * rename * test viscous forces * Update src/callbacks_step/analysis_surface_integral_2d.jl Co-authored-by: Andrew Winters * docstring * try to get test to run * Try nicer formatting * revert fmt --------- Co-authored-by: Andrew Winters --- ...xir_navierstokes_NACA0012airfoil_mach08.jl | 18 +- src/Trixi.jl | 3 +- src/callbacks_step/analysis.jl | 30 +- .../analysis_surface_integral_2d.jl | 263 +++++++++++++++--- test/test_parabolic_2d.jl | 29 +- 5 files changed, 296 insertions(+), 47 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl index 5bfa21c5322..1b485913ab2 100644 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl @@ -110,7 +110,7 @@ semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabol # ODE solvers # Run for a long time to reach a state where forces stabilize up to 3 digits -tspan = (0.0, 1.0) +tspan = (0.0, 10.0) ode = semidiscretize(semi, tspan) # Callbacks @@ -130,12 +130,26 @@ lift_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names, u_inf(equations), l_inf())) +drag_coefficient_shear_force = AnalysisSurfaceIntegral(semi, force_boundary_names, + DragCoefficientShearStress(aoa(), + rho_inf(), + u_inf(equations), + l_inf())) + +lift_coefficient_shear_force = AnalysisSurfaceIntegral(semi, force_boundary_names, + LiftCoefficientShearStress(aoa(), + rho_inf(), + u_inf(equations), + l_inf())) + analysis_callback = AnalysisCallback(semi, interval = analysis_interval, output_directory = "out", save_analysis = true, analysis_errors = Symbol[], analysis_integrals = (drag_coefficient, - lift_coefficient)) + lift_coefficient, + drag_coefficient_shear_force, + lift_coefficient_shear_force)) alive_callback = AliveCallback(analysis_interval = analysis_interval) diff --git a/src/Trixi.jl b/src/Trixi.jl index 883f8d66f07..3e7fc5ee519 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 62048853b53..8f89af755a2 100644 --- a/src/callbacks_step/analysis.jl +++ b/src/callbacks_step/analysis.jl @@ -9,11 +9,11 @@ # - analysis_interval part as PeriodicCallback called after a certain amount of simulation time """ AnalysisCallback(semi; interval=0, - save_analysis=false, - output_directory="out", - analysis_filename="analysis.dat", - extra_analysis_errors=Symbol[], - extra_analysis_integrals=()) + save_analysis=false, + output_directory="out", + analysis_filename="analysis.dat", + extra_analysis_errors=Symbol[], + extra_analysis_integrals=()) Analyze a numerical solution every `interval` time steps and print the results to the screen. If `save_analysis`, the results are also saved in @@ -634,9 +634,7 @@ pretty_form_utf(quantity) = get_name(quantity) pretty_form_ascii(quantity) = get_name(quantity) # Special analyze for `SemidiscretizationHyperbolicParabolic` such that -# precomputed gradients are available. For now only implemented for the `enstrophy` -#!!! warning "Experimental code" -# This code is experimental and may be changed or removed in any future release. +# precomputed gradients are available. function analyze(quantity::typeof(enstrophy), du, u, t, semi::SemidiscretizationHyperbolicParabolic) mesh, equations, solver, cache = mesh_equations_solver_cache(semi) @@ -695,3 +693,19 @@ include("analysis_surface_integral_2d.jl") include("analysis_dg2d_parallel.jl") include("analysis_dg3d.jl") include("analysis_dg3d_parallel.jl") + +# Special analyze for `SemidiscretizationHyperbolicParabolic` such that +# precomputed gradients are available. Required for `enstrophy` (see above) and viscous forces. +# Note that this needs to be included after `analysis_surface_integral_2d.jl` to +# have `VariableViscous` available. +function analyze(quantity::AnalysisSurfaceIntegral{Variable}, + du, u, t, + semi::SemidiscretizationHyperbolicParabolic) where { + Variable <: + VariableViscous} + mesh, equations, solver, cache = mesh_equations_solver_cache(semi) + equations_parabolic = semi.equations_parabolic + cache_parabolic = semi.cache_parabolic + analyze(quantity, du, u, t, mesh, equations, equations_parabolic, solver, cache, + cache_parabolic) +end diff --git a/src/callbacks_step/analysis_surface_integral_2d.jl b/src/callbacks_step/analysis_surface_integral_2d.jl index 902fb0b4e85..c5a3d885eb0 100644 --- a/src/callbacks_step/analysis_surface_integral_2d.jl +++ b/src/callbacks_step/analysis_surface_integral_2d.jl @@ -13,15 +13,19 @@ boundary_symbol_or_boundary_symbols, variable) - This struct is used to compute the surface integral of a quantity of interest `variable` alongside - the boundary/boundaries associated with particular name(s) given in `boundary_symbol` - or `boundary_symbols`. - For instance, this can be used to compute the lift [`LiftCoefficientPressure`](@ref) or - drag coefficient [`DragCoefficientPressure`](@ref) of e.g. an airfoil with the boundary - name `:Airfoil` in 2D. +This struct is used to compute the surface integral of a quantity of interest `variable` alongside +the boundary/boundaries associated with particular name(s) given in `boundary_symbol` +or `boundary_symbols`. +For instance, this can be used to compute the lift [`LiftCoefficientPressure`](@ref) or +drag coefficient [`DragCoefficientPressure`](@ref) of e.g. an airfoil with the boundary +name `:Airfoil` in 2D. + +- `semi::Semidiscretization`: Passed in to retrieve boundary condition information +- `boundary_symbol_or_boundary_symbols::Symbol|Vector{Symbol}`: Name(s) of the boundary/boundaries + where the quantity of interest is computed +- `variable::Variable`: Quantity of interest, like lift or drag """ -struct AnalysisSurfaceIntegral{Semidiscretization, Variable} - semi::Semidiscretization # passed in to retrieve boundary condition information +struct AnalysisSurfaceIntegral{Variable} indices::Vector{Int} # Indices in `boundary_condition_indices` where quantity of interest is computed variable::Variable # Quantity of interest, like lift or drag @@ -29,8 +33,7 @@ struct AnalysisSurfaceIntegral{Semidiscretization, Variable} @unpack boundary_symbol_indices = semi.boundary_conditions indices = boundary_symbol_indices[boundary_symbol] - return new{typeof(semi), typeof(variable)}(semi, indices, - variable) + return new{typeof(variable)}(indices, variable) end function AnalysisSurfaceIntegral(semi, boundary_symbols::Vector{Symbol}, variable) @@ -41,8 +44,7 @@ struct AnalysisSurfaceIntegral{Semidiscretization, Variable} end sort!(indices) - return new{typeof(semi), typeof(variable)}(semi, indices, - variable) + return new{typeof(variable)}(indices, variable) end end @@ -50,7 +52,7 @@ struct ForceState{RealT <: Real} psi::Tuple{RealT, RealT} # Unit vector normal or parallel to freestream rhoinf::RealT uinf::RealT - l_inf::RealT + linf::RealT end struct LiftCoefficientPressure{RealT <: Real} @@ -61,13 +63,25 @@ struct DragCoefficientPressure{RealT <: Real} force_state::ForceState{RealT} end +# Abstract base type used for dispatch of `analyze` for quantities +# requiring gradients of the velocity field. +abstract type VariableViscous end + +struct LiftCoefficientShearStress{RealT <: Real} <: VariableViscous + force_state::ForceState{RealT} +end + +struct DragCoefficientShearStress{RealT <: Real} <: VariableViscous + force_state::ForceState{RealT} +end + """ - LiftCoefficientPressure(aoa, rhoinf, uinf, l_inf) + LiftCoefficientPressure(aoa, rhoinf, uinf, linf) Compute the lift coefficient ```math C_{L,p} \\coloneqq \\frac{\\oint_{\\partial \\Omega} p \\boldsymbol n \\cdot \\psi_L \\, \\mathrm{d} S} - {0.5 \\cdot \\rho_{\\infty} \\cdot U_{\\infty}^2 \\cdot L_{\\infty}} + {0.5 \\rho_{\\infty} U_{\\infty}^2 L_{\\infty}} ``` based on the pressure distribution along a boundary. Supposed to be used in conjunction with [`AnalysisSurfaceIntegral`](@ref) @@ -76,25 +90,25 @@ which stores the boundary information and semidiscretization. - `aoa::Real`: Angle of attack in radians (for airfoils etc.) - `rhoinf::Real`: Free-stream density - `uinf::Real`: Free-stream velocity -- `l_inf::Real`: Reference length of geometry (e.g. airfoil chord length) +- `linf::Real`: Reference length of geometry (e.g. airfoil chord length) """ -function LiftCoefficientPressure(aoa, rhoinf, uinf, l_inf) +function LiftCoefficientPressure(aoa, rhoinf, uinf, linf) # psi_lift is the normal unit vector to the freestream direction. # Note: The choice of the normal vector psi_lift = (-sin(aoa), cos(aoa)) # leads to positive lift coefficients for positive angles of attack for airfoils. # One could also use psi_lift = (sin(aoa), -cos(aoa)) which results in the same # value, but with the opposite sign. psi_lift = (-sin(aoa), cos(aoa)) - return LiftCoefficientPressure(ForceState(psi_lift, rhoinf, uinf, l_inf)) + return LiftCoefficientPressure(ForceState(psi_lift, rhoinf, uinf, linf)) end """ - DragCoefficientPressure(aoa, rhoinf, uinf, l_inf) + DragCoefficientPressure(aoa, rhoinf, uinf, linf) Compute the drag coefficient ```math C_{D,p} \\coloneqq \\frac{\\oint_{\\partial \\Omega} p \\boldsymbol n \\cdot \\psi_D \\, \\mathrm{d} S} - {0.5 \\cdot \\rho_{\\infty} \\cdot U_{\\infty}^2 \\cdot L_{\\infty}} + {0.5 \\rho_{\\infty} U_{\\infty}^2 L_{\\infty}} ``` based on the pressure distribution along a boundary. Supposed to be used in conjunction with [`AnalysisSurfaceIntegral`](@ref) @@ -103,28 +117,140 @@ which stores the boundary information and semidiscretization. - `aoa::Real`: Angle of attack in radians (for airfoils etc.) - `rhoinf::Real`: Free-stream density - `uinf::Real`: Free-stream velocity -- `l_inf::Real`: Reference length of geometry (e.g. airfoil chord length) +- `linf::Real`: Reference length of geometry (e.g. airfoil chord length) +""" +function DragCoefficientPressure(aoa, rhoinf, uinf, linf) + # `psi_drag` is the unit vector tangent to the freestream direction + psi_drag = (cos(aoa), sin(aoa)) + return DragCoefficientPressure(ForceState(psi_drag, rhoinf, uinf, linf)) +end + +""" + LiftCoefficientShearStress(aoa, rhoinf, uinf, linf) + +Compute the lift coefficient +```math +C_{L,f} \\coloneqq \\frac{\\oint_{\\partial \\Omega} \\boldsymbol \\tau_w \\cdot \\psi_L \\, \\mathrm{d} S} + {0.5 \\rho_{\\infty} U_{\\infty}^2 L_{\\infty}} +``` +based on the wall shear stress vector ``\\tau_w`` along a boundary. +Supposed to be used in conjunction with [`AnalysisSurfaceIntegral`](@ref) +which stores the boundary information and semidiscretization. + +- `aoa::Real`: Angle of attack in radians (for airfoils etc.) +- `rhoinf::Real`: Free-stream density +- `uinf::Real`: Free-stream velocity +- `linf::Real`: Reference length of geometry (e.g. airfoil chord length) +""" +function LiftCoefficientShearStress(aoa, rhoinf, uinf, linf) + # psi_lift is the normal unit vector to the freestream direction. + # Note: The choice of the normal vector psi_lift = (-sin(aoa), cos(aoa)) + # leads to negative lift coefficients for airfoils. + # One could also use psi_lift = (sin(aoa), -cos(aoa)) which results in the same + # value, but with the opposite sign. + psi_lift = (-sin(aoa), cos(aoa)) + return LiftCoefficientShearStress(ForceState(psi_lift, rhoinf, uinf, linf)) +end + +""" + DragCoefficientShearStress(aoa, rhoinf, uinf, linf) + +Compute the drag coefficient +```math +C_{D,f} \\coloneqq \\frac{\\oint_{\\partial \\Omega} \\boldsymbol \\tau_w \\cdot \\psi_D \\, \\mathrm{d} S} + {0.5 \\rho_{\\infty} U_{\\infty}^2 L_{\\infty}} +``` +based on the wall shear stress vector ``\\tau_w`` along a boundary. +Supposed to be used in conjunction with [`AnalysisSurfaceIntegral`](@ref) +which stores the boundary information and semidiscretization. + +- `aoa::Real`: Angle of attack in radians (for airfoils etc.) +- `rhoinf::Real`: Free-stream density +- `uinf::Real`: Free-stream velocity +- `linf::Real`: Reference length of geometry (e.g. airfoil chord length) """ -function DragCoefficientPressure(aoa, rhoinf, uinf, l_inf) +function DragCoefficientShearStress(aoa, rhoinf, uinf, linf) # `psi_drag` is the unit vector tangent to the freestream direction psi_drag = (cos(aoa), sin(aoa)) - return DragCoefficientPressure(ForceState(psi_drag, rhoinf, uinf, l_inf)) + return DragCoefficientShearStress(ForceState(psi_drag, rhoinf, uinf, linf)) end function (lift_coefficient::LiftCoefficientPressure)(u, normal_direction, equations) p = pressure(u, equations) - @unpack psi, rhoinf, uinf, l_inf = lift_coefficient.force_state + @unpack psi, rhoinf, uinf, linf = lift_coefficient.force_state # Normalize as `normal_direction` is not necessarily a unit vector n = dot(normal_direction, psi) / norm(normal_direction) - return p * n / (0.5 * rhoinf * uinf^2 * l_inf) + return p * n / (0.5 * rhoinf * uinf^2 * linf) end function (drag_coefficient::DragCoefficientPressure)(u, normal_direction, equations) p = pressure(u, equations) - @unpack psi, rhoinf, uinf, l_inf = drag_coefficient.force_state + @unpack psi, rhoinf, uinf, linf = drag_coefficient.force_state # Normalize as `normal_direction` is not necessarily a unit vector n = dot(normal_direction, psi) / norm(normal_direction) - return p * n / (0.5 * rhoinf * uinf^2 * l_inf) + return p * n / (0.5 * rhoinf * uinf^2 * linf) +end + +# Compute the three components of the symmetric 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) + tau_11 = 4.0 / 3.0 * dv1dx - 2.0 / 3.0 * dv2dy + # ((v1)_y + (v2)_x) + # stress tensor is symmetric + tau_12 = dv1dy + dv2dx # = tau_21 + # (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] + + return (visc_stress_vector_1, visc_stress_vector_2) +end + +function (lift_coefficient::LiftCoefficientShearStress)(u, normal_direction, + equations_parabolic, + gradients_1, gradients_2) + visc_stress_vector = viscous_stress_vector(u, normal_direction, equations_parabolic, + gradients_1, gradients_2) + @unpack psi, rhoinf, uinf, linf = lift_coefficient.force_state + return (visc_stress_vector[1] * psi[1] + visc_stress_vector[2] * psi[2]) / + (0.5 * rhoinf * uinf^2 * linf) +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, + 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]) / + (0.5 * rhoinf * uinf^2 * linf) end function analyze(surface_variable::AnalysisSurfaceIntegral, du, u, t, @@ -169,21 +295,88 @@ function analyze(surface_variable::AnalysisSurfaceIntegral, du, u, t, return surface_integral end -function pretty_form_ascii(::AnalysisSurfaceIntegral{<:Any, - <:LiftCoefficientPressure{<:Any}}) +function analyze(surface_variable::AnalysisSurfaceIntegral{Variable}, + du, u, t, mesh::P4estMesh{2}, + equations, equations_parabolic, + dg::DGSEM, cache, + cache_parabolic) where {Variable <: VariableViscous} + @unpack boundaries = cache + @unpack surface_flux_values, node_coordinates, contravariant_vectors = cache.elements + @unpack weights = dg.basis + + @unpack indices, variable = surface_variable + + # Additions for parabolic + @unpack viscous_container = cache_parabolic + @unpack gradients = viscous_container + + gradients_x, gradients_y = gradients + + surface_integral = zero(eltype(u)) + index_range = eachnode(dg) + for boundary in indices + element = boundaries.neighbor_ids[boundary] + node_indices = boundaries.node_indices[boundary] + direction = indices2direction(node_indices) + + i_node_start, i_node_step = index_to_start_step_2d(node_indices[1], index_range) + j_node_start, j_node_step = index_to_start_step_2d(node_indices[2], index_range) + + i_node = i_node_start + j_node = j_node_start + for node_index in index_range + u_node = Trixi.get_node_vars(cache.boundaries.u, equations, dg, node_index, + boundary) + # Extract normal direction at nodes which points from the elements outwards, + # i.e., *into* the structure. + normal_direction = get_normal_direction(direction, contravariant_vectors, + i_node, j_node, + element) + + # L2 norm of normal direction (contravariant_vector) is the surface element + dS = weights[node_index] * norm(normal_direction) + + gradients_1 = get_node_vars(gradients_x, equations_parabolic, dg, i_node, + j_node, element) + gradients_2 = get_node_vars(gradients_y, equations_parabolic, dg, i_node, + j_node, element) + + # Integral over whole boundary surface + surface_integral += variable(u_node, normal_direction, equations_parabolic, + gradients_1, gradients_2) * dS + + i_node += i_node_step + j_node += j_node_step + end + end + return surface_integral +end + +function pretty_form_ascii(::AnalysisSurfaceIntegral{<:LiftCoefficientPressure{<:Any}}) "CL_p" end -function pretty_form_utf(::AnalysisSurfaceIntegral{<:Any, - <:LiftCoefficientPressure{<:Any}}) +function pretty_form_utf(::AnalysisSurfaceIntegral{<:LiftCoefficientPressure{<:Any}}) "CL_p" end -function pretty_form_ascii(::AnalysisSurfaceIntegral{<:Any, - <:DragCoefficientPressure{<:Any}}) +function pretty_form_ascii(::AnalysisSurfaceIntegral{<:DragCoefficientPressure{<:Any}}) "CD_p" end -function pretty_form_utf(::AnalysisSurfaceIntegral{<:Any, - <:DragCoefficientPressure{<:Any}}) +function pretty_form_utf(::AnalysisSurfaceIntegral{<:DragCoefficientPressure{<:Any}}) "CD_p" end + +function pretty_form_ascii(::AnalysisSurfaceIntegral{<:LiftCoefficientShearStress{<:Any}}) + "CL_f" +end +function pretty_form_utf(::AnalysisSurfaceIntegral{<:LiftCoefficientShearStress{<:Any}}) + "CL_f" +end + +function pretty_form_ascii(::AnalysisSurfaceIntegral{<:DragCoefficientShearStress{<:Any}}) + "CD_f" +end +function pretty_form_utf(::AnalysisSurfaceIntegral{<:DragCoefficientShearStress{<:Any}}) + "CD_f" +end end # muladd diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index edbdfe5a84f..a18dfb353b7 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -719,7 +719,10 @@ end 1.199362305026636, 0.9077214424040279, 5.666071182328691], tspan=(0.0, 0.001), - initial_refinement_level=0) + initial_refinement_level=0, + # With the default `maxiters = 1` in coverage tests, + # there would be no time steps after the restart. + coverage_override=(maxiters = 10_000,)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let @@ -728,6 +731,30 @@ end du_ode = similar(u_ode) @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end + + u_ode = copy(sol.u[end]) + du_ode = zero(u_ode) # Just a placeholder in this case + + u = Trixi.wrap_array(u_ode, semi) + du = Trixi.wrap_array(du_ode, semi) + + drag_p = Trixi.analyze(drag_coefficient, du, u, tspan[2], mesh, equations, solver, + semi.cache) + lift_p = Trixi.analyze(lift_coefficient, du, u, tspan[2], mesh, equations, solver, + semi.cache) + + drag_f = Trixi.analyze(drag_coefficient_shear_force, du, u, tspan[2], mesh, + equations, equations_parabolic, solver, + semi.cache, semi.cache_parabolic) + lift_f = Trixi.analyze(lift_coefficient_shear_force, du, u, tspan[2], mesh, + equations, equations_parabolic, solver, + semi.cache, semi.cache_parabolic) + + @test isapprox(drag_p, 0.17963843913309516, atol = 1e-13) + @test isapprox(lift_p, 0.26462588007949367, atol = 1e-13) + + @test isapprox(drag_f, 1.5427441885921553, atol = 1e-13) + @test isapprox(lift_f, 0.005621910087395724, atol = 1e-13) end end