Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
135 changes: 132 additions & 3 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -7,9 +7,11 @@ You tell the program the flow conditions (fluid, speed, chord), the angle of att
The project comes with two solvers that follow the same workflow:

- **In-house potential solver** – our own panel method, written in pure Python. No external programs needed. It gives lift and pressure distribution, but it ignores viscosity, so it cannot compute drag.
- **XFOIL viscous solver** – runs [XFOIL](https://web.mit.edu/drela/Public/web/xfoil/) in the background. Slower, but it includes the boundary layer, so it also gives drag and can respect a maximum drag coefficient.
- **XFOIL viscous solver** – runs [XFOIL](https://web.mit.edu/drela/Public/web/xfoil/) in the background. Slower, but it includes the boundary layer, so it also gives drag: it can find the lowest drag at a target lift, or the highest lift under a drag limit.

A third script compares the two, so you can see how far potential flow is from the viscous result.
A third script compares the two, so you can see how far potential flow is from the viscous result, and why.

What every symbol, message and number means, and why each value was chosen, is explained in [Reference](#reference).

---

Expand Down Expand Up @@ -239,7 +241,134 @@ It compares Cl and the full Cp distribution of NACA 0012, 2412 and 4412 at three
## Roadmap

- Boundary-layer model for the in-house solver, to estimate drag and transition without XFOIL
- One shared core package for both solvers

---

## Reference

This section explains everything you see on screen and every number the code uses. When a value comes from the original code and was never tuned, we say so.

### Symbols

| Symbol | Meaning |
|--------|---------|
| `m` | Maximum camber, as a fraction of the chord |
| `p` | Position of the maximum camber, as a fraction of the chord |
| `t` | Maximum thickness, as a fraction of the chord |
| NACA `mptt` | Name of the airfoil: first digit = m·100, second = p·10, last two = t·100, each rounded. The name is only a label: the geometry uses the exact values |
| Alpha (α) | Angle of attack, in degrees |
| Cl | Lift coefficient |
| Cd | Drag coefficient. 1 drag count = 0.0001 of Cd |
| Cp | Pressure coefficient. In the panel method Cp = 1 − (V/V∞)² |
| ΔCp | Cp,lower − Cp,upper, both taken at the same x/c. It shows where the lift is produced along the chord |
| x/c | Position along the chord: 0 = leading edge, 1 = trailing edge |
| Re | Reynolds number = speed · chord / kinematic viscosity |
| Mach | Speed / speed of sound |
| Ncrit | XFOIL transition parameter (eⁿ method): the higher it is, the later the boundary layer becomes turbulent |
| BB | Bounding box: the maximum total height the airfoil may have, in metres. It is the distance between the highest and the lowest point of the airfoil (thickness plus camber), measured at 0° angle of attack. Typical use: the space available inside a wing. `OUT` in the BB column means the airfoil is taller than this limit |
| Seed | Number that fixes the random part of the search. Same seed and same inputs give the same run |

### Terminal messages

| Prefix | Meaning |
|--------|---------|
| `[+]` | A step starts or ends |
| `[i]` | Information (computed conditions, exact parameters, seed, where files are saved) |
| `[!]` | Warning or error: read it |
| `PHASE 1 / 2 / 3` | Main steps: 1 = optimization, 2 = geometry and plots, 3 = pressure analysis of the final airfoil |
| `--> Phase 1 / Phase 2` | The two parts of the search inside PHASE 1: Genetic Algorithm, then SLSQP |

### Optimization table

```
| Eval | m | p | t | Cl | Cd | BB | Score |
```

| Column | Meaning |
|--------|---------|
| Eval | Number of the evaluation. The final message calls them "iterations": it is the total number of airfoils tried, including the ones outside the box or failed |
| m, p, t | Parameters of the airfoil being tried |
| Cl, Cd | Result of the analysis. Cd is `-` in the in-house solver (potential flow has no drag) |
| BB | `OUT` = the airfoil is taller than the bounding box. It is not analysed and gets a penalty |
| Score | The number the optimizer tries to make as small as possible (see below). Lower is better |

Special values in the Cl column:

| Value | Meaning |
|-------|---------|
| `Failed` | The analysis did not give a valid result. For XFOIL this includes "did not converge at the target angle" |
| `Timeout` | XFOIL took more than 30 s and was stopped |
| `-` | Not analysed (airfoil outside the bounding box) |

### Score

| Solver / objective | Score | Unit |
|--------------------|-------|------|
| In-house | (10 · (Cl − Cl_target))² | none |
| XFOIL, minimum Cd | Cd · 10⁴ + (10⁴ · Cl excess)², where Cl excess = how far Cl is outside `target ± 0.005` | drag counts |
| XFOIL, maximum Cl | −Cl + (10⁴ · Cd excess)², where Cd excess = how far Cd is above the limit | none (negative is normal) |
| Failed or timeout | 10⁹ + a small term that is lower for airfoils closer to t = 0.12 and m = 0.05 | none |
| Outside the box | 10⁹ + ((height − box) · 1000)² | none |

### Validation line

```
> Alpha = 4.0 deg ... [OK] In-House Cl: 0.7359 | XFOIL Cl: 0.6892 | dCl total +0.0467 = num -0.0078 | Mach -0.0098 | visc +0.0643
```

| Part | Meaning |
|------|---------|
| `[OK]` | All five analyses worked |
| `[FAILED <run>]`, `[TIMEOUT <run>]` | That XFOIL run failed. `<run>` is `inviscid M0`, `inviscid M`, `viscous M0` or `viscous M` (M0 = Mach 0, M = real Mach) |
| XFOIL Cl | The reference: XFOIL viscous at the real Mach |
| dCl total | In-house Cl − reference. Positive = the in-house solver overestimates Cl |
| num | Numerical error of the panel method (in-house vs XFOIL inviscid, both at Mach 0) |
| Mach | Part due to compressibility, which the in-house solver ignores |
| visc | Part due to the boundary layer, which the in-house solver ignores |
| `(Mach first)`, `(visc first)` | A single order was used because an intermediate run failed |
| `(Mach+visc combined)` | Mach and viscosity could only be computed together |
| `(breakdown N/A)` | Only the total is available |

### Values and why

| Value | Where | Why |
|-------|-------|-----|
| Air: ν = 1.46·10⁻⁵ m²/s, a = 340.3 m/s | Fluid choice 1 | Standard atmosphere at sea level (15 °C) |
| Water: ν = 1.00·10⁻⁶ m²/s, a = 1482 m/s | Fluid choice 2 | Water at about 20 °C |
| 50 m/s, 1 m, 4°, Cl 0.8 | Default inputs | Example values from the original code. In air they give Re ≈ 3.4·10⁶ and Mach ≈ 0.15 |
| Bounding box 0.3 m | Default input | From the original code. With a 1 m chord the tallest airfoil in the search range is 0.27 m high, so by default the box never cuts anything |
| 160 panels | Default input | Accurate enough (panel-method error about 1% of Cl against XFOIL inviscid) and fast |
| 0 ≤ m ≤ 0.09 | Search range | Keeps the first digit of the NACA name a single digit |
| 0.1 ≤ p ≤ 0.7, 0.05 ≤ t ≤ 0.25 | Search range | From the original code. They cover the usual NACA 4-digit airfoils; not tuned |
| NACA 2412 | Starting airfoil, always in the initial population | Common reference airfoil, from the original code |
| Population 5, 5 generations | Genetic Algorithm | From the original code: 15 airfoils × 6 generations = 90 evaluations. A short global search is enough because SLSQP refines the result |
| SLSQP: max 50 iterations, ftol 10⁻⁴ | Local search | From the original code; not tuned |
| Step 10⁻⁴ (in-house), 2·10⁻³ (XFOIL) | SLSQP finite differences | XFOIL writes Cl with 4 decimals. With a step of 10⁻⁴ the change of Cl due to p and t is below 0.0001, so the gradient came out as zero and the search stopped early. The panel method has no such limit |
| Weight 10 | In-house score | From the original code. It only scales the score: the best airfoil does not change |
| 10⁹ | Score of failed and out-of-box airfoils | Must be higher than any valid airfoil. The worst valid score is about 4·10⁸ (Cl 2.0 away from the target), so 10⁹ always ranks a valid airfoil first |
| Small term towards t = 0.12, m = 0.05 | Failed airfoils | From the original code: among failed airfoils it prefers usual shapes, so the search moves back towards shapes that converge |
| (excess · 1000)² | Out-of-box penalty | Grows with the excess height, so the search knows which way to go back |
| ± 0.005 | Cl tolerance, minimum-Cd objective | 0.6% of a Cl of 0.8: well above XFOIL's resolution (0.0001) and small enough to keep the target meaningful. An exact Cl would be impossible to hit numerically |
| 10⁴ | Cl penalty weight, minimum-Cd objective | 0.001 of Cl outside the tolerance costs 100 drag counts. Near the target, candidate airfoils differ by about 20 counts in our runs, so going outside the tolerance never pays off |
| 10⁴ | Cd penalty weight, maximum-Cl objective | 1 drag count above the limit costs as much as 1.0 of Cl |
| Cd 0.02 | Default limit, maximum-Cl objective | From the old "maximum tolerable Cd" prompt. At 4° it is too loose to matter (the best airfoil had Cd ≈ 0.007): still to be decided |
| Ncrit = 9 | XFOIL | XFOIL's default, the usual value for an average wind tunnel |
| ITER 500 | XFOIL | Maximum viscous iterations per angle, from the original code |
| 1° steps up to the target angle | XFOIL | XFOIL converges more easily when each angle starts from the previous solution. From the original code |
| 0.1° | Angle tolerance | If XFOIL's closest converged angle is further than this from the target, its Cl belongs to another angle and the evaluation counts as failed |
| 30 s | XFOIL timeout | From the original code. A normal XFOIL run takes about a second |
| 1 – 999 999 | Random seed | Any integer works; the range just keeps it short to type |
| 10⁻⁹ | Baseline tolerance (`tests/`) | Far above round-off (we see about 10⁻¹⁵) and far below a real change (a deliberate small error in the code changed Cp by 2·10⁻⁵) |
| NACA 0012, 2412, 4412 | Validation airfoils | Symmetric, mild camber, higher camber |
| −4° to 10°, step 2° | Validation sweep | From the original code |

### Conditions

- **No valid result:** if no airfoil was analysed successfully, the run stops with `[!] No airfoil could be analysed successfully` instead of reporting a meaningless airfoil.
- **Final checks (XFOIL):** after the search, XFOIL runs once more on the best airfoil. A warning appears if it does not converge, if it converges at another angle, if Cl is outside the tolerance (minimum-Cd objective) or if Cd is above the limit (maximum-Cl objective).
- **Missing libraries:** the script asks before installing them with pip. Answer no and it prints the command to run yourself.
- **Invalid input:** the optimizers stop. The validation script instead falls back to fixed values (Re 10⁶, Mach 0, alpha −4° to 10°), which are not the same as its prompt defaults.
- **Trailing edge:** the standard NACA formula leaves a small gap at the trailing edge (0.0025 of the chord for t = 0.12). The code keeps it as it is; whether to close it is still to be decided.

## Credits

Expand Down
43 changes: 31 additions & 12 deletions inhouse_potential_optimizer/README.md
Original file line number Diff line number Diff line change
@@ -1,20 +1,39 @@
# In-House Potential Optimizer
# In-house potential optimizer

This module represents the mathematical core of the project. The script performs aerodynamic optimization using a **Source-Vortex Panel Method** originally developed in MATLAB and ported here to Python to achieve extremely high performance.
Finds a NACA 4-digit airfoil that reaches a target lift coefficient, using our own source + vortex panel method (potential flow). No external programs are needed: only Python with `numpy`, `scipy` and `matplotlib`.

## Main Features
- **No External Dependencies**: It solves the potential flow field by calculating pressure (Cp) and lift (Cl) distributions completely autonomously, without needing to interface with XFOIL or other executables.
- **Hybrid Optimization**: The search for the ideal airfoil is a two-step process. It starts with a Genetic Algorithm to widely explore the design space and avoid getting stuck in local minima, followed by a sub-millimeter gradient-based refinement (SLSQP method).
- **Bounding Box**: You can impose a maximum wing height. The solver will penalize and discard geometries that, once scaled by the chord, turn out to be too thick for your project constraints.
## Run

## Usage
Simply run the main script from the terminal:
```bash
cd inhouse_potential_optimizer
python run.py
```
Upon launching, the pre-run checks will verify your environment and **automatically download and install** any missing Python dependencies (like `numpy`, `scipy`, or `matplotlib`) in the background.

You will be prompted for a few straightforward parameters (fluid, speed, chord, target angle of attack, target Cl, and maximum thickness constraint). The software will run automatically without interruptions, displaying the optimization progress directly in the terminal.
Press **Enter** to accept the default value of each prompt. If a Python library is missing, the script asks whether to install it with pip.

### Outputs
To keep your workspace clean, all outputs are isolated. The generated airfoil coordinates, pressure vector plots, and detailed textual data (CSV) are automatically saved into a categorized subfolder within a dedicated `Results/` directory.
## What it does

1. Asks for fluid, speed, chord, angle of attack, target Cl, bounding-box height, number of panels and random seed.
2. Searches `m`, `p` and `t` in two phases:
- **Genetic Algorithm** (scipy differential evolution): a random population of airfoils, always including a NACA 2412, explores the whole search range;
- **SLSQP**: a gradient-based method refines the best airfoil found.

Airfoils taller than the bounding box are not analysed (`OUT` in the table).
3. Analyses the best airfoil again and saves geometry, plots and data.

With the default inputs a full run takes a few seconds (about 3 s on our Linux test machine).

The same seed always gives the same run: the seed used is printed at the end and saved in the CSV.

## Limits

- **Potential flow:** no drag, no stall, no boundary layer. Reynolds and Mach are shown but not used. Cl is usually higher than the real (viscous) value: the validation script shows by how much and why.
- **Only Cl is matched:** many airfoils give the same Cl, so different seeds can return different airfoils. Minimising drag needs a boundary-layer model, which is planned.

## Output

`Results/Results_Re<Reynolds>_Alpha<angle>_Cl<target>/` with the airfoil coordinates (`.dat`), plots (`.svg`), the optimization history and the aerodynamic data (`.csv`). The full list of files is in the main README.

## More

What every symbol, message and value means, and why each value was chosen: [Reference](../README.md#reference) in the main README.
69 changes: 54 additions & 15 deletions validation_inhouse_vs_xfoil/README.md
Original file line number Diff line number Diff line change
@@ -1,23 +1,62 @@
# Validation: In-House vs XFOIL
# Validation: in-house vs XFOIL

This section of the project contains the tools and scripts needed to validate our custom potential flow solver.
Shows how far the in-house Cl is from XFOIL viscous, and **why**.

Aerodynamics is a field where mathematical approximations (such as ignoring viscosity in the potential method) must be understood and quantified. The goal of this folder is precisely to provide an isolated environment to compare the outputs of our "inhouse_potential_optimizer" (pressure distributions, theoretical lift) against data generated by more complex viscous solvers like "xfoil_viscous_optimizer" or wind tunnel results.
## Run

This way, anyone using the suite can have a clear idea of the margins of error and the operational limits of our panel model.

## Usage
Run the validation script from the terminal:
```bash
cd validation_inhouse_vs_xfoil
python run_validation.py
```

### Interactive Setup
The script is fully parameterized and interactive. Upon launch, it will prompt you for:
- Operating fluid, design speed, and chord length (which it uses to mathematically calculate the exact Reynolds and Mach numbers).
- The minimum, maximum, and step values for the Angle of Attack (Alpha) sweep.
- The panel density.
Press **Enter** to accept the default value of each prompt (fluid, angle-of-attack sweep, speed, chord, number of panels). XFOIL must be set up as described in `xfoil_viscous_optimizer/README.md`.

## What it computes

For NACA 0012, 2412 and 4412 at every angle, the script runs the in-house panel method (A) and XFOIL four times:

| | Run | Viscous | Mach |
|---|---|---|---|
| A | In-house | no | 0 |
| B | XFOIL | no | 0 |
| C | XFOIL | no | real |
| E | XFOIL | yes | 0 |
| D | XFOIL (reference) | yes | real |

The difference from the reference is split into three parts that always add up exactly to the total:

```
dCl total = A − D = num + Mach + visc
```

- **num** = A − B: numerical error of the panel method (same physics as XFOIL inviscid);
- **Mach**: compressibility, which the in-house solver ignores;
- **visc**: boundary layer, which the in-house solver ignores.

Mach and viscosity affect each other, so their parts are the average of the two possible orders. If an intermediate XFOIL run fails, the script uses the runs that worked and marks the line `(Mach first)`, `(visc first)` or `(Mach+visc combined)`; if XFOIL inviscid at Mach 0 fails it shows only the total, and if the reference run fails no error can be computed for that angle.

A positive dCl means the in-house solver overestimates Cl.

## Typical result

Default inputs (Re 3.4·10⁶, Mach 0.147, 160 panels), Windows, alpha = 8°:

| Airfoil | Total | Numerical | Mach | Viscosity |
|---|---|---|---|---|
| NACA 0012 | +0.054 | −0.008 | −0.015 | +0.077 |
| NACA 2412 | +0.083 | −0.012 | −0.018 | +0.112 |
| NACA 4412 | +0.127 | −0.016 | −0.020 | +0.162 |

The boundary layer explains most of the difference, and its share grows with camber and angle. The numerical error of the panel method is about 1% of Cl.

## Output

`Results/Validation_Re<Reynolds>_Mach<Mach>/`:

- `validation_results.csv`: all Cl values, the viscous Cd, the error parts and the breakdown method used;
- `validation_plot_cl.svg`: Cl vs alpha for each airfoil (in-house, XFOIL inviscid, XFOIL viscous);
- `validation_plot_error.svg`: the error parts vs alpha for each airfoil, with the total. Hollow markers = single-order breakdown.

## More

### Outputs
Once the comparative sweeps are completed, the script generates comparative lift coefficient plots (SVG) and raw data matrices (CSV) to help you visualize the potential flow vs viscous flow discrepancies.
To keep the main directory tidy, all validation exports are saved in a categorized subfolder inside the `Results/` directory (e.g., `Results/Validation_Re..._Mach.../`).
What every symbol, message and value means, and why each value was chosen: [Reference](../README.md#reference) in the main README.
Loading