Use Bezier curve to fit NACA 4-digit airfoil

[yungyuc]

In the current prototype of the naca app, the airfoil is described using sampled edges. I'd like to try to describe the airfoil using Bezier curves, so that generating meshes for the airfoil shapes ca use arbitrarily high resolution.

Progress

• PR Use cubic Bézier curve to sample NACA 4-digit airfoil #384 adds a prototype that samples the 4th-order polynomial profile of the airfoil using both cubic Bezier curves and line segments, but has not done fitting. More details in the comment: Use Bezier curve to fit NACA 4-digit airfoil  #320 (comment)
• PR Refactor the code related to NACA4 airfoil #406 refactors the code added in PR Use cubic Bézier curve to sample NACA 4-digit airfoil #384. More details in the comment: Use Bezier curve to fit NACA 4-digit airfoil  #320 (comment)
• Use cubic Bezier curves and line segments to fit the airfoil profile.

Open for discussions

Add GUI control for changing the airfoil drawing and fitting parameters.

[yungyuc]

PR #384 adds a prototype that samples the 4th-order polynomial profile of the airfoil using both cubic bezier curves:

modmesh/modmesh/gui/naca.py

Lines 113 to 127 in 6e634cd

	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

And straight lines:

modmesh/modmesh/gui/naca.py

Lines 106 to 110 in 6e634cd

	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

But fitting is not done yet.

[yungyuc]

@YenPeiChen07 Could you please help completing all tasks for fitting the airfoil? I may add more as we go. But if the work lasts too long I may create a new issue to track.

If you are not familiar with an optimization algorithm like linear least square, perhaps discussing with @Gene0315 can help.

[yungyuc]

PR #406 refactors the code added in PR #384 and moved the calculating and drawing code from modmesh/gui/naca.py to modmesh/pilot/airfoil/__init__.py. modmesh/gui/naca.py is left for adding interactive GUI code in the future.

A reference to the wikipedia page https://en.wikipedia.org/wiki/NACA_airfoil is added in the code comment:

modmesh/modmesh/pilot/airfoil/__init__.py

Line 49 in d9214b6

See also https://en.wikipedia.org/wiki/NACA_airfoil

The related equations are listed in the wikipedia page. I also copied them in the code comments. The symmetric profile:

modmesh/modmesh/pilot/airfoil/__init__.py

Lines 90 to 100 in d9214b6

	.. math::
	y_t = 5t(0.2969 \sqrt{x} - 0.1260 x - 0.3516 x^2 + 0.2843 x^3
	- 0.1015 x^4)
	Zero (close) trailing edge:
	.. math::
	y_t = 5t(0.2969 \sqrt{x} - 0.1260 x - 0.3516 x^2 + 0.2843 x^3
	- 0.1036x^4)

The asymmetric profile (camber):

modmesh/modmesh/pilot/airfoil/__init__.py

Lines 127 to 148 in d9214b6

	.. math:
	y_c = \frac{m}{p^2}(2px - x^2), \quad 0 \le x \le p \\
	y_c = \frac{m}{(1-p)^2}(1 - 2p + 2px - x^2), \quad p \le x \le 1
	The thickness needs to be perpendicular to the camber line, the upper
	and lower profile should be calculated by:
	.. math:
	x_u = x - y_t\sin\theta , \; y_u = y_c + y_t\cos\theta \\
	x_l = x + y_t\sin\theta , \; y_l = y_c - y_t\cos\theta
	The angle :math:`\theta` should be determined by:
	.. math:
	\theta = \arctan\frac{\mathrm{d}y_c}{\mathrm{d}x} \\
	\frac{\mathrm{d}y_c}{\mathrm{d}x} = \frac{2m}{p^2}(p-x),
	\quad 0 \le x \le p \\
	\frac{\mathrm{d}y_c}{\mathrm{d}x} = \frac{2m}{(1-p)^2}(p-x),
	\quad p \le x \le 1

[yungyuc]

Also in the future, we may consider to move other GUI code from modmesh/view.py to be under the package modmesh/pilot/. A package named as pilot may better suggest the geometry creation and mesh generation we are planning to add in modmesh. The old name view does not have the meaning.