Skip to content

Commit 0f7f62f

Browse files
danieljvickersdanieljvickers
authored
Smooth airfoil surfaces (#1924)
1 parent 1ececa3 commit 0f7f62f

11 files changed

Lines changed: 320 additions & 255 deletions

File tree

‎examples/2D_ibm_airfoil_surface_pressure/case.py‎

Lines changed: 116 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,116 @@
1+
#!/usr/bin/env python3
2+
"""
3+
Inviscid NACA 0012 at M = 0.3, alpha = 2 deg, slip-wall IB, 200 cells per chord.
4+
Surface Cp is compared to a vortex panel solution, Kuethe & Chow (1998) Sec. 5.10,
5+
corrected to M = 0.3 with the Karman-Tsien rule, von Karman (1941) and Tsien (1939).
6+
"""
7+
8+
import json
9+
import math
10+
import os
11+
12+
Ma = 0.3
13+
alpha_deg = 2.0
14+
gamma = 1.4
15+
rho_inf, U_inf, chord = 1.0, 1.0, 1.0
16+
P_inf = rho_inf * U_inf**2 / (gamma * Ma**2)
17+
18+
case = {
19+
# --- Output ---
20+
"run_time_info": "T",
21+
"format": 2,
22+
"precision": 2,
23+
"parallel_io": "T",
24+
"prim_vars_wrt": "T",
25+
"ib_state_wrt": "T", # static IBs only compute/write forces when this is on
26+
"ib_force_wrt": "T",
27+
"ib_force_stride": 20,
28+
# --- Domain: pre-stretch extents; the stretching maps them outward ---
29+
"x_domain%beg": -3.0,
30+
"x_domain%end": 4.0,
31+
"y_domain%beg": -3.0,
32+
"y_domain%end": 3.0,
33+
"m": 1399,
34+
"n": 1199,
35+
"p": 0,
36+
"cyl_coord": "F",
37+
"stretch_x": "T",
38+
"a_x": 15.0,
39+
"x_a": -0.8,
40+
"x_b": 1.8,
41+
"loops_x": 2,
42+
"stretch_y": "T",
43+
"a_y": 15.0,
44+
"y_a": -0.7,
45+
"y_b": 0.7,
46+
"loops_y": 2,
47+
# --- Time stepping ---
48+
"cfl_adap_dt": "T",
49+
"cfl_target": 0.5,
50+
"n_start": 0,
51+
"t_save": 5.0,
52+
"t_stop": 40.0,
53+
# --- Numerics ---
54+
"num_patches": 1,
55+
"num_fluids": 1,
56+
"model_eqns": 2,
57+
"alt_soundspeed": "F",
58+
"mpp_lim": "F",
59+
"mixture_err": "T",
60+
"time_stepper": 3,
61+
"weno_order": 5,
62+
"weno_eps": 1.0e-10,
63+
"weno_Re_flux": "F",
64+
"weno_avg": "T",
65+
"avg_state": 2,
66+
"mapped_weno": "T",
67+
"null_weights": "F",
68+
"mp_weno": "F",
69+
"riemann_solver": 2,
70+
"low_Mach": 2,
71+
"wave_speeds": 1,
72+
"viscous": "F",
73+
"fd_order": 4,
74+
# --- Uniform freestream ---
75+
"patch_icpp(1)%geometry": 3,
76+
"patch_icpp(1)%x_centroid": 0.0,
77+
"patch_icpp(1)%y_centroid": 0.0,
78+
"patch_icpp(1)%length_x": 1.0e3,
79+
"patch_icpp(1)%length_y": 1.0e3,
80+
"patch_icpp(1)%vel(1)": U_inf,
81+
"patch_icpp(1)%vel(2)": 0.0,
82+
"patch_icpp(1)%pres": P_inf,
83+
"patch_icpp(1)%alpha_rho(1)": rho_inf,
84+
"patch_icpp(1)%alpha(1)": 1.0,
85+
"fluid_pp(1)%gamma": 1.0 / (gamma - 1.0),
86+
"fluid_pp(1)%eos": "ideal_gas",
87+
# --- Characteristic far-field boundaries ---
88+
"bc_x%beg": -7,
89+
"bc_x%grcbc_in": "T",
90+
"bc_x%vel_in(1)": U_inf,
91+
"bc_x%vel_in(2)": 0.0,
92+
"bc_x%pres_in": P_inf,
93+
"bc_x%alpha_rho_in(1)": rho_inf,
94+
"bc_x%alpha_in(1)": 1.0,
95+
"bc_x%end": -8,
96+
"bc_x%grcbc_out": "T",
97+
"bc_x%pres_out": P_inf,
98+
"bc_y%beg": -9,
99+
"bc_y%end": -9,
100+
# --- IB: NACA 0012 (m must be > 0, so a negligible camber stands in for 0) ---
101+
"ib": "T",
102+
"num_ibs": 1,
103+
"patch_ib(1)%geometry": 4,
104+
"patch_ib(1)%x_centroid": 0.0,
105+
"patch_ib(1)%y_centroid": 0.0,
106+
"patch_ib(1)%airfoil_id": 1,
107+
"patch_ib(1)%angles(3)": -math.radians(alpha_deg),
108+
"patch_ib(1)%slip": "T",
109+
"patch_ib(1)%moving_ibm": 0,
110+
"ib_airfoil(1)%c": chord,
111+
"ib_airfoil(1)%t": 0.12,
112+
"ib_airfoil(1)%p": 0.4,
113+
"ib_airfoil(1)%m": 1.0e-9,
114+
}
115+
116+
print(json.dumps(case, indent=4))

‎examples/2D_ibm_airfoil_surface_pressure/convergence.png‎

94.8 KB
Loading

‎examples/2D_ibm_airfoil_surface_pressure/readme.md‎

Lines changed: 22 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,22 @@
1+
# Airfoil Surface Pressure
2+
3+
This test documents improvements made to the IBM airfoil patch.
4+
A 2D inviscid NACA 0012 at M = 0.3 and α = 2° is modeled as a slip-wall immersed boundary with 200 cells per chord.
5+
The surface pressure coefficient is compared to a vortex panel solution [1] corrected to M = 0.3 with the Kármán–Tsien rule [2, 3], using the NACA four-digit geometry [4].
6+
The flow is inviscid, attached, and subcritical, so the reference gives $C_l = 0.256$ and $C_d = 0$.
7+
8+
![Surface pressure coefficient compared to the panel reference](convergence.png)
9+
10+
To reproduce the figure, run the case and plot the last saved step:
11+
12+
```bash
13+
./mfc.sh run examples/2D_ibm_airfoil_surface_pressure/case.py # simulate to t = 40
14+
./build/venv/bin/python3 examples/2D_ibm_airfoil_surface_pressure/plot_cp.py # writes cp.png
15+
```
16+
17+
### References
18+
19+
1. A. M. Kuethe and C.-Y. Chow, *Foundations of Aerodynamics: Bases of Aerodynamic Design*, 5th ed., Wiley, 1998. Sec. 5.10.
20+
2. T. von Kármán, "Compressibility effects in aerodynamics," *Journal of the Aeronautical Sciences*, 8(9):337–356, 1941. doi:10.2514/8.10737
21+
3. H. S. Tsien, "Two-dimensional subsonic flow of compressible fluids," *Journal of the Aeronautical Sciences*, 6(10):399–407, 1939. doi:10.2514/8.916
22+
4. I. H. Abbott and A. E. von Doenhoff, *Theory of Wing Sections*, Dover, 1959. Ch. 6.

‎src/common/m_derived_types.fpp‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -326,9 +326,9 @@ module m_derived_types
326326

327327
!> Computed surface grid for a NACA airfoil (simulation-only, not in namelist)
328328
type ib_airfoil_grid
329-
integer :: Np = 0 !< number of surface grid points per surface
330-
type(vec3_dt), allocatable :: upper(:) !< upper surface grid points (1:Np)
331-
type(vec3_dt), allocatable :: lower(:) !< lower surface grid points (1:Np)
329+
integer :: Np = 0 !< number of surface grid points per surface
330+
real(wp), allocatable :: upper(:,:,:) !< upper segments (1:Np-1, vertex 1/vertex 2/normal, x/y), as STL boundary_v
331+
real(wp), allocatable :: lower(:,:,:) !< lower segments (1:Np-1, vertex 1/vertex 2/normal, x/y), as STL boundary_v
332332
end type ib_airfoil_grid
333333

334334
!> User-input parameters for an STL/OBJ immersed boundary model (namelist-safe: scalars + fixed arrays)

‎src/common/m_model.fpp‎

Lines changed: 16 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -782,30 +782,30 @@ contains
782782

783783
!> Determine the levelset distance and normals of 2D models by computing the exact closest point via projection onto boundary
784784
!! edges.
785-
subroutine s_distance_normals_2D(pid, boundary_edge_count, point, normals, distance)
785+
subroutine s_distance_normals_2D(boundary_v, boundary_edge_count, point, normals, distance)
786786

787787
$:GPU_ROUTINE(parallelism='[seq]')
788788

789-
integer, intent(in) :: pid
790-
integer, intent(in) :: boundary_edge_count
791-
real(wp), dimension(1:3), intent(in) :: point
792-
real(wp), dimension(1:3), intent(out) :: normals
793-
real(wp), intent(out) :: distance
794-
integer :: i
795-
real(wp) :: dist_min, dist, t
796-
real(wp) :: v1(1:2), v2(1:2), edge(1:2), pv(1:2)
797-
real(wp) :: edge_len_sq, proj(1:2), norm(1:2)
789+
real(wp), dimension(:,:,:), intent(in) :: boundary_v !< edges (edge, vertex 1/vertex 2/normal, x/y)
790+
integer, intent(in) :: boundary_edge_count
791+
real(wp), dimension(1:3), intent(in) :: point
792+
real(wp), dimension(1:3), intent(out) :: normals
793+
real(wp), intent(out) :: distance
794+
integer :: i
795+
real(wp) :: dist_min, dist, t
796+
real(wp) :: v1(1:2), v2(1:2), edge(1:2), pv(1:2)
797+
real(wp) :: edge_len_sq, proj(1:2), norm(1:2)
798798

799799
dist_min = initial_distance_buffer
800800
normals = 0._wp
801801
norm = 0._wp
802802

803803
do i = 1, boundary_edge_count
804804
! Edge endpoints
805-
v1(1) = gpu_boundary_v(i, 1, 1, pid)
806-
v1(2) = gpu_boundary_v(i, 1, 2, pid)
807-
v2(1) = gpu_boundary_v(i, 2, 1, pid)
808-
v2(2) = gpu_boundary_v(i, 2, 2, pid)
805+
v1(1) = boundary_v(i, 1, 1)
806+
v1(2) = boundary_v(i, 1, 2)
807+
v2(1) = boundary_v(i, 2, 1)
808+
v2(2) = boundary_v(i, 2, 2)
809809

810810
! Edge vector and point-to-v1 vector
811811
edge = v2 - v1
@@ -824,8 +824,8 @@ contains
824824
if (t >= 0._wp .and. t <= 1._wp) then
825825
proj = v1 + t*edge
826826
dist = sqrt((point(1) - proj(1))**2 + (point(2) - proj(2))**2)
827-
norm(1) = gpu_boundary_v(i, 3, 1, pid)
828-
norm(2) = gpu_boundary_v(i, 3, 2, pid)
827+
norm(1) = boundary_v(i, 3, 1)
828+
norm(2) = boundary_v(i, 3, 2)
829829
else if (t < 0._wp) then ! negative t means that v1 is the closest point on the edge
830830
dist = sqrt((point(1) - v1(1))**2 + (point(2) - v1(2))**2)
831831
norm(1) = v1(1) - point(1)

‎src/common/m_patch_geometries.fpp‎

Lines changed: 15 additions & 30 deletions
Original file line numberDiff line numberDiff line change
@@ -78,7 +78,8 @@ contains
7878
integer, intent(in) :: airfoil_id
7979
logical :: is_inside
8080
integer :: k
81-
real(wp) :: f
81+
real(wp) :: segment_fraction, surface_height
82+
real(wp) :: vertex_1(1:2), vertex_2(1:2)
8283

8384
is_inside = .false.
8485

@@ -88,43 +89,27 @@ contains
8889
! if we are in 3D, we must also check the z axis
8990
if (num_dims == 3 .and. (.not. (-0.5_wp*length <= z .and. 0.5_wp*length >= z))) return
9091

91-
! our check branches for the upper and lower half of the airfoil
92+
! find the segment spanning x on the surface bounding this half
93+
k = 1
9294
if (y >= 0._wp) then
93-
! increment the iterator so we know where in the airfoil arrays to look
94-
k = 1
95-
do while (ib_airfoil_grids(airfoil_id)%upper(k)%x < x)
95+
do while (ib_airfoil_grids(airfoil_id)%upper(k, 2, 1) < x)
9696
k = k + 1
9797
end do
98-
99-
! If the values are approximately equivalent, skip the next check
100-
if (f_approx_equal(ib_airfoil_grids(airfoil_id)%upper(k)%x, x)) then
101-
if (y <= ib_airfoil_grids(airfoil_id)%upper(k)%y) is_inside = .true.
102-
else
103-
! check if the y value is below the upper edge of the airfoil
104-
f = (ib_airfoil_grids(airfoil_id)%upper(k)%x - x)/(ib_airfoil_grids(airfoil_id)%upper(k)%x &
105-
& - ib_airfoil_grids(airfoil_id)%upper(k - 1)%x)
106-
if (y <= ((1._wp - f)*ib_airfoil_grids(airfoil_id)%upper(k)%y + f*ib_airfoil_grids(airfoil_id)%upper(k - 1)%y)) &
107-
& is_inside = .true.
108-
end if
98+
vertex_1 = ib_airfoil_grids(airfoil_id)%upper(k, 1,:)
99+
vertex_2 = ib_airfoil_grids(airfoil_id)%upper(k, 2,:)
109100
else
110-
! increment the iterator so we know where in the airfoil arrays to look
111-
k = 1
112-
do while (ib_airfoil_grids(airfoil_id)%lower(k)%x < x)
101+
do while (ib_airfoil_grids(airfoil_id)%lower(k, 2, 1) < x)
113102
k = k + 1
114103
end do
115-
116-
! If the values are approximately equivalent, skip the next check
117-
if (f_approx_equal(ib_airfoil_grids(airfoil_id)%lower(k)%x, x)) then
118-
if (y >= ib_airfoil_grids(airfoil_id)%lower(k)%y) is_inside = .true.
119-
else
120-
! check if the y value is above the lower edge of the airfoil
121-
f = (ib_airfoil_grids(airfoil_id)%lower(k)%x - x)/(ib_airfoil_grids(airfoil_id)%lower(k)%x &
122-
& - ib_airfoil_grids(airfoil_id)%lower(k - 1)%x)
123-
if (y >= ((1._wp - f)*ib_airfoil_grids(airfoil_id)%lower(k)%y + f*ib_airfoil_grids(airfoil_id)%lower(k - 1)%y)) &
124-
& is_inside = .true.
125-
end if
104+
vertex_1 = ib_airfoil_grids(airfoil_id)%lower(k, 1,:)
105+
vertex_2 = ib_airfoil_grids(airfoil_id)%lower(k, 2,:)
126106
end if
127107

108+
! inside if y is between the chord side and the surface height interpolated along the segment
109+
segment_fraction = (x - vertex_1(1))/(vertex_2(1) - vertex_1(1))
110+
surface_height = vertex_1(2) + segment_fraction*(vertex_2(2) - vertex_1(2))
111+
is_inside = merge(y <= surface_height, y >= surface_height, y >= 0._wp)
112+
128113
end function f_is_inside_airfoil
129114

130115
function f_is_inside_ellipse(x, y, length) result(is_inside)

0 commit comments

Comments
 (0)