Skip to content

Commit

Permalink
Merge pull request #406 from yungyuc/refactor/naca4-airfoil
Browse files Browse the repository at this point in the history
Refactor the code related to NACA4 airfoil
  • Loading branch information
yungyuc authored Aug 11, 2024
2 parents 7fe55a9 + 8007f61 commit d9214b6
Show file tree
Hide file tree
Showing 5 changed files with 373 additions and 110 deletions.
119 changes: 9 additions & 110 deletions modmesh/gui/naca.py
Original file line number Diff line number Diff line change
Expand Up @@ -29,125 +29,24 @@
Show NACA airfoil shape
"""

import numpy as np

from modmesh import core
from modmesh import view


def calc_naca4_points(number, npoint, open_trailing_edge=False,
cosine_spacing=False):
"""
Returns points in [0 1] for the given 4 digit NACA number string.
Reference from
https://ntrs.nasa.gov/api/citations/19930091108/downloads/19930091108.pdf
and https://github.com/dgorissen/naca
"""

camber = float(number[0]) / 100.0
pos = float(number[1]) / 10.0
thick = float(number[2:]) / 100.0

a0 = +0.2969
a1 = -0.1260
a2 = -0.3516
a3 = +0.2843
a4 = -0.1015 if open_trailing_edge else -0.1036

if cosine_spacing:
xc = 0.5 * (1.0 - np.cos(np.linspace(0.0, np.pi, npoint + 1)))
else:
xc = np.linspace(0.0, 1.0, npoint + 1)

xh = np.sqrt(xc)
x2 = np.power(xc, 2)
x3 = np.power(xc, 3)
x4 = np.power(xc, 4)
yt = 5.0 * thick * (a0 * xh + a1 * xc + a2 * x2 + a3 * x3 + a4 * x4)

if pos == 0:
xu = xc
yu = yt
xl = xc
yl = -yt
else:
xc1 = xc[xc <= pos]
yc1 = camber / np.power(pos, 2) * xc1 * (2 * pos - xc1)
dyc1_dx = camber / np.power(pos, 2) * (2 * pos - 2 * xc1)

xc2 = xc[xc > pos]
yc2 = camber / np.power(1 - pos, 2) * (1 - 2 * pos + xc2) * (1 - xc2)
dyc2_dx = camber / np.power(1 - pos, 2) * (2 * pos - 2 * xc2)

yc = np.concatenate([yc1, yc2])
dyc_dx = np.concatenate([dyc1_dx, dyc2_dx])
theta = np.arctan(dyc_dx)

xu = xc - yt * np.sin(theta)
yu = yc + yt * np.cos(theta)
xl = xc + yt * np.sin(theta)
yl = yc - yt * np.cos(theta)

xr = np.concatenate([xu[::-1], xl[1:]])
yr = np.concatenate([yu[::-1], yl[1:]])
ret = np.vstack([xr, yr]).T.copy()
return ret


def calc_sample_points(points):
# determine the number of sample points for a cubic bezier curve
spacing = 0.01
segments = np.hypot((points[:-1, 0] - points[1:, 0]),
(points[:-1, 1] - points[1:, 1])) // spacing
nsample = np.where(segments > 2, segments - 1, 2).astype(int)
return nsample


def linear_approximate(world, points):
for it in range(points.shape[0] - 1):
world.add_edge(points[it, 0], points[it, 1], 0,
points[it + 1, 0], points[it + 1, 1], 0)
return


def cubic_bezier_approximate(world, points, nsample=None):
Vector = core.Vector3dFp64

if nsample is None:
nsample = calc_sample_points(points)
for it in range(points.shape[0] - 1):
p1 = points[it] + (1 / 3) * (points[it + 1] - points[it])
p2 = points[it] + (2 / 3) * (points[it + 1] - points[it])

b = world.add_bezier([Vector(points[it, 0], points[it, 1], 0),
Vector(p1[0], p1[1], 0),
Vector(p2[0], p2[1], 0),
Vector(points[it + 1, 0], points[it + 1, 1], 0)])
b.sample(nsample[it])
return


def draw_naca4(world, number, npoint, fac, off_x, off_y,
use_bezier=False, **kw):
crds = calc_naca4_points(number=number, npoint=npoint, **kw)
crds *= fac # scaling factor
crds[:, 0] += off_x # offset in x
crds[:, 1] += off_y # offset in y

if use_bezier:
cubic_bezier_approximate(world, crds)
else:
linear_approximate(world, crds)
return crds
from modmesh.pilot import airfoil


def runmain():
"""
A simple example for drawing a couple of cubic Bezier curves.
"""
w = core.WorldFp64()
draw_naca4(w, number='0012', npoint=101, fac=5.0, off_x=0.0, off_y=2.0,
open_trailing_edge=False, use_bezier=True, cosine_spacing=True)
naca4 = airfoil.Naca4(number='0012', open_trailing_edge=False,
cosine_spacing=False)
sampler = airfoil.Naca4Sampler(w, naca4)
sampler.populate_points(npoint=101, fac=5.0, off_x=0.0, off_y=2.0)
if False:
sampler.draw_line()
else:
sampler.draw_cbc()
wid = view.mgr.add3DWidget()
wid.updateWorld(w)
wid.showMark()
Expand Down
38 changes: 38 additions & 0 deletions modmesh/pilot/__init__.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,38 @@
# Copyright (c) 2024, Yung-Yu Chen <[email protected]>
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
# - Redistributions of source code must retain the above copyright notice, this
# list of conditions and the following disclaimer.
# - Redistributions in binary form must reproduce the above copyright notice,
# this list of conditions and the following disclaimer in the documentation
# and/or other materials provided with the distribution.
# - Neither the name of the copyright holder nor the names of its contributors
# may be used to endorse or promote products derived from this software
# without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.


"""
Drawing and visualization sub-system of modmesh.
"""

from . import airfoil # noqa: F401

__all__ = [
'airfoil',
]

# vim: set ff=unix fenc=utf8 et sw=4 ts=4 sts=4:
Loading

0 comments on commit d9214b6

Please sign in to comment.