Skip to content

Commit

Permalink
Some fixes
Browse files Browse the repository at this point in the history
  • Loading branch information
Arpit-Babbar committed Mar 28, 2024
1 parent 7779eb1 commit 16fc47f
Show file tree
Hide file tree
Showing 3 changed files with 30 additions and 30 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -87,11 +87,11 @@ l_inf = 1.0 # Length of airfoil
force_boundary_names = [:AirfoilBottom, :AirfoilTop]
drag_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names,
DragCoefficientPressure(aoa(), rho_inf(),
u_inf(equations), linf))
u_inf(equations), l_inf))

lift_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names,
LiftCoefficientPressure(aoa(), rho_inf(),
u_inf(equations), linf))
u_inf(equations), l_inf))

analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
output_directory = "out",
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@ using Trixi
###############################################################################
# semidiscretization of the compressible Euler equations

# Laminar transonic flow around an airfoil.
# Laminar transonic flow around a NACA0012 airfoil.

# This test is taken from the paper of Swanson and Langer. The values for the drag and lift coefficients
# from Case 5 in Table 3 are used to validate the scheme and computation of surface forces.
Expand Down Expand Up @@ -113,7 +113,7 @@ semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabol
###############################################################################
# ODE solvers

# Run for a long time to reach a steady state
# Run for a long time to reach a state where forces stabilize up to 3 digits
tspan = (0.0, 10.0)
ode = semidiscretize(semi, tspan)

Expand Down
52 changes: 26 additions & 26 deletions src/callbacks_step/analysis_surface_integral_2d.jl
Original file line number Diff line number Diff line change
Expand Up @@ -9,14 +9,14 @@
# surface forces

"""
AnalysisSurfaceIntegral{Semidiscretization, Variable}(semi,
boundary_symbol_or_boundary_symbols,
AnalysisSurfaceIntegral{Semidiscretization, Variable}(semi,
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`
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
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.
"""
Expand Down Expand Up @@ -50,7 +50,7 @@ struct ForceState{RealT <: Real}
psi::Tuple{RealT, RealT} # Unit vector normal or parallel to freestream
rhoinf::RealT
uinf::RealT
linf::RealT
l_inf::RealT
end

struct LiftCoefficientPressure{RealT <: Real}
Expand All @@ -61,70 +61,70 @@ struct DragCoefficientPressure{RealT <: Real}
force_state::ForceState{RealT}
end

"""
LiftCoefficientPressure(aoa, rhoinf, uinf, linf)
"""
LiftCoefficientPressure(aoa, rhoinf, uinf, l_inf)
Compute the lift coefficient
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}}
```
based on the pressure distribution along a boundary.
Supposed to be used in conjunction with [`AnalysisSurfaceIntegral`](@ref)
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)
- `l_inf::Real`: Reference length of geometry (e.g. airfoil chord length)
"""
function LiftCoefficientPressure(aoa, rhoinf, uinf, linf)
function LiftCoefficientPressure(aoa, rhoinf, uinf, l_inf)
# 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
# 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, linf))
return LiftCoefficientPressure(ForceState(psi_lift, rhoinf, uinf, l_inf))
end

"""
DragCoefficientPressure(aoa, rhoinf, uinf, linf)
"""
DragCoefficientPressure(aoa, rhoinf, uinf, l_inf)
Compute the drag coefficient
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}}
```
based on the pressure distribution along a boundary.
Supposed to be used in conjunction with [`AnalysisSurfaceIntegral`](@ref)
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)
- `l_inf::Real`: Reference length of geometry (e.g. airfoil chord length)
"""
function DragCoefficientPressure(aoa, rhoinf, uinf, linf)
function DragCoefficientPressure(aoa, rhoinf, uinf, l_inf)
# `psi_drag` is the unit vector in direction of the freestream.
psi_drag = (cos(aoa), sin(aoa))
return DragCoefficientPressure(ForceState(psi_drag, rhoinf, uinf, linf))
return DragCoefficientPressure(ForceState(psi_drag, rhoinf, uinf, l_inf))
end

function (lift_coefficient::LiftCoefficientPressure)(u, normal_direction, equations)
p = pressure(u, equations)
@unpack psi, rhoinf, uinf, linf = lift_coefficient.force_state
@unpack psi, rhoinf, uinf, l_inf = 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 * linf)
return p * n / (0.5 * rhoinf * uinf^2 * l_inf)
end

function (drag_coefficient::DragCoefficientPressure)(u, normal_direction, equations)
p = pressure(u, equations)
@unpack psi, rhoinf, uinf, linf = drag_coefficient.force_state
@unpack psi, rhoinf, uinf, l_inf = 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 * linf)
return p * n / (0.5 * rhoinf * uinf^2 * l_inf)
end

function analyze(surface_variable::AnalysisSurfaceIntegral, du, u, t,
Expand Down Expand Up @@ -152,7 +152,7 @@ function analyze(surface_variable::AnalysisSurfaceIntegral, du, u, t,
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.
# i.e., *into* the structure.
normal_direction = get_normal_direction(direction, contravariant_vectors,
i_node, j_node,
element)
Expand Down

0 comments on commit 16fc47f

Please sign in to comment.