From 38d365fded6775ec21bbc470573fcd53f5f8a39e Mon Sep 17 00:00:00 2001 From: Bart Date: Sat, 18 Jul 2026 12:41:23 +0200 Subject: [PATCH 1/6] Add cp/cf for both solvers --- docs/src/functions.md | 18 ++- docs/src/private_functions.md | 11 +- src/VortexStepMethod.jl | 9 +- src/airfoil_aero/AirfoilAero.jl | 6 +- src/airfoil_aero/airfoil_solvers/common.jl | 61 +++++----- .../airfoil_solvers/neuralfoil_solver.jl | 27 ++++- .../airfoil_solvers/xfoil_solver.jl | 14 ++- src/airfoil_aero/cp_gen.jl | 66 ----------- src/airfoil_aero/section_aero_gen.jl | 80 +++++++++++++ src/body_aerodynamics.jl | 2 +- src/cp_polars.jl | 106 ------------------ src/cp_types.jl | 38 ------- src/panel.jl | 29 +++-- src/section_aero.jl | 96 ++++++++++++++++ src/solver.jl | 47 +------- src/wing_geometry.jl | 78 +++++++------ src/yaml_geometry.jl | 16 +-- test/airfoil_aero/test_airfoil_aero.jl | 91 +++++++-------- 18 files changed, 362 insertions(+), 433 deletions(-) delete mode 100644 src/airfoil_aero/cp_gen.jl create mode 100644 src/airfoil_aero/section_aero_gen.jl delete mode 100644 src/cp_polars.jl delete mode 100644 src/cp_types.jl create mode 100644 src/section_aero.jl diff --git a/docs/src/functions.md b/docs/src/functions.md index aebf3568..3462b56e 100644 --- a/docs/src/functions.md +++ b/docs/src/functions.md @@ -9,16 +9,14 @@ calculate_span calculate_projected_area ``` -## Surface-pressure (Cp) tables -The Cp table types and IO live in the core package; the generation of Cp tables -(from XFoil/NeuralFoil) lives in `AirfoilAero`. +## Surface aero (contour + Cp + cf) tables +The per-section surface aero table type and its npz IO live in the core package; the +generation (from XFoil/NeuralFoil) lives in `AirfoilAero`. ```@docs -CpData -CpPolar -read_cp_data -write_cp_data -cp_distribution -delta_cp +SectionAero +section_surface +read_section_aero +write_section_aero ``` ## Airfoil aerodynamics (AirfoilAero) @@ -47,7 +45,7 @@ neuralfoil_aero generate_aero_matrices generate_polar_from_coordinates generate_polar_from_dat -generate_cp_polar +generate_section_aero ``` ## OBJ mesh conversion (ObjAdapter) diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index fe6963d6..d4bfffc1 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -80,9 +80,8 @@ remove_vector_nans generate_polar_data extract_literature_polar_data parse_literature_column -interpolate_cp_to_refined! -prepare_cp_output! -validate_cp_sections +interpolate_section_aero_to_refined! +validate_section_aero ``` ### Examples @@ -115,9 +114,7 @@ flood_outside trace_level_set largest_linking_gap resample_arc -resample_x smooth_turning! -split_surfaces ``` ### NeuralFoil network @@ -143,7 +140,9 @@ read_dat_coordinates write_dat write_polar_csv write_polar_matrix_csv -fill_slice_nans! +flat_plate_cf +neuralfoil_contour_solution +fill_node_nans! ``` ## OBJ mesh conversion (ObjAdapter) diff --git a/src/VortexStepMethod.jl b/src/VortexStepMethod.jl index 2df50a85..e5bd957e 100644 --- a/src/VortexStepMethod.jl +++ b/src/VortexStepMethod.jl @@ -13,6 +13,7 @@ using NonlinearSolve using SciMLBase import NonlinearSolve: solve, solve! using Interpolations +using NPZ using Parameters using Serialization using Timers @@ -42,9 +43,8 @@ export FAILURE, FEASIBLE, INFEASIBLE, SolverStatus export LOOP, NONLIN, SolverType export load_polar_data -# Surface-pressure (Cp) table types + IO (generation lives in AirfoilAero) -export CpData, CpPolar, cp_distribution, delta_cp -export read_cp_data, write_cp_data +# Per-section surface aero (contour + Cp + cf) + npz IO (generation lives in AirfoilAero) +export SectionAero, section_surface, read_section_aero, write_section_aero export plot_combined_analysis, plot_distribution, plot_geometry, plot_polar_data, plot_polars, plot_section_polars, save_plot, show_plot @@ -414,7 +414,7 @@ end # Include core functionality include("settings.jl") -include("cp_types.jl") +include("section_aero.jl") include("wing_geometry.jl") include("polars.jl") include("yaml_geometry.jl") @@ -423,7 +423,6 @@ include("panel.jl") include("body_aerodynamics.jl") include("wake.jl") include("solver.jl") -include("cp_polars.jl") include("plotting_helpers.jl") diff --git a/src/airfoil_aero/AirfoilAero.jl b/src/airfoil_aero/AirfoilAero.jl index 055b0696..6ac100e6 100644 --- a/src/airfoil_aero/AirfoilAero.jl +++ b/src/airfoil_aero/AirfoilAero.jl @@ -7,7 +7,7 @@ using Interpolations using NPZ using Xfoil using Printf: @sprintf -using ..VortexStepMethod: CpData, CpPolar, interpolate_matrix_nans!, write_aero_matrix +using ..VortexStepMethod: SectionAero, interpolate_matrix_nans!, write_aero_matrix include("kulfan.jl") include("shrink_wrap.jl") @@ -19,7 +19,7 @@ include("airfoil_solvers/neuralfoil_solver.jl") include("airfoil_io.jl") include("polar_gen.jl") include("polar_export.jl") -include("cp_gen.jl") +include("section_aero_gen.jl") export KulfanParameters, KulfanFitMethod, LeastSquaresFit export ShrinkWrap, shrink_wrap @@ -28,7 +28,7 @@ export NeuralFoilModel, NeuralFoilResult, load_neuralfoil_model export neuralfoil_aero, neuralfoil_section export AbstractAirfoilSolver, XFoilSolver, NeuralFoilSolver export SectionSolution, DeformedSection, deform_section, analyze_section, analyze_sweep -export create_2d_polars, generate_aero_matrices, generate_cp_polar, lei_poly_coeffs +export create_2d_polars, generate_aero_matrices, generate_section_aero, lei_poly_coeffs export turn_trailing_edge!, get_lower_upper export read_dat_coordinates, write_dat, write_polar_csv export generate_polar_from_coordinates, generate_polar_from_dat, resolve_airfoil diff --git a/src/airfoil_aero/airfoil_solvers/common.jl b/src/airfoil_aero/airfoil_solvers/common.jl index 0ded169d..8bbde69b 100644 --- a/src/airfoil_aero/airfoil_solvers/common.jl +++ b/src/airfoil_aero/airfoil_solvers/common.jl @@ -11,11 +11,13 @@ abstract type AbstractAirfoilSolver end """ SectionSolution -Result of a single 2D analysis: integrated coefficients plus the surface pressure -as `(x/c, Cp)` per surface (each single-valued in x, leading edge to trailing +Result of a single 2D analysis: integrated coefficients plus the full closed surface +as node arrays `(x, y)` with surface pressure `cp` and skin friction `cf` per node +(one continuous contour, trailing edge → upper → leading edge → lower → trailing edge). `confidence` is `1.0`/`NaN` for XFoil (converged/not) or NeuralFoil's -`analysis_confidence`. The generator resamples the `(x, Cp)` pairs onto the shared -chord slices. +`analysis_confidence`. Non-converged angles carry empty node arrays. `cf` is the +tangential skin-friction coefficient (exact from XFoil `bldump`, approximate from a +flat-plate closure for NeuralFoil). """ struct SectionSolution alpha::Float64 @@ -23,10 +25,10 @@ struct SectionSolution cd::Float64 cm::Float64 confidence::Float64 - x_upper::Vector{Float64} - cp_upper::Vector{Float64} - x_lower::Vector{Float64} - cp_lower::Vector{Float64} + x::Vector{Float64} + y::Vector{Float64} + cp::Vector{Float64} + cf::Vector{Float64} end """ @@ -113,33 +115,6 @@ function deform_section(x, y, delta; crease_frac=0.9, thickness_frac=1.0, return DeformedSection(kulfan, xd, yd) end -""" - split_surfaces(x, cp) -> (x_upper, cp_upper, x_lower, cp_lower) - -Split a Selig-ordered `(x, cp)` node set at the leading edge (`argmin(x)`) into -upper and lower surfaces, each ordered leading edge → trailing edge (x increasing). -""" -function split_surfaces(x, cp) - le = argmin(x) - upper = le:-1:1 - lower = le:length(x) - return x[upper], cp[upper], x[lower], cp[lower] -end - -""" - resample_x(xs, values, targets) -> Vector{Float64} - -Linearly interpolate `values` (given at `xs`) onto `targets`, sorting by `x` and -dropping zero-length segments so the breakpoints are strictly increasing. -""" -function resample_x(xs, values, targets) - p = sortperm(xs) - xs, values = xs[p], values[p] - keep = [true; diff(xs) .> 0] - itp = linear_interpolation(xs[keep], values[keep]; extrapolation_bc=Line()) - return itp.(targets) -end - """ analyze_sweep(solver, def, alpha_range, Re) -> Vector{SectionSolution} @@ -150,3 +125,19 @@ NeuralFoil vectorization). function analyze_sweep(solver::AbstractAirfoilSolver, def::DeformedSection, alpha_range, Re) return [analyze_section(solver, def, a, Re) for a in alpha_range] end + +""" + flat_plate_cf(xc, Re) -> Float64 + +Approximate local skin-friction coefficient at chord fraction `xc` (0..1) for +Reynolds number `Re` (`Re_x = Re·xc`), taking the larger of the two standard +flat-plate correlations: the Blasius laminar solution `cf = 0.664·Re_x^(-1/2)` and +Prandtl's one-seventh-power-law turbulent estimate `cf = 0.027·Re_x^(-1/7)` (see +e.g. White, *Viscous Fluid Flow*; Schlichting & Gersten, *Boundary-Layer Theory*). +A stand-in per-node `cf` for backends that do not expose one (NeuralFoil); XFoil +returns the exact distribution via `bldump`. +""" +function flat_plate_cf(xc, Re) + rex = max(Re * max(xc, 1e-3), 1.0) + return max(0.664 / sqrt(rex), 0.027 / rex^(1 / 7)) +end diff --git a/src/airfoil_aero/airfoil_solvers/neuralfoil_solver.jl b/src/airfoil_aero/airfoil_solvers/neuralfoil_solver.jl index 7eaf4e91..693beb5f 100644 --- a/src/airfoil_aero/airfoil_solvers/neuralfoil_solver.jl +++ b/src/airfoil_aero/airfoil_solvers/neuralfoil_solver.jl @@ -24,15 +24,34 @@ end analyze_sweep(solver::NeuralFoilSolver, def, alpha_range, Re) -> Vector{SectionSolution} Evaluate all angles (radians) in one vectorized NeuralFoil call on the deformed -Kulfan parameters. Cp comes at NeuralFoil's fixed station x/c. +Kulfan parameters, then map the station `Cp` onto the deformed contour nodes +(`def.x`, `def.y`) — upper nodes from `cp_upper`, lower from `cp_lower`. `cf` is a +flat-plate closure ([`flat_plate_cf`](@ref)), NeuralFoil not exposing skin friction. """ function analyze_sweep(solver::NeuralFoilSolver, def::DeformedSection, alpha_range, Re) res = neuralfoil_section(def.kulfan, rad2deg.(collect(alpha_range)), Re; model_size=solver.model_size, weights_dir=solver.weights_dir, n_crit=solver.n_crit, xtr_upper=solver.xtr_upper, xtr_lower=solver.xtr_lower) - return [SectionSolution(alpha_range[i], res.cl[i], res.cd[i], res.cm[i], - res.confidence[i], res.x, res.cp_upper[:, i], - res.x, res.cp_lower[:, i]) for i in eachindex(alpha_range)] + x, y = collect(float.(def.x)), collect(float.(def.y)) + le = argmin(x) + cf = [flat_plate_cf(clamp(xk, 0.0, 1.0), Re) for xk in x] + return [neuralfoil_contour_solution(alpha_range[i], res, i, x, y, le, cf) + for i in eachindex(alpha_range)] +end + +""" + neuralfoil_contour_solution(alpha, res, i, x, y, le, cf) -> SectionSolution + +Assemble a full-contour [`SectionSolution`](@ref) for case `i` of a NeuralFoil sweep +`res`: interpolate `res.cp_upper`/`res.cp_lower` (at `res.x`) onto the contour nodes +`(x, y)` split at the leading edge `le`, carrying the precomputed `cf`. +""" +function neuralfoil_contour_solution(alpha, res, i, x, y, le, cf) + up = linear_interpolation(res.x, res.cp_upper[:, i]; extrapolation_bc=Line()) + lo = linear_interpolation(res.x, res.cp_lower[:, i]; extrapolation_bc=Line()) + cp = [(k <= le ? up : lo)(clamp(x[k], 0.0, 1.0)) for k in eachindex(x)] + return SectionSolution(alpha, res.cl[i], res.cd[i], res.cm[i], + res.confidence[i], x, y, cp, cf) end """ diff --git a/src/airfoil_aero/airfoil_solvers/xfoil_solver.jl b/src/airfoil_aero/airfoil_solvers/xfoil_solver.jl index 8ee23319..28ff516e 100644 --- a/src/airfoil_aero/airfoil_solvers/xfoil_solver.jl +++ b/src/airfoil_aero/airfoil_solvers/xfoil_solver.jl @@ -33,8 +33,10 @@ end Set the deformed coordinates once (repaneling if `solver.repanel`), then solve `alpha_range` (radians) sweeping negative and positive angles outward from zero with a -reinit at each side for convergence. Non-converged angles yield empty Cp and `NaN` -confidence. +reinit at each side for convergence. Each converged angle reads the surface pressure +(`Xfoil.cpdump`) and the boundary layer (`Xfoil.bldump`, giving `cf` and the node +coordinates) at the same panel nodes. Non-converged angles yield empty node arrays and +`NaN` confidence. """ function analyze_sweep(solver::XFoilSolver, def::DeformedSection, alpha_range, Re) Xfoil.set_coordinates(def.x, def.y) @@ -50,9 +52,11 @@ function analyze_sweep(solver::XFoilSolver, def::DeformedSection, alpha_range, R xtrip=solver.xtrip, ncrit=solver.ncrit) reinit = false if converged - xc, cp = Xfoil.cpdump() - xu, cu, xl, cl2 = split_surfaces(xc, cp) - sols[ia] = SectionSolution(alpha_range[ia], cl, cd, cm, 1.0, xu, cu, xl, cl2) + _, cp = Xfoil.cpdump() + _, xb, yb, _, _, _, cf = Xfoil.bldump() + n = min(length(cp), length(cf)) + sols[ia] = SectionSolution(alpha_range[ia], cl, cd, cm, 1.0, + xb[1:n], yb[1:n], cp[1:n], cf[1:n]) else sols[ia] = SectionSolution(alpha_range[ia], NaN, NaN, NaN, NaN, Float64[], Float64[], Float64[], Float64[]) diff --git a/src/airfoil_aero/cp_gen.jl b/src/airfoil_aero/cp_gen.jl deleted file mode 100644 index 7dc7e44c..00000000 --- a/src/airfoil_aero/cp_gen.jl +++ /dev/null @@ -1,66 +0,0 @@ -""" - generate_cp_polar(solver, base; alpha_range, delta_range, n_chord, - reynolds_number, crease_frac=0.9, remove_nan=true) -> CpPolar - -Build a chord-slice `Cp(alpha, delta)` table for a base airfoil (Kulfan -parameters) using the given [`AbstractAirfoilSolver`](@ref). For each `delta` the -section is deformed ([`deform_section`](@ref)), analysed over `alpha_range`, and each -surface's `Cp(x/c)` resampled onto the `n_chord` chord slices. Non-converged points -stay `NaN` and, when `remove_nan`, are filled per slice with -`interpolate_matrix_nans!`. -""" -function generate_cp_polar(solver::AbstractAirfoilSolver, base::KulfanParameters; - alpha_range, delta_range, n_chord, reynolds_number, - crease_frac=0.9, remove_nan=true) - x0, y0 = kulfan_to_coordinates(base) - chord_x = [(i - 0.5) / n_chord for i in 1:n_chord] - n_alpha, n_delta = length(alpha_range), length(delta_range) - cp_upper = fill(NaN, n_chord, n_alpha, n_delta) - cp_lower = fill(NaN, n_chord, n_alpha, n_delta) - - for (jd, delta) in enumerate(delta_range) - def = deform_section(x0, y0, delta; crease_frac) - sols = analyze_sweep(solver, def, alpha_range, reynolds_number) - for (ia, sol) in enumerate(sols) - isempty(sol.cp_upper) && continue - cp_upper[:, ia, jd] = resample_x(sol.x_upper, sol.cp_upper, chord_x) - cp_lower[:, ia, jd] = resample_x(sol.x_lower, sol.cp_lower, chord_x) - end - end - - if remove_nan - fill_slice_nans!(cp_upper) - fill_slice_nans!(cp_lower) - end - data = CpData(n_chord, chord_x, collect(float.(alpha_range)), - collect(float.(delta_range)), cp_upper, cp_lower) - return CpPolar(data) -end - -""" - generate_cp_polar(solver, x::Vector, y::Vector; kwargs...) -> CpPolar - -Convenience: [`shrink_wrap`](@ref) the coordinates into a clean airfoil and fit base -Kulfan parameters ([`LeastSquaresFit`](@ref)) first. -""" -function generate_cp_polar(solver::AbstractAirfoilSolver, x::Vector, y::Vector; kwargs...) - xw, yw = shrink_wrap(x, y, ShrinkWrap()) - return generate_cp_polar(solver, fit_kulfan_parameters(xw, yw, LeastSquaresFit()); - kwargs...) -end - -""" - fill_slice_nans!(grid) - -Fill `NaN`s in each chord slice's `(alpha, delta)` matrix with -`interpolate_matrix_nans!`; all-`NaN` slices are left untouched. -""" -function fill_slice_nans!(grid) - for i in axes(grid, 1) - slice = grid[i, :, :] - all(isnan, slice) && continue - interpolate_matrix_nans!(slice; prn=false) - grid[i, :, :] = slice - end - return nothing -end diff --git a/src/airfoil_aero/section_aero_gen.jl b/src/airfoil_aero/section_aero_gen.jl new file mode 100644 index 00000000..4d86c7a2 --- /dev/null +++ b/src/airfoil_aero/section_aero_gen.jl @@ -0,0 +1,80 @@ +""" + generate_section_aero(solver, base; alpha_range, delta_range, reynolds_number, + crease_frac=0.9, remove_nan=true) -> SectionAero + +Build a native-resolution surface table (closed contour + `Cp` + `cf` per node over the +`(alpha, delta)` grid) for a base airfoil (Kulfan parameters) using the given +[`AbstractAirfoilSolver`](@ref). For each `delta` the section is deformed +([`deform_section`](@ref)) and swept over `alpha_range` (radians); each converged angle +contributes its full-contour `(x, y, cp, cf)`. Non-converged points stay `NaN` and, +when `remove_nan`, are filled per node with `interpolate_matrix_nans!`. +""" +function generate_section_aero(solver::AbstractAirfoilSolver, base::KulfanParameters; + alpha_range, delta_range, reynolds_number, crease_frac=0.9, remove_nan=true) + x0, y0 = kulfan_to_coordinates(base) + n_alpha, n_delta = length(alpha_range), length(delta_range) + sols = Vector{Vector{SectionSolution}}(undef, n_delta) + for (jd, delta) in enumerate(delta_range) + def = deform_section(x0, y0, delta; crease_frac) + sols[jd] = analyze_sweep(solver, def, alpha_range, reynolds_number) + end + + n_node = 0 + for col in sols, s in col + if !isempty(s.cp) + n_node = length(s.cp) + break + end + end + n_node == 0 && throw(ErrorException("generate_section_aero: no converged angle")) + + x = fill(NaN, n_node, n_delta) + y = fill(NaN, n_node, n_delta) + cp = fill(NaN, n_node, n_alpha, n_delta) + cf = fill(NaN, n_node, n_alpha, n_delta) + for jd in 1:n_delta, (ia, s) in enumerate(sols[jd]) + (isempty(s.cp) || length(s.cp) != n_node) && continue + cp[:, ia, jd] .= s.cp + cf[:, ia, jd] .= s.cf + if any(isnan, view(x, :, jd)) + x[:, jd] .= s.x + y[:, jd] .= s.y + end + end + + if remove_nan + for i in 1:n_node + fill_node_nans!(cp, i) + fill_node_nans!(cf, i) + end + end + return SectionAero(collect(float.(alpha_range)), collect(float.(delta_range)), + x, y, cp, cf) +end + +""" + generate_section_aero(solver, x::Vector, y::Vector; kwargs...) -> SectionAero + +Convenience: [`shrink_wrap`](@ref) the coordinates and fit base Kulfan parameters +([`LeastSquaresFit`](@ref)) first. +""" +function generate_section_aero(solver::AbstractAirfoilSolver, x::Vector, y::Vector; + kwargs...) + xw, yw = shrink_wrap(x, y, ShrinkWrap()) + return generate_section_aero(solver, fit_kulfan_parameters(xw, yw, LeastSquaresFit()); + kwargs...) +end + +""" + fill_node_nans!(grid, i) + +Fill `NaN`s in node `i`'s `(alpha, delta)` matrix `grid[i, :, :]` with +`interpolate_matrix_nans!`; an all-`NaN` node is left untouched. +""" +function fill_node_nans!(grid, i) + m = grid[i, :, :] + all(isnan, m) && return grid + interpolate_matrix_nans!(m; prn=false) + grid[i, :, :] = m + return grid +end diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl index b384162c..08c6f0b8 100644 --- a/src/body_aerodynamics.jl +++ b/src/body_aerodynamics.jl @@ -249,7 +249,7 @@ function reinit!(body_aero::BodyAerodynamics{P, W, T}; vec = zeros(MVector{3, T}) for wing in body_aero.wings reinit!(wing) - validate_cp_sections(wing.refined_sections) + validate_section_aero(wing.refined_sections) panel_props = wing.panel_props wing_init_aero = init_aero && !_can_skip_panel_aero_reinit(wing, body_aero.panels, idx) diff --git a/src/cp_polars.jl b/src/cp_polars.jl deleted file mode 100644 index de48bbb4..00000000 --- a/src/cp_polars.jl +++ /dev/null @@ -1,106 +0,0 @@ -""" - CpPolar(data::CpData) -> CpPolar - -Build the per-slice `cp(alpha, delta)` interpolants (linear, `Line()` -extrapolation) from raw [`CpData`](@ref). -""" -function CpPolar(data::CpData) - build(grid) = [linear_interpolation((data.alpha_range, data.delta_range), - grid[i, :, :]; extrapolation_bc=Line()) - for i in 1:data.n_chord] - return CpPolar(data, build(data.cp_upper), build(data.cp_lower)) -end - -""" - CpData(polar::CpPolar) -> CpData - -The raw samples underlying a [`CpPolar`](@ref) (strips the interpolants). -""" -CpData(polar::CpPolar) = polar.data - -""" - read_cp_data(path) -> Union{Nothing, CpData} - -Load a section surface-pressure table from a single `cp.csv`; returns `nothing` if -`path` does not exist. Header: `alpha, delta, up@…, lo@…`. `alpha`/`delta` -are in degrees (converted to radians); `chord_x`/`n_chord` come from the `@` -column headers (upper and lower must share the same slices). Each row is one -`(alpha, delta)` point; the rows must form a full `alpha × delta` grid. -""" -function read_cp_data(path::AbstractString) - isfile(path) || return nothing - lines = readlines(String(path)) - header = strip.(split(lines[1], ',')) - xof(prefix) = [parse(Float64, String(h[(length(prefix) + 1):end])) - for h in header if startswith(h, prefix)] - chord_x = xof("up@") - chord_x == xof("lo@") || - throw(ArgumentError("cp.csv upper/lower chord slices differ: $path")) - n_chord = length(chord_x) - n_chord > 0 || throw(ArgumentError("cp.csv has no up@/lo@ columns: $path")) - - rows = [strip.(split(l, ',')) for l in lines[2:end] if !isempty(strip(l))] - alphas = [deg2rad(parse(Float64, String(r[1]))) for r in rows] - deltas = [deg2rad(parse(Float64, String(r[2]))) for r in rows] - alpha_range = sort(unique(alphas)) - delta_range = sort(unique(deltas)) - length(rows) == length(alpha_range) * length(delta_range) || - throw(ArgumentError("cp.csv rows do not form a full alpha × delta grid: $path")) - - cp_upper = fill(NaN, n_chord, length(alpha_range), length(delta_range)) - cp_lower = fill(NaN, n_chord, length(alpha_range), length(delta_range)) - for (k, r) in enumerate(rows) - length(r) == 2 + 2n_chord || - throw(ArgumentError("cp.csv row $(k) has wrong column count: $path")) - ia = searchsortedfirst(alpha_range, alphas[k]) - jd = searchsortedfirst(delta_range, deltas[k]) - vals = parse.(Float64, String.(r[3:end])) - cp_upper[:, ia, jd] .= vals[1:n_chord] - cp_lower[:, ia, jd] .= vals[(n_chord + 1):2n_chord] - end - return CpData(n_chord, chord_x, alpha_range, delta_range, cp_upper, cp_lower) -end - -""" - write_cp_data(path, data::CpData) -> path - -Write a [`CpData`](@ref) to `path` as a `cp.csv` in the [`read_cp_data`](@ref) -format: header `alpha, delta, up@…, lo@…`, one row per `(alpha, delta)` -with `alpha`/`delta` in degrees. -""" -function write_cp_data(path::AbstractString, data::CpData) - xhdr(prefix) = [string(prefix, round(x; digits=5)) for x in data.chord_x] - header = ["alpha"; "delta"; xhdr("up@"); xhdr("lo@")] - open(String(path), "w") do io - println(io, join(header, ",")) - for (jd, delta) in enumerate(data.delta_range), - (ia, alpha) in enumerate(data.alpha_range) - row = [string(round(rad2deg(alpha); digits=6)), - string(round(rad2deg(delta); digits=6))] - append!(row, string.(data.cp_upper[:, ia, jd])) - append!(row, string.(data.cp_lower[:, ia, jd])) - println(io, join(row, ",")) - end - end - return path -end - -""" - cp_distribution(polar::CpPolar, alpha, delta) -> (cp_upper, cp_lower) - -Upper and lower surface `Cp` at every chord slice for `alpha` and `delta` (radians). -""" -function cp_distribution(polar::CpPolar, alpha, delta) - return ([itp(alpha, delta) for itp in polar.cp_upper_interp], - [itp(alpha, delta) for itp in polar.cp_lower_interp]) -end - -""" - delta_cp(polar::CpPolar, alpha, delta) -> Vector{Float64} - -Chordwise load `ΔCp = Cp_lower - Cp_upper` at every chord slice. -""" -function delta_cp(polar::CpPolar, alpha, delta) - upper, lower = cp_distribution(polar, alpha, delta) - return lower .- upper -end diff --git a/src/cp_types.jl b/src/cp_types.jl deleted file mode 100644 index c690718c..00000000 --- a/src/cp_types.jl +++ /dev/null @@ -1,38 +0,0 @@ -""" - CpData - -Raw chord-slice surface-pressure samples for one 2D section (no interpolants) — the -CSV/section-level counterpart of [`CpPolar`](@ref). Sampled at `n_chord` **vertical -chord slices** `chord_x` (x/c), over angle of attack `alpha_range` and trailing-edge -deflection `delta_range` (both radians). - -Fields: -- `n_chord`, `chord_x`: number and x/c of the chord slices. -- `alpha_range`, `delta_range`: grid axes (radians). -- `cp_upper`, `cp_lower`: samples, `n_chord × n_alpha × n_delta`. -""" -struct CpData - n_chord::Int - chord_x::Vector{Float64} - alpha_range::Vector{Float64} - delta_range::Vector{Float64} - cp_upper::Array{Float64,3} - cp_lower::Array{Float64,3} -end - -""" - CpPolar - -Interpolatable surface-pressure table for a 2D section: a [`CpData`](@ref) plus one -`cp(alpha, delta)` interpolant per chord slice for each surface. `ΔCp = Cp_lower - -Cp_upper` is the chordwise load. Build one with `CpPolar(::CpData)`. - -Fields: -- `data`: the underlying [`CpData`](@ref). -- `cp_upper_interp`, `cp_lower_interp`: one `cp(alpha, delta)` interpolant per slice. -""" -struct CpPolar{I} - data::CpData - cp_upper_interp::Vector{I} - cp_lower_interp::Vector{I} -end diff --git a/src/panel.jl b/src/panel.jl index fc3da45e..7207c1b4 100644 --- a/src/panel.jl +++ b/src/panel.jl @@ -19,7 +19,7 @@ Represents a panel in a vortex step method simulation. All points and vectors ar - cl_interp::CL = nothing: lift interpolation (its type is a struct parameter) - cd_interp::CD = nothing: drag interpolation - cm_interp::CM = nothing: moment interpolation -- cp_polar::CP = nothing: optional surface-pressure table, see [CpPolar](@ref) +- section_aero::SA = nothing: optional surface aero table, see [SectionAero](@ref) - `control_point`::Vector{MVec3}: Panel control point - `bound_point_1`::Vector{MVec3}: First bound point - `bound_point_2`::Vector{MVec3}: Second bound point @@ -35,7 +35,7 @@ Represents a panel in a vortex step method simulation. All points and vectors ar SemiInfiniteFilament() ): Panel filaments, see: [BoundFilament](@ref) """ -@with_kw mutable struct Panel{T, CL, CD, CM, CP} +@with_kw mutable struct Panel{T, CL, CD, CM, SA} TE_point_1::MVector{3, T} = zeros(MVector{3, T}) LE_point_1::MVector{3, T} = zeros(MVector{3, T}) TE_point_2::MVector{3, T} = zeros(MVector{3, T}) @@ -50,7 +50,7 @@ Represents a panel in a vortex step method simulation. All points and vectors ar cl_interp::CL = nothing cd_interp::CD = nothing cm_interp::CM = nothing - cp_polar::CP = nothing + section_aero::SA = nothing aero_center::MVector{3, T} = zeros(MVector{3, T}) control_point::MVector{3, T} = zeros(MVector{3, T}) bound_point_1::MVector{3, T} = zeros(MVector{3, T}) @@ -123,11 +123,11 @@ function init_pos!( end """ - build_interps(section_1, section_2, remove_nan) -> (cl, cd, cm, cp) + build_interps(section_1, section_2, remove_nan) -> (cl, cd, cm, section_aero) Build the averaged aerodynamic interpolations for the panel between two sections. -Returns `(cl_interp, cd_interp, cm_interp, cp_polar)`, each `nothing` for models that -do not use it (INVISCID, POLY). `cl`/`cm` clamp (`Flat`) past the alpha range but +Returns `(cl_interp, cd_interp, cm_interp, section_aero)`, each `nothing` for models +that do not use it (INVISCID, POLY). `cl`/`cm` clamp (`Flat`) past the alpha range but extrapolate linearly (`Line`) over delta; `cd` extrapolates linearly in both. The concrete return types parameterise [`Panel`](@ref) — see [`panel_interp_types`](@ref). """ @@ -179,15 +179,14 @@ function build_interps(section_1::Section, section_2::Section, remove_nan) cm_i = linear_interpolation((alphas, deltas), cm; extrapolation_bc=cm_bc) end end - cp = nothing - if section_1.cp_data !== nothing - cp_1, cp_2 = section_1.cp_data, section_2.cp_data - cp_up = (cp_1.cp_upper .+ cp_2.cp_upper) ./ 2 - cp_low = (cp_1.cp_lower .+ cp_2.cp_lower) ./ 2 - cp = CpPolar(CpData(cp_1.n_chord, cp_1.chord_x, - cp_1.alpha_range, cp_1.delta_range, cp_up, cp_low)) + aero = nothing + if section_1.section_aero !== nothing + a1, a2 = section_1.section_aero, section_2.section_aero + aero = SectionAero(a1.alpha_range, a1.delta_range, + (a1.x .+ a2.x) ./ 2, (a1.y .+ a2.y) ./ 2, + (a1.cp .+ a2.cp) ./ 2, (a1.cf .+ a2.cf) ./ 2) end - return cl_i, cd_i, cm_i, cp + return cl_i, cd_i, cm_i, aero end """ @@ -218,7 +217,7 @@ function init_aero!(panel::Panel, section_1::Section, section_2::Section; elseif !(panel.aero_model in (POLAR_VECTORS, POLAR_MATRICES, INVISCID)) throw(ArgumentError("Unsupported aero model: $(panel.aero_model)")) end - panel.cl_interp, panel.cd_interp, panel.cm_interp, panel.cp_polar = + panel.cl_interp, panel.cd_interp, panel.cm_interp, panel.section_aero = build_interps(section_1, section_2, remove_nan) return nothing end diff --git a/src/section_aero.jl b/src/section_aero.jl new file mode 100644 index 00000000..0cdbf047 --- /dev/null +++ b/src/section_aero.jl @@ -0,0 +1,96 @@ +""" + SectionAero + +Per-2D-section surface aerodynamics at native airfoil resolution: the closed contour +`(x, y)` per trailing-edge deflection, plus surface pressure `cp` and skin friction +`cf` per contour node over an `(alpha, delta)` grid (both radians). Replaces the old +chord-slice Cp table. Use [`section_surface`](@ref) to get the interpolated contour + +`cp`/`cf` at any `(alpha, delta)`. + +Fields: +- `alpha_range`, `delta_range`: grid axes (radians). +- `x`, `y`: contour node coordinates, `n_node × n_delta`. +- `cp`, `cf`: surface pressure and skin friction, `n_node × n_alpha × n_delta`. +- `cp_interp`, `cf_interp`, `x_interp`, `y_interp`: per-node interpolants. +""" +struct SectionAero{A,B,C,D} + alpha_range::Vector{Float64} + delta_range::Vector{Float64} + x::Matrix{Float64} + y::Matrix{Float64} + cp::Array{Float64,3} + cf::Array{Float64,3} + cp_interp::Vector{A} + cf_interp::Vector{B} + x_interp::Vector{C} + y_interp::Vector{D} +end + +""" + SectionAero(alpha_range, delta_range, x, y, cp, cf) -> SectionAero + +Build the per-node interpolants (linear, `Line()` extrapolation) from the raw grid. A +single-element `delta_range` is padded to two constant slices so every interpolant +shares one 2D `(alpha, delta)` form. +""" +function SectionAero(alpha_range, delta_range, x, y, cp, cf) + alpha_range = collect(float.(alpha_range)) + delta_range = collect(float.(delta_range)) + if length(delta_range) == 1 + delta_range = [delta_range[1], delta_range[1] + 1.0] + x = hcat(x, x) + y = hcat(y, y) + cp = cat(cp, cp; dims=3) + cf = cat(cf, cf; dims=3) + end + n_node = size(cp, 1) + cp_interp = [linear_interpolation((alpha_range, delta_range), cp[i, :, :]; + extrapolation_bc=Line()) for i in 1:n_node] + cf_interp = [linear_interpolation((alpha_range, delta_range), cf[i, :, :]; + extrapolation_bc=Line()) for i in 1:n_node] + x_interp = [linear_interpolation(delta_range, x[i, :]; extrapolation_bc=Line()) + for i in 1:n_node] + y_interp = [linear_interpolation(delta_range, y[i, :]; extrapolation_bc=Line()) + for i in 1:n_node] + return SectionAero(alpha_range, delta_range, x, y, cp, cf, + cp_interp, cf_interp, x_interp, y_interp) +end + +""" + section_surface(aero::SectionAero, alpha, delta) -> (x, y, cp, cf) + +Interpolated closed contour `(x, y)` and per-node surface pressure `cp` and skin +friction `cf` at `(alpha, delta)` (radians). +""" +function section_surface(aero::SectionAero, alpha, delta) + x = [itp(delta) for itp in aero.x_interp] + y = [itp(delta) for itp in aero.y_interp] + cp = [itp(alpha, delta) for itp in aero.cp_interp] + cf = [itp(alpha, delta) for itp in aero.cf_interp] + return x, y, cp, cf +end + +""" + write_section_aero(path, aero::SectionAero) -> path + +Write a [`SectionAero`](@ref) to `path` as an `.npz` (arrays `alpha_range`, +`delta_range`, `x`, `y`, `cp`, `cf`). +""" +function write_section_aero(path::AbstractString, aero::SectionAero) + npzwrite(String(path), Dict("alpha_range" => aero.alpha_range, + "delta_range" => aero.delta_range, "x" => aero.x, "y" => aero.y, + "cp" => aero.cp, "cf" => aero.cf)) + return path +end + +""" + read_section_aero(path) -> Union{Nothing, SectionAero} + +Load a [`SectionAero`](@ref) from an `.npz`; `nothing` if `path` does not exist. +""" +function read_section_aero(path::AbstractString) + isfile(path) || return nothing + d = npzread(String(path)) + return SectionAero(d["alpha_range"], d["delta_range"], + d["x"], d["y"], d["cp"], d["cf"]) +end diff --git a/src/solver.jl b/src/solver.jl index 0e9f0915..283b45bf 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -49,9 +49,6 @@ Struct for storing the solution of the [solve!](@ref) function. Must contain all cl_dist::Vector{T} = zeros(T, P) cd_dist::Vector{T} = zeros(T, P) cm_dist::Vector{T} = zeros(T, P) - cp_chord_x::Vector{T} = T[] - cp_upper_dist::Matrix{T} = zeros(T, 0, P) - cp_lower_dist::Matrix{T} = zeros(T, 0, P) lift_dist::Vector{T} = zeros(T, P) drag_dist::Vector{T} = zeros(T, P) panel_moment_dist::Vector{T} = zeros(T, P) @@ -257,41 +254,6 @@ function solve!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics, gamma_dist return calc_forces!(solver, body_aero; reference_point, moment_frac) end -""" - prepare_cp_output!(sol::VSMSolution, panels) - -Size and reset the surface-pressure output of `sol` for the current `panels`. When any -panel carries a `cp_polar`, allocate `n_chord × n_panels` matrices filled with `NaN` -and copy the chord slices; otherwise leave the Cp output empty. -""" -function prepare_cp_output!(sol::VSMSolution{P, U, T}, panels) where {P, U, T} - idx = findfirst(p -> p.cp_polar !== nothing, panels) - if idx === nothing - isempty(sol.cp_chord_x) || (sol.cp_chord_x = T[]) - size(sol.cp_upper_dist, 1) == 0 || - (sol.cp_upper_dist = zeros(T, 0, P); sol.cp_lower_dist = zeros(T, 0, P)) - return nothing - end - n_chord = panels[idx].cp_polar.data.n_chord - if size(sol.cp_upper_dist) != (n_chord, P) - sol.cp_upper_dist = fill(T(NaN), n_chord, P) - sol.cp_lower_dist = fill(T(NaN), n_chord, P) - else - fill!(sol.cp_upper_dist, T(NaN)) - fill!(sol.cp_lower_dist, T(NaN)) - end - sol.cp_chord_x = T.(panels[idx].cp_polar.data.chord_x) - return nothing -end - -""" - delta_cp(sol::VSMSolution) -> Matrix - -Chordwise load `ΔCp = Cp_lower - Cp_upper` per chord slice (rows) and panel (columns) -from the last [`solve!`](@ref). Empty when the wing carried no Cp table. -""" -delta_cp(sol::VSMSolution) = sol.cp_lower_dist .- sol.cp_upper_dist - """ calc_forces!(solver::Solver, body_aero::BodyAerodynamics; reference_point=solver.reference_point, moment_frac=0.1) @@ -322,18 +284,11 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics; density = solver.density aerodynamic_model_type = solver.aerodynamic_model_type - prepare_cp_output!(solver.sol, panels) - # Calculate coefficients for each panel for (i, panel) in enumerate(panels) # zero bytes cl_dist[i] = calculate_cl(panel, alpha_dist[i]) cd_dist[i], cm_dist[i] = calculate_cd_cm(panel, alpha_dist[i]) width_dist[i] = panel.width - if panel.cp_polar !== nothing - up, low = cp_distribution(panel.cp_polar, alpha_dist[i], panel.delta) - solver.sol.cp_upper_dist[:, i] .= up - solver.sol.cp_lower_dist[:, i] .= low - end # Geometric AoA using panel-local axes and prescribed # freestream — scalar ops to avoid allocations @@ -1190,7 +1145,7 @@ function _section_with_eltype(section::Section, ::Type{TD}) where TD MVector{3, TD}(section.TE_point), section.aero_model, section.aero_data, - section.cp_data, + section.section_aero, ) end diff --git a/src/wing_geometry.jl b/src/wing_geometry.jl index 6eb95d7d..c10cf495 100644 --- a/src/wing_geometry.jl +++ b/src/wing_geometry.jl @@ -8,19 +8,20 @@ Represents a wing section with leading edge, trailing edge, and aerodynamic prop - `TE_point::MVector{3, T}`: Trailing edge point coordinates - `aero_model::AeroModel`: [AeroModel](@ref) - `aero_data::AeroData`: See: [AeroData](@ref) -- `cp_data::Union{Nothing, CpData}`: optional surface-pressure table, see [CpData](@ref) +- `section_aero::Union{Nothing, SectionAero}`: optional surface aero table, see + [SectionAero](@ref) """ mutable struct Section{T} LE_point::MVector{3, T} TE_point::MVector{3, T} aero_model::AeroModel aero_data::AeroData - cp_data::Union{Nothing, CpData} + section_aero::Union{Nothing, SectionAero} end Section{T}(; LE_point=zeros(MVector{3, T}), TE_point=zeros(MVector{3, T}), - aero_model=INVISCID, aero_data=nothing, cp_data=nothing) where {T} = - Section{T}(LE_point, TE_point, aero_model, aero_data, cp_data) + aero_model=INVISCID, aero_data=nothing, section_aero=nothing) where {T} = + Section{T}(LE_point, TE_point, aero_model, aero_data, section_aero) Section() = Section{Float64}() @@ -43,9 +44,9 @@ function Section(LE_point, TE_point, aero_model) aero_model, nothing, nothing) end -function Section(LE_point, TE_point, aero_model, aero_data, cp_data=nothing) +function Section(LE_point, TE_point, aero_model, aero_data, section_aero=nothing) return Section{Float64}(MVector{3,Float64}(LE_point), MVector{3,Float64}(TE_point), - aero_model, aero_data, cp_data) + aero_model, aero_data, section_aero) end """ @@ -56,7 +57,7 @@ Function to update a [Section](@ref) in place. @inline _section_sort_key(s::Section) = s.LE_point[2] function reinit!(section::Section, LE_point, TE_point, aero_model=nothing, - aero_data=nothing, cp_data=nothing) + aero_data=nothing, section_aero=nothing) section.LE_point .= LE_point section.TE_point .= TE_point (!isnothing(aero_model)) && (section.aero_model = aero_model) @@ -69,7 +70,7 @@ function reinit!(section::Section, LE_point, TE_point, aero_model=nothing, section.aero_data .= aero_data end end - (!isnothing(cp_data)) && (section.cp_data = cp_data) + (!isnothing(section_aero)) && (section.section_aero = section_aero) nothing end @@ -80,27 +81,27 @@ function reinit!(refined_section::Section{Tr}, section::Section) where {Tr} section.TE_point, section.aero_model, section.aero_data, - section.cp_data, + section.section_aero, ) end """ - validate_cp_sections(sections) + validate_section_aero(sections) -Enforce the all-or-none Cp rule and a uniform chord resolution: either every section -in `sections` carries [`CpData`](@ref) or none does, and all present tables share the -same `chord_x` slices. Throws `ArgumentError` on a violation. +Enforce the all-or-none surface-aero rule and a uniform node resolution: either every +section in `sections` carries [`SectionAero`](@ref) or none does, and all present tables +share the same node count. Throws `ArgumentError` on a violation. """ -function validate_cp_sections(sections) - have = [!isnothing(s.cp_data) for s in sections] +function validate_section_aero(sections) + have = [!isnothing(s.section_aero) for s in sections] all(have) || !any(have) || throw(ArgumentError( - "Cp data must be present on all sections or none " * + "Section aero must be present on all sections or none " * "(found $(count(have))/$(length(have))).")) any(have) || return nothing - reference = sections[findfirst(have)].cp_data.chord_x + reference = size(sections[findfirst(have)].section_aero.cp, 1) for s in sections - s.cp_data.chord_x == reference || - throw(ArgumentError("All sections must share the same Cp chord_x slices.")) + size(s.section_aero.cp, 1) == reference || + throw(ArgumentError("All sections must share the same surface node count.")) end return nothing end @@ -701,7 +702,7 @@ end """ add_section!(wing::Wing, LE_point::PosVector, TE_point::PosVector, - aero_model, aero_data::AeroData=nothing, cp_data=nothing) + aero_model, aero_data::AeroData=nothing, section_aero=nothing) Add a new section to the wing. @@ -711,18 +712,19 @@ Add a new section to the wing. - TE_point::PosVector: [PosVector](@ref) of the point on the side of the trailing edge - `aero_model`::AeroModel: [AeroModel](@ref) - `aero_data`::AeroData: See [AeroData](@ref) -- `cp_data`::Union{Nothing, CpData}: optional surface-pressure table, see [CpData](@ref) +- `section_aero`::Union{Nothing, SectionAero}: optional surface aero table, see + [SectionAero](@ref) """ function add_section!(wing::Wing{P, T}, LE_point, TE_point, aero_model::AeroModel, aero_data::AeroData=nothing, - cp_data::Union{Nothing, CpData}=nothing) where {P, T} + section_aero::Union{Nothing, SectionAero}=nothing) where {P, T} if aero_model == POLAR_VECTORS && wing.remove_nan aero_data = remove_vector_nans(aero_data) elseif aero_model == POLAR_MATRICES && wing.remove_nan interpolate_polar_matrix_nans!(aero_data) end push!(wing.unrefined_sections, Section{T}(MVector{3,T}(LE_point), - MVector{3,T}(TE_point), aero_model, aero_data, cp_data)) + MVector{3,T}(TE_point), aero_model, aero_data, section_aero)) wing.n_unrefined_sections = Int16(length(wing.unrefined_sections)) return nothing end @@ -1125,30 +1127,32 @@ function compute_refined_section_interpolation!(wing::AbstractWing{T}) where {T} wing.refined_section_left_idx[n_sections] = Int16(n_unref - 1) wing.refined_section_weight[n_sections] = zero(T) - interpolate_cp_to_refined!(wing) + interpolate_section_aero_to_refined!(wing) return nothing end """ - interpolate_cp_to_refined!(wing) + interpolate_section_aero_to_refined!(wing) -Set each refined section's [`CpData`](@ref) by spanwise-interpolating the unrefined -sections' Cp tables using the refined→unrefined mapping (`refined_section_left_idx`, -`refined_section_weight`). No-op when the unrefined sections carry no Cp. +Set each refined section's [`SectionAero`](@ref) by spanwise-interpolating the unrefined +sections' surface tables (contour, Cp, cf) using the refined→unrefined mapping +(`refined_section_left_idx`, `refined_section_weight`). No-op when the unrefined +sections carry no surface aero. """ -function interpolate_cp_to_refined!(wing::AbstractWing) +function interpolate_section_aero_to_refined!(wing::AbstractWing) unref = wing.unrefined_sections - all(s -> isnothing(s.cp_data), unref) && return nothing + all(s -> isnothing(s.section_aero), unref) && return nothing for i in eachindex(wing.refined_sections) left = Int(wing.refined_section_left_idx[i]) weight = Float64(wing.refined_section_weight[i]) - cp_l = unref[left].cp_data - cp_r = unref[left + 1].cp_data - (isnothing(cp_l) || isnothing(cp_r)) && continue - cp_up = weight .* cp_l.cp_upper .+ (1 - weight) .* cp_r.cp_upper - cp_low = weight .* cp_l.cp_lower .+ (1 - weight) .* cp_r.cp_lower - wing.refined_sections[i].cp_data = CpData(cp_l.n_chord, cp_l.chord_x, - cp_l.alpha_range, cp_l.delta_range, cp_up, cp_low) + aero_l = unref[left].section_aero + aero_r = unref[left + 1].section_aero + (isnothing(aero_l) || isnothing(aero_r)) && continue + blend(a, b) = weight .* a .+ (1 - weight) .* b + wing.refined_sections[i].section_aero = SectionAero( + aero_l.alpha_range, aero_l.delta_range, + blend(aero_l.x, aero_r.x), blend(aero_l.y, aero_r.y), + blend(aero_l.cp, aero_r.cp), blend(aero_l.cf, aero_r.cf)) end return nothing end diff --git a/src/yaml_geometry.jl b/src/yaml_geometry.jl index 0499b00b..5b3ab5d9 100644 --- a/src/yaml_geometry.jl +++ b/src/yaml_geometry.jl @@ -2,7 +2,7 @@ @with_kw struct WingAirfoilInfo csv_file_path::String dat_file::String = "" - cp_file_path::String = "" + aero_file_path::String = "" cl_file_path::String = "" cd_file_path::String = "" cm_file_path::String = "" @@ -266,7 +266,7 @@ function Wing( info_dict = WingAirfoilInfo( csv_file_path = get(airfoil_dict["info_dict"], "csv_file_path", ""), dat_file = get(airfoil_dict["info_dict"], "dat_file", ""), - cp_file_path = get(airfoil_dict["info_dict"], "cp_file_path", ""), + aero_file_path = get(airfoil_dict["info_dict"], "aero_file_path", ""), cl_file_path = get(airfoil_dict["info_dict"], "cl_file_path", ""), cd_file_path = get(airfoil_dict["info_dict"], "cd_file_path", ""), cm_file_path = get(airfoil_dict["info_dict"], "cm_file_path", "")) @@ -275,14 +275,14 @@ function Wing( # Create CSV file mapping from airfoils airfoil_csv_map = Dict{Int64, String}() - airfoil_cp_map = Dict{Int64, String}() + airfoil_aero_map = Dict{Int64, String}() airfoil_matrix_map = Dict{Int64, NTuple{3, String}}() for airfoil in airfoils if !isempty(airfoil.info_dict.csv_file_path) airfoil_csv_map[airfoil.airfoil_id] = airfoil.info_dict.csv_file_path end - if !isempty(airfoil.info_dict.cp_file_path) - airfoil_cp_map[airfoil.airfoil_id] = airfoil.info_dict.cp_file_path + if !isempty(airfoil.info_dict.aero_file_path) + airfoil_aero_map[airfoil.airfoil_id] = airfoil.info_dict.aero_file_path end if !isempty(airfoil.info_dict.cl_file_path) airfoil_matrix_map[airfoil.airfoil_id] = (airfoil.info_dict.cl_file_path, @@ -326,12 +326,12 @@ function Wing( aero_data, aero_model = load_polar_data(csv_file_path) end - cp_file_path = resolve(get(airfoil_cp_map, section.airfoil_id, "")) - cp_data = isempty(cp_file_path) ? nothing : read_cp_data(cp_file_path) + aero_file_path = resolve(get(airfoil_aero_map, section.airfoil_id, "")) + section_aero = isempty(aero_file_path) ? nothing : read_section_aero(aero_file_path) prn && println("Section airfoil_id $(section.airfoil_id): Using $aero_model model") - add_section!(wing, le_coord, te_coord, aero_model, aero_data, cp_data) + add_section!(wing, le_coord, te_coord, aero_model, aero_data, section_aero) end refine!(wing; sort_sections) diff --git a/test/airfoil_aero/test_airfoil_aero.jl b/test/airfoil_aero/test_airfoil_aero.jl index a1ec0a8a..946f140e 100644 --- a/test/airfoil_aero/test_airfoil_aero.jl +++ b/test/airfoil_aero/test_airfoil_aero.jl @@ -5,8 +5,8 @@ using VortexStepMethod.AirfoilAero: KulfanParameters, LeastSquaresFit, ShrinkWra shrink_wrap, fit_kulfan_parameters, kulfan_to_coordinates, neuralfoil_aero, class_function, bernstein_basis, leading_edge_basis, normalize_airfoil -using VortexStepMethod: CpData, CpPolar, read_cp_data, write_cp_data, - cp_distribution, delta_cp +using VortexStepMethod: SectionAero, section_surface, read_section_aero, + write_section_aero read_dat_coords(path) = begin x = Float64[]; y = Float64[] @@ -116,60 +116,55 @@ end @test size(cm) == (length(alpha_range), length(delta_range)) end -@testset "Cp table round-trip and interpolation" begin - n_chord = 5 - chord_x = [0.1, 0.3, 0.5, 0.7, 0.9] +@testset "SectionAero round-trip and interpolation" begin alpha_range = deg2rad.([-5.0, 0.0, 5.0, 10.0]) delta_range = deg2rad.([-3.0, 0.0, 3.0]) + xc = [1.0, 0.5, 0.0, 0.5, 1.0] + yc = [0.0, 0.06, 0.0, -0.04, 0.0] + n_node = length(xc) + x = repeat(xc, 1, length(delta_range)) + y = repeat(yc, 1, length(delta_range)) val(i, ia, jd) = 100i + 10ia + jd - cp_upper = [float(val(i, ia, jd)) for i in 1:n_chord, - ia in eachindex(alpha_range), jd in eachindex(delta_range)] - cp_lower = -cp_upper - data = CpData(n_chord, chord_x, alpha_range, delta_range, cp_upper, cp_lower) - - path = joinpath(mktempdir(), "cp.csv") - write_cp_data(path, data) - back = read_cp_data(path) - @test back.n_chord == n_chord - @test back.chord_x ≈ chord_x - @test back.alpha_range ≈ alpha_range - @test back.delta_range ≈ delta_range - @test back.cp_upper ≈ cp_upper - @test back.cp_lower ≈ cp_lower - @test read_cp_data(joinpath(@__DIR__, "does_not_exist.csv")) === nothing - - polar = CpPolar(data) - @test CpData(polar) === data - up, lo = cp_distribution(polar, alpha_range[2], delta_range[1]) - @test up ≈ cp_upper[:, 2, 1] - @test lo ≈ cp_lower[:, 2, 1] - @test delta_cp(polar, alpha_range[3], delta_range[2]) ≈ - cp_lower[:, 3, 2] .- cp_upper[:, 3, 2] + cp = [float(val(i, ia, jd)) for i in 1:n_node, + ia in eachindex(alpha_range), jd in eachindex(delta_range)] + cf = cp ./ 1000 + aero = SectionAero(alpha_range, delta_range, x, y, cp, cf) + + path = joinpath(mktempdir(), "aero.npz") + write_section_aero(path, aero) + back = read_section_aero(path) + xb, yb, cpb, cfb = section_surface(back, alpha_range[2], delta_range[1]) + @test xb ≈ xc + @test yb ≈ yc + @test cpb ≈ cp[:, 2, 1] + @test cfb ≈ cf[:, 2, 1] + @test read_section_aero(joinpath(@__DIR__, "does_not_exist.npz")) === nothing + + _, _, cpa, _ = section_surface(aero, alpha_range[3], delta_range[2]) + @test cpa ≈ cp[:, 3, 2] end -@testset "generate_cp_polar builds a Cp table" begin +@testset "generate_section_aero builds a surface table" begin truth = KulfanParameters(fill(0.15, 8), fill(-0.15, 8), 0.1, 0.0) alpha_range = deg2rad.(-4.0:2.0:4.0) delta_range = deg2rad.([0.0, 5.0]) - polar = generate_cp_polar(NeuralFoilSolver(model_size="medium"), truth; - alpha_range, delta_range, n_chord=6, reynolds_number=5e5) - @test polar isa CpPolar - @test CpData(polar).n_chord == 6 - up, lo = cp_distribution(polar, alpha_range[2], delta_range[1]) - @test length(up) == 6 - @test all(isfinite, up) && all(isfinite, lo) - @test all(isfinite, delta_cp(polar, alpha_range[end], delta_range[end])) - - up_lo_a, _ = cp_distribution(polar, alpha_range[1], delta_range[1]) - up_hi_a, _ = cp_distribution(polar, alpha_range[end], delta_range[1]) - @test maximum(abs.(up_lo_a .- up_hi_a)) > 0.1 - up_hi_d, _ = cp_distribution(polar, alpha_range[2], delta_range[end]) - @test maximum(abs.(up .- up_hi_d)) > 0.005 - - x, y = read_dat_coords(joinpath(@__DIR__, "data", "test_airfoil.dat")) - polar2 = generate_cp_polar(NeuralFoilSolver(model_size="medium"), x, y; - alpha_range, delta_range=deg2rad.([0.0, 5.0]), n_chord=4, reynolds_number=5e5) - @test CpData(polar2).n_chord == 4 + aero = generate_section_aero(NeuralFoilSolver(model_size="medium"), truth; + alpha_range, delta_range, reynolds_number=5e5) + @test aero isa SectionAero + x, y, cp, cf = section_surface(aero, alpha_range[2], delta_range[1]) + @test length(cp) == size(aero.cp, 1) + @test all(isfinite, cp) && all(isfinite, cf) + @test all(cf .>= 0) + + _, _, cp_lo_a, _ = section_surface(aero, alpha_range[1], delta_range[1]) + _, _, cp_hi_a, _ = section_surface(aero, alpha_range[end], delta_range[1]) + @test maximum(abs.(cp_lo_a .- cp_hi_a)) > 0.1 + + xdat, ydat = read_dat_coords(joinpath(@__DIR__, "data", "test_airfoil.dat")) + aero2 = generate_section_aero(NeuralFoilSolver(model_size="medium"), xdat, ydat; + alpha_range, delta_range=deg2rad.([0.0, 5.0]), reynolds_number=5e5) + @test aero2 isa SectionAero + @test size(aero2.cp, 1) > 0 end @testset "NeuralFoil physical invariants" begin From 19fa8e31eba8e08c91b941fab03cf2a4f661f466 Mon Sep 17 00:00:00 2001 From: Bart Date: Sat, 18 Jul 2026 18:41:48 +0200 Subject: [PATCH 2/6] Use human readable files and only generate/write in submodule --- src/VortexStepMethod.jl | 6 +- src/airfoil_aero/AirfoilAero.jl | 3 +- src/airfoil_aero/polar_export.jl | 21 +++++++ src/airfoil_aero/section_aero_gen.jl | 38 +++++++++++++ src/obj_adapter/obj_to_yaml.jl | 12 +++- src/polars.jl | 36 ------------ src/section_aero.jl | 76 +++++++++++++++++++++----- src/yaml_geometry.jl | 18 +++--- test/airfoil_aero/test_airfoil_aero.jl | 18 +++--- 9 files changed, 154 insertions(+), 74 deletions(-) diff --git a/src/VortexStepMethod.jl b/src/VortexStepMethod.jl index e5bd957e..0d6cdde4 100644 --- a/src/VortexStepMethod.jl +++ b/src/VortexStepMethod.jl @@ -13,7 +13,6 @@ using NonlinearSolve using SciMLBase import NonlinearSolve: solve, solve! using Interpolations -using NPZ using Parameters using Serialization using Timers @@ -43,8 +42,9 @@ export FAILURE, FEASIBLE, INFEASIBLE, SolverStatus export LOOP, NONLIN, SolverType export load_polar_data -# Per-section surface aero (contour + Cp + cf) + npz IO (generation lives in AirfoilAero) -export SectionAero, section_surface, read_section_aero, write_section_aero +# Per-section surface aero (contour + Cp + cf). Human-readable CSV/dat IO: reading here, +# writing in AirfoilAero. +export SectionAero, section_surface, read_section_aero export plot_combined_analysis, plot_distribution, plot_geometry, plot_polar_data, plot_polars, plot_section_polars, save_plot, show_plot diff --git a/src/airfoil_aero/AirfoilAero.jl b/src/airfoil_aero/AirfoilAero.jl index 6ac100e6..25f63899 100644 --- a/src/airfoil_aero/AirfoilAero.jl +++ b/src/airfoil_aero/AirfoilAero.jl @@ -7,7 +7,7 @@ using Interpolations using NPZ using Xfoil using Printf: @sprintf -using ..VortexStepMethod: SectionAero, interpolate_matrix_nans!, write_aero_matrix +using ..VortexStepMethod: SectionAero, interpolate_matrix_nans! include("kulfan.jl") include("shrink_wrap.jl") @@ -29,6 +29,7 @@ export neuralfoil_aero, neuralfoil_section export AbstractAirfoilSolver, XFoilSolver, NeuralFoilSolver export SectionSolution, DeformedSection, deform_section, analyze_section, analyze_sweep export create_2d_polars, generate_aero_matrices, generate_section_aero, lei_poly_coeffs +export write_section_aero export turn_trailing_edge!, get_lower_upper export read_dat_coordinates, write_dat, write_polar_csv export generate_polar_from_coordinates, generate_polar_from_dat, resolve_airfoil diff --git a/src/airfoil_aero/polar_export.jl b/src/airfoil_aero/polar_export.jl index 82b13e26..40b15b02 100644 --- a/src/airfoil_aero/polar_export.jl +++ b/src/airfoil_aero/polar_export.jl @@ -1,3 +1,24 @@ +""" + write_aero_matrix(filepath, matrix, alpha_range, delta_range, label) -> filepath + +Write an `(alpha × delta)` coefficient matrix to a labelled CSV: the header row holds +the flap deflections (`δ=…°`), the first column the angles of attack (`α=…°`), both in +degrees. `alpha_range`/`delta_range` are radians; `label` (e.g. `"C_l"`) names the +coefficient. The submodule-side writer; `read_aero_matrix` (main package) reads it back. +""" +function write_aero_matrix(filepath::AbstractString, matrix::Matrix{Float64}, + alpha_range::Vector{Float64}, delta_range::Vector{Float64}, label::AbstractString) + open(String(filepath), "w") do io + deltas_str = join(("δ=$(round(rad2deg(δ), digits=1))°" for δ in delta_range), ",") + println(io, string(label, "/delta,", deltas_str)) + for i in eachindex(alpha_range) + coeffs_str = join((round(v; digits=4) for v in matrix[i, :]), ",") + println(io, string("α=$(round(rad2deg(alpha_range[i]), digits=1))°,", coeffs_str)) + end + end + return filepath +end + """ generate_aero_matrices(solver, x, y; alpha_range, delta_range, Re, crease_frac=0.75, remove_nan=true, on_deform=nothing) diff --git a/src/airfoil_aero/section_aero_gen.jl b/src/airfoil_aero/section_aero_gen.jl index 4d86c7a2..a1a09be5 100644 --- a/src/airfoil_aero/section_aero_gen.jl +++ b/src/airfoil_aero/section_aero_gen.jl @@ -78,3 +78,41 @@ function fill_node_nans!(grid, i) grid[i, :, :] = m return grid end + +""" + write_node_table(path, aero, values) -> path + +Write a per-node aero table (`Cp` or `cf`, shaped `n_node × n_alpha × n_delta`) as a +human-readable CSV: header `alpha, delta, n0, n1, …` (node columns in the contour node +order), one row per `(alpha, delta)` with angles in degrees. +""" +function write_node_table(path::AbstractString, aero::SectionAero, values) + open(String(path), "w") do io + println(io, "alpha,delta," * join(("n$(k - 1)" for k in 1:size(values, 1)), ",")) + for jd in eachindex(aero.delta_range), ia in eachindex(aero.alpha_range) + row = [rad2deg(aero.alpha_range[ia]), rad2deg(aero.delta_range[jd])] + append!(row, values[:, ia, jd]) + println(io, join(row, ",")) + end + end + return path +end + +""" + write_section_aero(prefix, aero::SectionAero) -> (dat, cp_csv, cf_csv) + +Write a [`SectionAero`](@ref) as human-readable files sharing `prefix`: `{prefix}.dat` +(contour at `delta=0`, `{prefix}_d{deg}.dat` per non-zero deflection) plus +`{prefix}_cp.csv` and `{prefix}_cf.csv` (per-node tables in the contour node order). +`read_section_aero` reads them back. The single writer for surface aero (submodule +side); loading lives in the main package. +""" +function write_section_aero(prefix::AbstractString, aero::SectionAero) + for (jd, d) in enumerate(aero.delta_range) + path = iszero(d) ? "$prefix.dat" : "$(prefix)_d$(round(Int, rad2deg(d))).dat" + write_dat(path, "section", aero.x[:, jd], aero.y[:, jd]) + end + write_node_table("$(prefix)_cp.csv", aero, aero.cp) + write_node_table("$(prefix)_cf.csv", aero, aero.cf) + return "$prefix.dat", "$(prefix)_cp.csv", "$(prefix)_cf.csv" +end diff --git a/src/obj_adapter/obj_to_yaml.jl b/src/obj_adapter/obj_to_yaml.jl index 6ed318b7..33afcb5c 100644 --- a/src/obj_adapter/obj_to_yaml.jl +++ b/src/obj_adapter/obj_to_yaml.jl @@ -127,7 +127,6 @@ function obj_to_yaml(obj_path::String, output_dir::String; ok = Int[] for j in unique(ids) s = stations[j] - dat_rel = joinpath("airfoils", "$j.dat") raw_rel = joinpath("airfoils", "$(j)_raw.dat") csv_rel = joinpath("polars", "$j.csv") try @@ -137,11 +136,16 @@ function obj_to_yaml(obj_path::String, output_dir::String; dat_prefix=joinpath(output_dir, "airfoils", "$j")) clvals = res isa AbstractVector ? collect(sol.cl for sol in res) : vec(res[1]) all(isnan, clvals) && error("solver produced no converged points") - write_dat(joinpath(output_dir, dat_rel), "section_$j", s.x_fit, s.y_fit) + aero = generate_section_aero(aero_solver, s.x_fit, s.y_fit; + alpha_range=deg2rad.(alpha_range), + delta_range=(isnothing(delta_range) ? [0.0] : deg2rad.(collect(delta_range))), + reynolds_number=Float64(Re), crease_frac) + paths = write_section_aero(joinpath(output_dir, "airfoils", "$j"), aero) + dat_rel, cp_rel, cf_rel = (relpath(p, output_dir) for p in paths) write_dat(joinpath(output_dir, raw_rel), "section_$(j)_raw", s.xa, s.ya) push!(airfoil_rows, Any[j, "polar_vectors", Dict("dat_file" => dat_rel, "raw_dat_file" => raw_rel, - "csv_file_path" => csv_rel)]) + "csv_file_path" => csv_rel, "cp_file" => cp_rel, "cf_file" => cf_rel)]) push!(ok, j) finite_cl = [v for v in clvals if !isnan(v)] cl_max = isempty(finite_cl) ? NaN : maximum(finite_cl) @@ -165,6 +169,8 @@ function obj_to_yaml(obj_path::String, output_dir::String; s.TE_point[1], s.TE_point[2], s.TE_point[3]]) end + sort!(section_rows; by = row -> row[3]) # clean spanwise order (by LE_y) + sort!(airfoil_rows; by = row -> row[1]) # airfoils by id yaml_path = joinpath(output_dir, "geometry.yaml") write_geometry_yaml(yaml_path, section_rows, airfoil_rows) verbose && @info "Wrote geometry to $yaml_path ($(length(section_rows)) sections)" diff --git a/src/polars.jl b/src/polars.jl index 8c97ca9f..c04c331a 100644 --- a/src/polars.jl +++ b/src/polars.jl @@ -52,42 +52,6 @@ function interpolate_matrix_nans!(matrix::Matrix{Float64}; prn=true) end -""" - write_aero_matrix(filepath::String, matrix::Matrix{Float64}, - alpha_range::Vector{Float64}, delta_range::Vector{Float64}; - label::String="C_l") - -Write an aerodynamic coefficient matrix to CSV with angle labels. -The first row contains flap deflection angles, first column contains angles of attack. - -# Arguments -- `filepath`: Path to output CSV file -- `matrix`: Matrix of aerodynamic coefficients -- `alpha_range`: Vector of angle of attack values in radians -- `delta_range`: Vector of flap deflection angles in radians -- `label`: Coefficient label for the header -""" -function write_aero_matrix(filepath::AbstractString, matrix::Matrix{Float64}, - alpha_range::Vector{Float64}, delta_range::Vector{Float64}, - label::AbstractString) - open(String(filepath), "w") do io - # Write header with delta values - delta_labels = ["δ=$(round(rad2deg(δ), digits=1))°" for δ in delta_range] - deltas_str = join(delta_labels, ",") - deltas_str isa String || throw(ArgumentError("Failed to serialize delta header labels.")) - header = string(label, "/delta,", deltas_str) - println(io, header) - - # Write data rows with alpha values and coefficients (4 decimals is plenty) - for i in eachindex(alpha_range) - coeffs_str = join((round(v; digits=4) for v in matrix[i, :]), ",") - coeffs_str isa String || throw(ArgumentError("Failed to serialize coefficient row.")) - row = string("α=$(round(rad2deg(alpha_range[i]), digits=1))°,", coeffs_str) - println(io, row) - end - end -end - """ read_aero_matrix(filepath::AbstractString) -> (Matrix{Float64}, Vector{Float64}, Vector{Float64}) diff --git a/src/section_aero.jl b/src/section_aero.jl index 0cdbf047..328d09d7 100644 --- a/src/section_aero.jl +++ b/src/section_aero.jl @@ -71,26 +71,72 @@ function section_surface(aero::SectionAero, alpha, delta) end """ - write_section_aero(path, aero::SectionAero) -> path + read_dat(path) -> (x, y) -Write a [`SectionAero`](@ref) to `path` as an `.npz` (arrays `alpha_range`, -`delta_range`, `x`, `y`, `cp`, `cf`). +Read Selig `.dat` airfoil coordinates (two whitespace-separated columns), skipping the +name/header line and any non-numeric lines. """ -function write_section_aero(path::AbstractString, aero::SectionAero) - npzwrite(String(path), Dict("alpha_range" => aero.alpha_range, - "delta_range" => aero.delta_range, "x" => aero.x, "y" => aero.y, - "cp" => aero.cp, "cf" => aero.cf)) - return path +function read_dat(path::AbstractString) + x = Float64[] + y = Float64[] + for ln in eachline(String(path)) + p = split(strip(ln)) + length(p) >= 2 || continue + xv, yv = tryparse(Float64, p[1]), tryparse(Float64, p[2]) + (xv === nothing || yv === nothing) && continue + push!(x, xv) + push!(y, yv) + end + return x, y +end + +""" + read_node_table(path) -> (alpha, delta, values) + +Read a per-node aero CSV (header `alpha, delta, n0, n1, …`; angles in degrees) into +radian `alpha`/`delta` vectors (one entry per row) and a `nrow × n_node` value matrix. +""" +function read_node_table(path::AbstractString) + lines = [l for l in readlines(String(path)) if !isempty(strip(l))] + rows = [parse.(Float64, split(l, ',')) for l in lines[2:end]] + alpha = deg2rad.([r[1] for r in rows]) + delta = deg2rad.([r[2] for r in rows]) + values = reduce(vcat, (permutedims(r[3:end]) for r in rows)) + return alpha, delta, values end """ - read_section_aero(path) -> Union{Nothing, SectionAero} + read_section_aero(dat_file, cp_file, cf_file) -> Union{Nothing, SectionAero} -Load a [`SectionAero`](@ref) from an `.npz`; `nothing` if `path` does not exist. +Assemble a [`SectionAero`](@ref) from the human-readable files: the airfoil contour +(`dat_file`, plus `{stem}_d{deg}.dat` per non-zero deflection) and the per-node `Cp` +and `cf` tables (`alpha, delta, n0…` in the `.dat` node order). Returns `nothing` if any +file is missing. This is the single loader for both provided and generated aero. """ -function read_section_aero(path::AbstractString) - isfile(path) || return nothing - d = npzread(String(path)) - return SectionAero(d["alpha_range"], d["delta_range"], - d["x"], d["y"], d["cp"], d["cf"]) +function read_section_aero(dat_file::AbstractString, cp_file::AbstractString, + cf_file::AbstractString) + (isfile(dat_file) && isfile(cp_file) && isfile(cf_file)) || return nothing + alpha, delta, cp_rows = read_node_table(cp_file) + _, _, cf_rows = read_node_table(cf_file) + alpha_range = sort(unique(alpha)) + delta_range = sort(unique(delta)) + n_node = size(cp_rows, 2) + stem = replace(String(dat_file), r"\.dat$" => "") + x = fill(NaN, n_node, length(delta_range)) + y = fill(NaN, n_node, length(delta_range)) + for (jd, d) in enumerate(delta_range) + xd, yd = read_dat(iszero(d) ? String(dat_file) : + "$(stem)_d$(round(Int, rad2deg(d))).dat") + x[:, jd] .= xd + y[:, jd] .= yd + end + cp = fill(NaN, n_node, length(alpha_range), length(delta_range)) + cf = fill(NaN, n_node, length(alpha_range), length(delta_range)) + for r in eachindex(alpha) + ia = searchsortedfirst(alpha_range, alpha[r]) + jd = searchsortedfirst(delta_range, delta[r]) + cp[:, ia, jd] .= cp_rows[r, :] + cf[:, ia, jd] .= cf_rows[r, :] + end + return SectionAero(alpha_range, delta_range, x, y, cp, cf) end diff --git a/src/yaml_geometry.jl b/src/yaml_geometry.jl index 5b3ab5d9..4cda79d2 100644 --- a/src/yaml_geometry.jl +++ b/src/yaml_geometry.jl @@ -2,7 +2,8 @@ @with_kw struct WingAirfoilInfo csv_file_path::String dat_file::String = "" - aero_file_path::String = "" + cp_file::String = "" + cf_file::String = "" cl_file_path::String = "" cd_file_path::String = "" cm_file_path::String = "" @@ -266,7 +267,8 @@ function Wing( info_dict = WingAirfoilInfo( csv_file_path = get(airfoil_dict["info_dict"], "csv_file_path", ""), dat_file = get(airfoil_dict["info_dict"], "dat_file", ""), - aero_file_path = get(airfoil_dict["info_dict"], "aero_file_path", ""), + cp_file = get(airfoil_dict["info_dict"], "cp_file", ""), + cf_file = get(airfoil_dict["info_dict"], "cf_file", ""), cl_file_path = get(airfoil_dict["info_dict"], "cl_file_path", ""), cd_file_path = get(airfoil_dict["info_dict"], "cd_file_path", ""), cm_file_path = get(airfoil_dict["info_dict"], "cm_file_path", "")) @@ -275,14 +277,15 @@ function Wing( # Create CSV file mapping from airfoils airfoil_csv_map = Dict{Int64, String}() - airfoil_aero_map = Dict{Int64, String}() + airfoil_surface_map = Dict{Int64, NTuple{3, String}}() airfoil_matrix_map = Dict{Int64, NTuple{3, String}}() for airfoil in airfoils if !isempty(airfoil.info_dict.csv_file_path) airfoil_csv_map[airfoil.airfoil_id] = airfoil.info_dict.csv_file_path end - if !isempty(airfoil.info_dict.aero_file_path) - airfoil_aero_map[airfoil.airfoil_id] = airfoil.info_dict.aero_file_path + if !isempty(airfoil.info_dict.cp_file) && !isempty(airfoil.info_dict.cf_file) + airfoil_surface_map[airfoil.airfoil_id] = (airfoil.info_dict.dat_file, + airfoil.info_dict.cp_file, airfoil.info_dict.cf_file) end if !isempty(airfoil.info_dict.cl_file_path) airfoil_matrix_map[airfoil.airfoil_id] = (airfoil.info_dict.cl_file_path, @@ -326,8 +329,9 @@ function Wing( aero_data, aero_model = load_polar_data(csv_file_path) end - aero_file_path = resolve(get(airfoil_aero_map, section.airfoil_id, "")) - section_aero = isempty(aero_file_path) ? nothing : read_section_aero(aero_file_path) + surface = get(airfoil_surface_map, section.airfoil_id, nothing) + section_aero = isnothing(surface) ? nothing : + read_section_aero(resolve(surface[1]), resolve(surface[2]), resolve(surface[3])) prn && println("Section airfoil_id $(section.airfoil_id): Using $aero_model model") diff --git a/test/airfoil_aero/test_airfoil_aero.jl b/test/airfoil_aero/test_airfoil_aero.jl index 946f140e..b0d8ab5a 100644 --- a/test/airfoil_aero/test_airfoil_aero.jl +++ b/test/airfoil_aero/test_airfoil_aero.jl @@ -5,8 +5,8 @@ using VortexStepMethod.AirfoilAero: KulfanParameters, LeastSquaresFit, ShrinkWra shrink_wrap, fit_kulfan_parameters, kulfan_to_coordinates, neuralfoil_aero, class_function, bernstein_basis, leading_edge_basis, normalize_airfoil -using VortexStepMethod: SectionAero, section_surface, read_section_aero, - write_section_aero +using VortexStepMethod: SectionAero, section_surface, read_section_aero +using VortexStepMethod.AirfoilAero: write_section_aero read_dat_coords(path) = begin x = Float64[]; y = Float64[] @@ -116,7 +116,7 @@ end @test size(cm) == (length(alpha_range), length(delta_range)) end -@testset "SectionAero round-trip and interpolation" begin +@testset "SectionAero CSV round-trip and interpolation" begin alpha_range = deg2rad.([-5.0, 0.0, 5.0, 10.0]) delta_range = deg2rad.([-3.0, 0.0, 3.0]) xc = [1.0, 0.5, 0.0, 0.5, 1.0] @@ -130,15 +130,15 @@ end cf = cp ./ 1000 aero = SectionAero(alpha_range, delta_range, x, y, cp, cf) - path = joinpath(mktempdir(), "aero.npz") - write_section_aero(path, aero) - back = read_section_aero(path) + prefix = joinpath(mktempdir(), "af") + dat, cp_csv, cf_csv = write_section_aero(prefix, aero) + back = read_section_aero(dat, cp_csv, cf_csv) xb, yb, cpb, cfb = section_surface(back, alpha_range[2], delta_range[1]) - @test xb ≈ xc - @test yb ≈ yc + @test isapprox(xb, xc; atol=1e-6) # contour via .dat (8-decimal) + @test isapprox(yb, yc; atol=1e-6) @test cpb ≈ cp[:, 2, 1] @test cfb ≈ cf[:, 2, 1] - @test read_section_aero(joinpath(@__DIR__, "does_not_exist.npz")) === nothing + @test read_section_aero("no.dat", "no_cp.csv", "no_cf.csv") === nothing _, _, cpa, _ = section_surface(aero, alpha_range[3], delta_range[2]) @test cpa ≈ cp[:, 3, 2] From ee75dcd5160b999bffeb37f5a1591fd3e61b931b Mon Sep 17 00:00:00 2001 From: Bart Date: Mon, 20 Jul 2026 00:43:09 +0200 Subject: [PATCH 3/6] Submodule for generating/wringing and main module just reads --- docs/src/functions.md | 8 ++-- docs/src/private_functions.md | 5 +- examples/V3_neuralfoil.jl | 41 +++++----------- src/airfoil_aero/AirfoilAero.jl | 4 +- src/airfoil_aero/airfoil_io.jl | 2 +- src/airfoil_aero/section_aero_gen.jl | 40 +++++++++++----- src/obj_adapter/ObjAdapter.jl | 2 - src/obj_adapter/obj_to_yaml.jl | 23 +++++---- src/obj_adapter/polar_generation.jl | 65 -------------------------- test/airfoil_aero/test_airfoil_aero.jl | 4 +- test/obj_adapter/test_obj_adapter.jl | 25 +++++----- 11 files changed, 80 insertions(+), 139 deletions(-) delete mode 100644 src/obj_adapter/polar_generation.jl diff --git a/docs/src/functions.md b/docs/src/functions.md index 3462b56e..a3257487 100644 --- a/docs/src/functions.md +++ b/docs/src/functions.md @@ -10,13 +10,12 @@ calculate_projected_area ``` ## Surface aero (contour + Cp + cf) tables -The per-section surface aero table type and its npz IO live in the core package; the -generation (from XFoil/NeuralFoil) lives in `AirfoilAero`. +The per-section surface aero table type and its loader live in the core package; the +generation and CSV/dat writing (from XFoil/NeuralFoil) live in `AirfoilAero`. ```@docs SectionAero section_surface read_section_aero -write_section_aero ``` ## Airfoil aerodynamics (AirfoilAero) @@ -45,7 +44,9 @@ neuralfoil_aero generate_aero_matrices generate_polar_from_coordinates generate_polar_from_dat +generate_airfoil_aero generate_section_aero +write_section_aero ``` ## OBJ mesh conversion (ObjAdapter) @@ -57,7 +58,6 @@ CurrentModule = VortexStepMethod.ObjAdapter ``` ```@docs obj_to_yaml -generate_section_polars write_yaml perpendicular_sections plot_airfoils diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index d4bfffc1..610ebe0a 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -74,7 +74,8 @@ calculate_new_aero_data assemble_polar_matrix load_matrix_polar_data read_aero_matrix -write_aero_matrix +read_dat +read_node_table interpolate_matrix_nans! remove_vector_nans generate_polar_data @@ -140,6 +141,8 @@ read_dat_coordinates write_dat write_polar_csv write_polar_matrix_csv +write_aero_matrix +write_node_table flat_plate_cf neuralfoil_contour_solution fill_node_nans! diff --git a/examples/V3_neuralfoil.jl b/examples/V3_neuralfoil.jl index df76ac1e..4e2f2127 100644 --- a/examples/V3_neuralfoil.jl +++ b/examples/V3_neuralfoil.jl @@ -59,13 +59,13 @@ fig_slices_3d = plot_slices_3d(OBJ_PATH; n_slices=N_SLICES, rotation=ROTATION, n_bins=N_BINS, wrap_method=WRAP) GLMakie.save("V3_slices_3d.png", fig_slices_3d) -# Generate per-slice NeuralFoil polars from the shrink-wrapped sections. XFoil is also -# supported (pass `solver=XFoilSolver(...)`) but its panel method does not converge on -# the V3's thin, tube-nosed membrane airfoils, so NeuralFoil is used here. -println("Generating NeuralFoil polars from OBJ mesh...") -generate_section_polars(OBJ_PATH, POLARS_DIR; n_slices=N_SLICES, Re=RE, - rotation=ROTATION, n_bins=N_BINS, wrap_method=WRAP, - solver=NF_SOLVER, delta_range=DELTA_RANGE, verbose=true) +# `obj_to_yaml` generates the whole aero geometry from the OBJ: shape (.dat), polar CSV, +# and Cp/cf tables plus a geometry.yaml. XFoil (`aero_solver=XFoilSolver(...)`) also works +# but doesn't converge on the V3's thin membrane airfoils, so NeuralFoil is used here. +println("Generating V3 aero geometry from OBJ (obj_to_yaml)...") +gen_dir = joinpath("data", "TUDELFT_V3_KITE", "generated_neuralfoil") +nf_yaml = obj_to_yaml(OBJ_PATH, gen_dir; n_sections=N_SLICES, Re=RE, + rotation=ROTATION, wrap_method=WRAP, aero_solver=NF_SOLVER, verbose=true) # Flight conditions v_a = 10.0 @@ -82,29 +82,12 @@ body_cfd = BodyAerodynamics([wing_cfd]) VortexStepMethod.reinit!(body_cfd) solver_cfd = Solver(body_cfd, settings_cfd) -# Build a wing whose every airfoil points at a generated POLAR_VECTORS CSV under -# `polars_subdir`. Repoints a copy of the awesIO geometry and loads it as a normal -# uniform POLAR_VECTORS wing — no manual section surgery needed. -using YAML -yaml_path = joinpath("data", "TUDELFT_V3_KITE", "aero_geometry.yaml") -function build_polar_wing(polars_subdir, out_name) - geom = YAML.load_file(yaml_path) - for (i, airfoil) in enumerate(geom["wing_airfoils"]["data"]) - airfoil[2] = "polars" - airfoil[3] = Dict("csv_file_path" => "$polars_subdir/$i.csv") - end - out_yaml = joinpath("data", "TUDELFT_V3_KITE", out_name) - write_yaml(out_yaml, geom) - wing = Wing(out_yaml; n_panels=50, spanwise_distribution=LINEAR) - refine!(wing) - body = BodyAerodynamics([wing]) - VortexStepMethod.reinit!(body) - s = Solver(body, settings_cfd) - return s, body -end - println("Creating wing with NeuralFoil polars...") -solver_nf, body_nf = build_polar_wing("polars_neuralfoil", "aero_geometry_neuralfoil.yaml") +wing_nf = Wing(nf_yaml; n_panels=50, spanwise_distribution=LINEAR) +refine!(wing_nf) +body_nf = BodyAerodynamics([wing_nf]) +VortexStepMethod.reinit!(body_nf) +solver_nf = Solver(body_nf, settings_cfd) # Compare CFD-polar and NeuralFoil-polar wings against published references # (Poland 2025 RANS CFD and wind tunnel). `plot_polars` sweeps each solver over the diff --git a/src/airfoil_aero/AirfoilAero.jl b/src/airfoil_aero/AirfoilAero.jl index 25f63899..44f815ff 100644 --- a/src/airfoil_aero/AirfoilAero.jl +++ b/src/airfoil_aero/AirfoilAero.jl @@ -29,9 +29,9 @@ export neuralfoil_aero, neuralfoil_section export AbstractAirfoilSolver, XFoilSolver, NeuralFoilSolver export SectionSolution, DeformedSection, deform_section, analyze_section, analyze_sweep export create_2d_polars, generate_aero_matrices, generate_section_aero, lei_poly_coeffs -export write_section_aero +export generate_airfoil_aero, write_section_aero export turn_trailing_edge!, get_lower_upper -export read_dat_coordinates, write_dat, write_polar_csv +export read_dat_coordinates, write_dat, write_polar_csv, write_polar_matrix_csv export generate_polar_from_coordinates, generate_polar_from_dat, resolve_airfoil end diff --git a/src/airfoil_aero/airfoil_io.jl b/src/airfoil_aero/airfoil_io.jl index 1c653efc..ad8cbc9a 100644 --- a/src/airfoil_aero/airfoil_io.jl +++ b/src/airfoil_aero/airfoil_io.jl @@ -30,7 +30,7 @@ function write_dat(filepath::String, name::String, x::Vector, y::Vector) open(filepath, "w") do io println(io, name) for i in eachindex(x) - println(io, @sprintf("%.8f %.8f", x[i], y[i])) + println(io, @sprintf("%.5f %.5f", x[i], y[i])) end end return filepath diff --git a/src/airfoil_aero/section_aero_gen.jl b/src/airfoil_aero/section_aero_gen.jl index a1a09be5..cbfbe666 100644 --- a/src/airfoil_aero/section_aero_gen.jl +++ b/src/airfoil_aero/section_aero_gen.jl @@ -1,15 +1,15 @@ """ - generate_section_aero(solver, base; alpha_range, delta_range, reynolds_number, - crease_frac=0.9, remove_nan=true) -> SectionAero + generate_airfoil_aero(solver, base; alpha_range, delta_range, reynolds_number, + crease_frac=0.9, remove_nan=true) -> (SectionAero, sols) -Build a native-resolution surface table (closed contour + `Cp` + `cf` per node over the -`(alpha, delta)` grid) for a base airfoil (Kulfan parameters) using the given -[`AbstractAirfoilSolver`](@ref). For each `delta` the section is deformed -([`deform_section`](@ref)) and swept over `alpha_range` (radians); each converged angle -contributes its full-contour `(x, y, cp, cf)`. Non-converged points stay `NaN` and, -when `remove_nan`, are filled per node with `interpolate_matrix_nans!`. +Run **one** solver sweep of a base airfoil (Kulfan) over the `(alpha, delta)` grid +(radians) and return both the [`SectionAero`](@ref) (contour + `Cp` + `cf` per node) and +the raw `sols::Vector{Vector{SectionSolution}}` (one inner vector per delta). The `sols` +also carry `cl/cd/cm`, so a caller can write the polar from the same sweep — this is how +`obj_to_yaml` avoids a second sweep. Non-converged points stay `NaN` and, when +`remove_nan`, are filled per node with `interpolate_matrix_nans!`. """ -function generate_section_aero(solver::AbstractAirfoilSolver, base::KulfanParameters; +function generate_airfoil_aero(solver::AbstractAirfoilSolver, base::KulfanParameters; alpha_range, delta_range, reynolds_number, crease_frac=0.9, remove_nan=true) x0, y0 = kulfan_to_coordinates(base) n_alpha, n_delta = length(alpha_range), length(delta_range) @@ -26,7 +26,7 @@ function generate_section_aero(solver::AbstractAirfoilSolver, base::KulfanParame break end end - n_node == 0 && throw(ErrorException("generate_section_aero: no converged angle")) + n_node == 0 && throw(ErrorException("generate_airfoil_aero: no converged angle")) x = fill(NaN, n_node, n_delta) y = fill(NaN, n_node, n_delta) @@ -48,16 +48,32 @@ function generate_section_aero(solver::AbstractAirfoilSolver, base::KulfanParame fill_node_nans!(cf, i) end end - return SectionAero(collect(float.(alpha_range)), collect(float.(delta_range)), + aero = SectionAero(collect(float.(alpha_range)), collect(float.(delta_range)), x, y, cp, cf) + return aero, sols end """ - generate_section_aero(solver, x::Vector, y::Vector; kwargs...) -> SectionAero + generate_airfoil_aero(solver, x::Vector, y::Vector; kwargs...) -> (SectionAero, sols) Convenience: [`shrink_wrap`](@ref) the coordinates and fit base Kulfan parameters ([`LeastSquaresFit`](@ref)) first. """ +function generate_airfoil_aero(solver::AbstractAirfoilSolver, x::Vector, y::Vector; + kwargs...) + xw, yw = shrink_wrap(x, y, ShrinkWrap()) + return generate_airfoil_aero(solver, fit_kulfan_parameters(xw, yw, LeastSquaresFit()); + kwargs...) +end + +""" + generate_section_aero(solver, base_or_coords; kwargs...) -> SectionAero + +Just the [`SectionAero`](@ref) from [`generate_airfoil_aero`](@ref) (drops the raw sols). +""" +generate_section_aero(solver::AbstractAirfoilSolver, base::KulfanParameters; kwargs...) = + generate_airfoil_aero(solver, base; kwargs...)[1] + function generate_section_aero(solver::AbstractAirfoilSolver, x::Vector, y::Vector; kwargs...) xw, yw = shrink_wrap(x, y, ShrinkWrap()) diff --git a/src/obj_adapter/ObjAdapter.jl b/src/obj_adapter/ObjAdapter.jl index 698ca9a0..492ca228 100644 --- a/src/obj_adapter/ObjAdapter.jl +++ b/src/obj_adapter/ObjAdapter.jl @@ -13,7 +13,6 @@ import ..VortexStepMethod include("obj_geometry.jl") include("obj_slice.jl") -include("polar_generation.jl") include("obj_to_yaml.jl") """ @@ -46,7 +45,6 @@ function plot_slices_3d end export obj_to_yaml, resolve_aero_geometry export write_geometry_yaml, write_yaml, airfoils_from_yaml export perpendicular_sections -export generate_section_polars export read_faces, center_to_com!, calculate_inertia_tensor, calc_inertia_y_rotation export plot_airfoil_fit, plot_airfoils, plot_slices_3d diff --git a/src/obj_adapter/obj_to_yaml.jl b/src/obj_adapter/obj_to_yaml.jl index 33afcb5c..c54248ec 100644 --- a/src/obj_adapter/obj_to_yaml.jl +++ b/src/obj_adapter/obj_to_yaml.jl @@ -130,16 +130,21 @@ function obj_to_yaml(obj_path::String, output_dir::String; raw_rel = joinpath("airfoils", "$(j)_raw.dat") csv_rel = joinpath("polars", "$j.csv") try - res = generate_polar_from_coordinates(s.x_fit, s.y_fit, - joinpath(output_dir, csv_rel); Re=Float64(Re), alpha_range, - solver=aero_solver, delta_range, crease_frac, - dat_prefix=joinpath(output_dir, "airfoils", "$j")) - clvals = res isa AbstractVector ? collect(sol.cl for sol in res) : vec(res[1]) + alphas = deg2rad.(collect(Float64, alpha_range)) + deltas = isnothing(delta_range) ? [0.0] : deg2rad.(collect(Float64, delta_range)) + # one solver sweep -> both the SectionAero and the polar + aero, sols = generate_airfoil_aero(aero_solver, s.x_fit, s.y_fit; + alpha_range=alphas, delta_range=deltas, reynolds_number=Float64(Re), + crease_frac) + clvals = collect(sol.cl for sol in sols[1]) all(isnan, clvals) && error("solver produced no converged points") - aero = generate_section_aero(aero_solver, s.x_fit, s.y_fit; - alpha_range=deg2rad.(alpha_range), - delta_range=(isnothing(delta_range) ? [0.0] : deg2rad.(collect(delta_range))), - reynolds_number=Float64(Re), crease_frac) + if isnothing(delta_range) + write_polar_csv(joinpath(output_dir, csv_rel), sols[1]) + else + coeff(f) = [f(sols[jd][ia]) for ia in eachindex(alphas), jd in eachindex(deltas)] + write_polar_matrix_csv(joinpath(output_dir, csv_rel), alphas, deltas, + coeff(s -> s.cl), coeff(s -> s.cd), coeff(s -> s.cm)) + end paths = write_section_aero(joinpath(output_dir, "airfoils", "$j"), aero) dat_rel, cp_rel, cf_rel = (relpath(p, output_dir) for p in paths) write_dat(joinpath(output_dir, raw_rel), "section_$(j)_raw", s.xa, s.ya) diff --git a/src/obj_adapter/polar_generation.jl b/src/obj_adapter/polar_generation.jl deleted file mode 100644 index 917dc6fd..00000000 --- a/src/obj_adapter/polar_generation.jl +++ /dev/null @@ -1,65 +0,0 @@ -""" -Obj-based polar generation: slice a mesh and drive AirfoilAero per section. The -airfoil-general `.dat`/coordinate → CSV entry points live in AirfoilAero. -""" - -""" - generate_section_polars(obj_path, output_dir; n_slices, Re, - alpha_range=-180:1:180, delta_range=nothing, rotation=I, - n_bins=60, wrap_method=ShrinkWrap(), - solver=NeuralFoilSolver(), verbose=true) - -Slice `obj_path` into `n_slices` airfoils ([`perpendicular_sections`](@ref)), -[`shrink_wrap`](@ref) each slice's point cloud into a clean airfoil, and write both its -coordinates (`output_dir/{i}.dat`) and polar CSV (`output_dir/{i}.csv`, via -`generate_polar_from_coordinates`). With a `delta_range`, each deflected shape -is also written to `output_dir/{i}_d{deg}.dat`. The shrink-wrap wraps an open -single-membrane canopy slice (no closed airfoil) into a valid thin cambered airfoil. -`solver` picks the backend — [`NeuralFoilSolver`](@ref) (full ±180° range) or -[`XFoilSolver`](@ref) (use a narrow `alpha_range` where XFoil converges). Pass a -`delta_range` of trailing-edge deflections to write long-format `(alpha, delta)` -`POLAR_MATRICES` CSVs instead of `POLAR_VECTORS`. Pass a `3×3` `rotation` matrix to -reorient the mesh first. -""" -function generate_section_polars(obj_path::String, output_dir::String; - n_slices::Int, Re::Real, - alpha_range=-180:1:180, - delta_range=nothing, - rotation=I, - n_bins=60, - wingtip_distance=0.0, - wrap_method::ShrinkWrap=ShrinkWrap(), - solver::AbstractAirfoilSolver=NeuralFoilSolver(), - verbose=true) - isfile(obj_path) || error("OBJ file not found: $obj_path") - mkpath(output_dir) - verbose && @info "Slicing OBJ mesh at $n_slices positions..." - vertices, faces = read_faces(obj_path) - sections = perpendicular_sections(vertices, faces, n_slices; - n_bins, rotation, wingtip_distance) - verbose && @info "Processing $(length(sections)) airfoil sections..." - for (i, sec) in enumerate(sections) - x, y = sec.x_airfoil, sec.y_airfoil - if isempty(x) || length(x) < 5 - @warn "Slice $i: Insufficient points, skipping" - continue - end - verbose && print(" Slice $i (y=$(round(sec.LE_point[2], digits=3)))... ") - try - xw, yw = shrink_wrap(x, y, wrap_method) - write_dat(joinpath(output_dir, "$i.dat"), "slice_$i", xw, yw) - res = generate_polar_from_coordinates(xw, yw, - joinpath(output_dir, "$i.csv"); Re=Float64(Re), alpha_range, - solver, delta_range, dat_prefix=joinpath(output_dir, "$i")) - clvals = res isa AbstractVector ? (s.cl for s in res) : vec(res[1]) - finite_cl = [v for v in clvals if !isnan(v)] - cl_max = isempty(finite_cl) ? NaN : maximum(finite_cl) - verbose && println("done (CL_max=$(round(cl_max, digits=2)))") - catch e - @warn "Slice $i failed: $e" - continue - end - end - verbose && @info "Polars written to $output_dir" - return nothing -end diff --git a/test/airfoil_aero/test_airfoil_aero.jl b/test/airfoil_aero/test_airfoil_aero.jl index b0d8ab5a..a187ebd7 100644 --- a/test/airfoil_aero/test_airfoil_aero.jl +++ b/test/airfoil_aero/test_airfoil_aero.jl @@ -134,8 +134,8 @@ end dat, cp_csv, cf_csv = write_section_aero(prefix, aero) back = read_section_aero(dat, cp_csv, cf_csv) xb, yb, cpb, cfb = section_surface(back, alpha_range[2], delta_range[1]) - @test isapprox(xb, xc; atol=1e-6) # contour via .dat (8-decimal) - @test isapprox(yb, yc; atol=1e-6) + @test isapprox(xb, xc; atol=1e-4) # contour via .dat (5-decimal) + @test isapprox(yb, yc; atol=1e-4) @test cpb ≈ cp[:, 2, 1] @test cfb ≈ cf[:, 2, 1] @test read_section_aero("no.dat", "no_cp.csv", "no_cf.csv") === nothing diff --git a/test/obj_adapter/test_obj_adapter.jl b/test/obj_adapter/test_obj_adapter.jl index 8afdc4c7..e51368eb 100644 --- a/test/obj_adapter/test_obj_adapter.jl +++ b/test/obj_adapter/test_obj_adapter.jl @@ -3,6 +3,7 @@ using VortexStepMethod.ObjAdapter using VortexStepMethod using VortexStepMethod.AirfoilAero: NeuralFoilSolver using LinearAlgebra +import YAML obj_path = normpath(joinpath(@__DIR__, "..", "..", "data", "ram_air_kite", "ram_air_kite_body.obj")) @@ -45,19 +46,19 @@ obj_path = normpath(joinpath(@__DIR__, "..", "..", @test all(af -> !isempty(af.x), airfoils) end - @testset "generate_section_polars writes per-slice files" begin - secdir = mktempdir() - generate_section_polars(obj_path, secdir; - n_slices=3, Re=5e5, alpha_range=-4:2:4, delta_range=-2:2:2, - solver=NeuralFoilSolver(model_size="medium"), verbose=false) - @test isfile(joinpath(secdir, "1.dat")) - @test isfile(joinpath(secdir, "1.csv")) - # a delta_range writes long-format POLAR_MATRICES + a .dat per deflection - @test occursin("delta", lowercase(readline(joinpath(secdir, "1.csv")))) - @test isfile(joinpath(secdir, "1_d2.dat")) + @testset "obj_to_yaml writes shape + polar + Cp/cf per airfoil" begin + outdir = mktempdir() + yaml = obj_to_yaml(obj_path, outdir; n_sections=3, Re=5e5, + alpha_range=-4:2:4, aero_solver=NeuralFoilSolver(model_size="medium"), + verbose=false) + @test isfile(yaml) + info = Dict(YAML.load_file(yaml)["wing_airfoils"]["data"][1][3]) + @test all(haskey(info, k) for k in ("dat_file", "csv_file_path", "cp_file", "cf_file")) + @test isfile(joinpath(outdir, info["dat_file"])) + @test isfile(joinpath(outdir, info["cp_file"])) + @test isfile(joinpath(outdir, info["cf_file"])) - @test_throws ErrorException generate_section_polars("missing.obj", secdir; - n_slices=3, Re=5e5) + @test_throws ErrorException obj_to_yaml("missing.obj", outdir; n_sections=3, Re=5e5) end @testset "center_to_com! rejects non-triangular faces" begin From 37683e10d595c588ab3df30bd6d566c6e7552f27 Mon Sep 17 00:00:00 2001 From: Bart Date: Mon, 20 Jul 2026 00:43:23 +0200 Subject: [PATCH 4/6] Go back to high fidelity --- src/airfoil_aero/airfoil_io.jl | 2 +- test/airfoil_aero/test_airfoil_aero.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/airfoil_aero/airfoil_io.jl b/src/airfoil_aero/airfoil_io.jl index ad8cbc9a..1c653efc 100644 --- a/src/airfoil_aero/airfoil_io.jl +++ b/src/airfoil_aero/airfoil_io.jl @@ -30,7 +30,7 @@ function write_dat(filepath::String, name::String, x::Vector, y::Vector) open(filepath, "w") do io println(io, name) for i in eachindex(x) - println(io, @sprintf("%.5f %.5f", x[i], y[i])) + println(io, @sprintf("%.8f %.8f", x[i], y[i])) end end return filepath diff --git a/test/airfoil_aero/test_airfoil_aero.jl b/test/airfoil_aero/test_airfoil_aero.jl index a187ebd7..b0d8ab5a 100644 --- a/test/airfoil_aero/test_airfoil_aero.jl +++ b/test/airfoil_aero/test_airfoil_aero.jl @@ -134,8 +134,8 @@ end dat, cp_csv, cf_csv = write_section_aero(prefix, aero) back = read_section_aero(dat, cp_csv, cf_csv) xb, yb, cpb, cfb = section_surface(back, alpha_range[2], delta_range[1]) - @test isapprox(xb, xc; atol=1e-4) # contour via .dat (5-decimal) - @test isapprox(yb, yc; atol=1e-4) + @test isapprox(xb, xc; atol=1e-6) # contour via .dat (8-decimal) + @test isapprox(yb, yc; atol=1e-6) @test cpb ≈ cp[:, 2, 1] @test cfb ≈ cf[:, 2, 1] @test read_section_aero("no.dat", "no_cp.csv", "no_cf.csv") === nothing From 17446ff308206d73531c94c1bb1c487938c7f23c Mon Sep 17 00:00:00 2001 From: Bart Date: Mon, 20 Jul 2026 17:02:57 +0200 Subject: [PATCH 5/6] Derive from spanwise dir --- src/wing_geometry.jl | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/wing_geometry.jl b/src/wing_geometry.jl index c10cf495..3733f869 100644 --- a/src/wing_geometry.jl +++ b/src/wing_geometry.jl @@ -1470,17 +1470,21 @@ function refine_mesh_by_splitting_provided_sections!( section_pair; endpoints=false, reuse_aero_data) - # Apply billowing by rotating chords around LE + # Apply billowing by rotating chords around LE. if billowing_percentage > 0 && idx > start_idx s_l = sections[li]; s_r = sections[li + 1] - le_l = s_l.LE_point; le_r = s_r.LE_point - diff_vec = le_l - le_r + if dot(s_l.LE_point - s_r.LE_point, wing.spanwise_direction) >= 0 + le_neg, le_pos = s_r.LE_point, s_l.LE_point + else + le_neg, le_pos = s_l.LE_point, s_r.LE_point + end + diff_vec = le_pos - le_neg span_len = norm(diff_vec) y_hat = diff_vec / span_len apply_billowing_to_pair!( wing.refined_sections, start_idx, idx - 1, - y_hat, span_len, le_r, + y_hat, span_len, le_neg, s_l.TE_point, s_r.TE_point, billowing_percentage) end From d8ee8f472a46cab3c594bb7a5de8f744e15f57b4 Mon Sep 17 00:00:00 2001 From: Bart Date: Wed, 22 Jul 2026 13:14:01 +0100 Subject: [PATCH 6/6] Work with both spanwise orders --- src/body_aerodynamics.jl | 3 ++- src/panel.jl | 23 ++++++++++++++++++++++- src/wing_geometry.jl | 1 + 3 files changed, 25 insertions(+), 2 deletions(-) diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl index 08c6f0b8..c0b20d03 100644 --- a/src/body_aerodynamics.jl +++ b/src/body_aerodynamics.jl @@ -273,7 +273,8 @@ function reinit!(body_aero::BodyAerodynamics{P, W, T}; panel_props.y_airf[i, :], panel_props.z_airf[i, :], delta, - vec; + vec, + wing.spanwise_direction; remove_nan=wing.remove_nan, init_aero=wing_init_aero ) diff --git a/src/panel.jl b/src/panel.jl index 7207c1b4..7a4d8269 100644 --- a/src/panel.jl +++ b/src/panel.jl @@ -222,6 +222,17 @@ function init_aero!(panel::Panel, section_1::Section, section_2::Section; return nothing end +""" + reinit!(panel, section_1, section_2, aero_center, control_point, bound_point_1, + bound_point_2, x_airf, y_airf, z_airf, delta, vec, spanwise_direction; kwargs...) + +Reinitialize a panel's geometry, horseshoe filaments and aerodynamic interpolations. + +The panel is oriented so its `y_airf` (and the bound vortex `bound_2 -> bound_1`) points +along `+spanwise_direction`, with `z_airf` pointing to the airfoil upper surface. This +makes the aero independent of section ordering: a reversed order would otherwise flip the +normal and make the panel look up its polar at a negated angle of attack. +""" function reinit!( panel::Panel, section_1::Section, @@ -234,12 +245,22 @@ function reinit!( y_airf, z_airf, delta, - vec; + vec, + spanwise_direction; init_aero = true, remove_nan = true ) + flip = dot(y_airf, spanwise_direction) < 0 + if flip + section_1, section_2 = section_2, section_1 + bound_point_1, bound_point_2 = bound_point_2, bound_point_1 + end init_pos!(panel, section_1, section_2, aero_center, control_point, bound_point_1, bound_point_2, x_airf, y_airf, z_airf, delta, vec) + if flip + panel.y_airf .*= -1 + panel.z_airf .*= -1 + end init_aero && init_aero!(panel, section_1, section_2; remove_nan) return nothing end diff --git a/src/wing_geometry.jl b/src/wing_geometry.jl index 3733f869..6b720998 100644 --- a/src/wing_geometry.jl +++ b/src/wing_geometry.jl @@ -895,6 +895,7 @@ function refine!(wing::AbstractWing{T}; recompute_mapping=true, sort_sections=tr sorted || sort!(wing.unrefined_sections; by=_section_sort_key, rev=true) end + n_sections = wing.n_panels + 1 reuse_aero_data = _can_reuse_prior_refined_polar_data(wing, n_sections)