pyPanair Tutorial#3 Wing/body model

In this tutorial we will perform an analysis of a wing/body configuration.

The model is based on the AGARD-B model1
(The sting is not included, and the aft-body is slightly modified for this tutorial.)

1.Defining the geometry

The model will be composed from 12 networks:

  1. Wing
  2. Nose+Fore-body (upper surface)
  3. Nose+Fore-body (lower surface)
  4. Mid-body (upper surface)
  5. Mid-body (lower surface)
  6. Aft-body (upper surface)
  7. Aft-body (lower surface)
  8. Body base
  9. Wing wake
  10. Body base wake (upper)
  11. Body base wake (lower)
  12. Body-wing wake

The geometry of the model is defined in terms of the body diameter $D$.
For simplicity, it will be set at $D=1$.

1.1 Wing

First, initialize a LaWGS object.


In [1]:
from pyPanair.preprocess import wgs_creator
wgs = wgs_creator.LaWGS("agardb_mod")

Next, we will create a Line object that defines the airfoil at the wing root.
The airfoil of the AGARD-B model is a circular-arc airfoil of 4% thickness.
A csv file has been prepared in advance. The coordinates in circular_arc.csv are normalized, so when using the read_airfoil function, we shall set the expantion_ratio at 2.598. By doing so, the x,z-coordinates of the resulting airfoil will be multiplied by 2.598.


In [2]:
root_airfoil = wgs_creator.read_airfoil("circular_arc.csv", y_coordinate=0.5, expansion_ratio=2.598)

In [3]:
%matplotlib notebook
import matplotlib.pyplot as plt
plt.plot(root_airfoil[:,0], root_airfoil[:,2], "s", mfc="None", mec="b")
plt.xlabel("x")
plt.ylabel("z")
plt.grid()
plt.show()


Since the apex of the nose is the origin, we'll need to shift the coordinates of the root_airfoil, so that the x-coordinate of its LE becomes 4.5. This will be done with the shift method.


In [4]:
root_airfoil = root_airfoil.shift((4.5, 0., 0.))

The shift method will create a copy of root_airfoil and add (4.5, 0., 0.) to the xyz-coordinates of each point in the copy.

Next, we define the tip of the wing. Even though the tip of the wing is a point, when defining it, we need to use the same number of points we used to define the wing root. This means that we need to define a Line with 61 identical points.

This can be done easily with the replace method. By typing


In [5]:
tip_airfoil = root_airfoil.replace(x=7.098, y=2., z=0.)

a copy of root_airfoil will be created. The xyz-coordinates of each point in the copy will be replaced by 2.598, 2., and 0., respectively.

The wing Network will be created and registered to wgs in the same way as tutorials 1 and 2.


In [6]:
wing = root_airfoil.linspace(tip_airfoil, num=30)
wgs.append_network("wing", wing, 1)

In [7]:
wing.plot_wireframe()


1.2 Nose+Fore-body

The nose profile is described by $$r=\frac{x}{3}[1-\frac{1}{9}(\frac{x}{D})^2+\frac{1}{54}(\frac{x}{D})^3]$$ $r$ is the radius and $x$ is the x-coordinate.

The fore-body is a cylinder with a radius of $D$.

First, we define the Line for the nose and fore-body at z=0..


In [8]:
import numpy as np

n_nose = 8  # number of points in the line
x_nose = np.linspace(0., 3., num=n_nose)  # x-coordinates of the line
y_nose = x_nose/3 * (1 - 1/9*(x_nose)**2 + 1/54*(x_nose)**3)  # y-coordinates of the line
nose_line = wgs_creator.Line(np.zeros((n_nose, 3)))  # create an array of zeros
nose_line[:,0] = x_nose
nose_line[:,1] = y_nose

fbody_p1 = wgs_creator.Point(nose_line[-1])  # point at the end of the nose_line
fbody_p2 = fbody_p1.replace(x=4.5)
fbody_line = fbody_p1.linspace(fbody_p2, num=4)  # interpolate fbody_p1 and fbody_p2

nose_fbody_line = nose_line.concat(fbody_line)  # concatenate nose_line & fbody_line

Other Lines that define the nose and fore-body will be created by rotating fbody_line.
The upper/lower surfaces can be created by typing


In [9]:
nose_fbody_up = list()
for i in np.linspace(0, 90, 7):
    line = nose_fbody_line.rotx(rotcenter=nose_fbody_line[0], angle=i)  # create a copy of nose_fbody_line and rotate it around the x-axis
    nose_fbody_up.append(line)
nose_fbody_up = wgs_creator.Network(nose_fbody_up)  # create a Network object from a list of Lines

nose_fbody_low = list()
for i in np.linspace(-90, 0, 7):
    line = nose_fbody_line.rotx(rotcenter=nose_fbody_line[0], angle=i)
    nose_fbody_low.append(line)
nose_fbody_low = wgs_creator.Network(nose_fbody_low)

wgs.append_network("n_fb_up", nose_fbody_up, 1)
wgs.append_network("n_fb_low", nose_fbody_low, 1)

In [10]:
nose_fbody_up.plot_wireframe()



In [11]:
nose_fbody_low.plot_wireframe()


1.3 Mid-body

The mid-body is a cylinder with a hole at the section where the body and wing intersect.

First, we define the upper surface.


In [12]:
wingroot_up, wingroot_low = root_airfoil.split_half()
mbody_line = wingroot_up.replace(z=0.)

mbody_up = list()
for i in np.linspace(90, 0, 7)[:-1]:
    line = mbody_line.rotx((0,0,0), angle=i)
    mbody_up.append(line)
mbody_up.append(wingroot_up)
mbody_up = wgs_creator.Network(mbody_up)
wgs.append_network("mbody_up", mbody_up, 1)

In [13]:
mbody_up.plot_wireframe()


Next, we define the lower surface.


In [14]:
wingroot_low = wingroot_low.flip()
mbody_low = list()
mbody_low.append(wingroot_low)
for i in np.linspace(0, -90, 7)[1:]:
    line = mbody_line.rotx((0,0,0), angle=i)
    mbody_low.append(line)
mbody_low = wgs_creator.Network(mbody_low)
wgs.append_network("mbody_low", mbody_low, 1)

In [15]:
mbody_low.plot_wireframe()


In the original AGARD-B model, the aft-body is a cylinder.
However, in this tutorial it will be a circular truncated cone, for the sake of learning how to attach body wakes.

First, we define the Line for the aft-body at z=0..


In [16]:
aft_body_p1 = wgs_creator.Point(root_airfoil[0])  # TE of the root airfoil
aft_body_p2 = aft_body_p1.shift((1.402, -0.05, 0.))
aft_body_line = aft_body_p1.linspace(aft_body_p2, num=6)

After that, we define the Networks for the upper and lower surfaces.


In [17]:
aft_body_up = list()
for i in np.linspace(0, 90, num=7):
    line = aft_body_line.rotx((0,0,0), angle=i)
    aft_body_up.append(line)
aft_body_up = wgs_creator.Network(aft_body_up)
wgs.append_network("abody_up", aft_body_up, 1)

In [18]:
aft_body_up.plot_wireframe()



In [19]:
aft_body_low = list()
for i in np.linspace(-90, 0, num=7):
    line = aft_body_line.rotx((0,0,0), angle=i)
    aft_body_low.append(line)
aft_body_low = wgs_creator.Network(aft_body_low)
wgs.append_network("abody_low", aft_body_low, 1)

In [20]:
aft_body_low.plot_wireframe()


The body base is a circle. The arc is defined by edge3 from the Networks aft_body_up and aft_body_low
We'll use the edge method to create Line objects from these Networks.


In [21]:
body_base_line_up = aft_body_up.edge(3)  # create a `Line` that corresponds to edge3 of the `Network` `body_base_up`
body_base_line_up = body_base_line_up.flip()  # the order of points are reversed
body_base_line_low = aft_body_low.edge(3)
body_base_line_low = body_base_line_low.flip()
body_base_line = body_base_line_up.concat(body_base_line_low)
body_base_line2 = body_base_line.replace(x=8.5, y=0., z=0.)  # create a `Line` with 13 identical points (the center of the circle)
body_base = body_base_line.linspace(body_base_line2, num=3)
wgs.append_network("bodybase", body_base, boun_type=5)

In [22]:
body_base.plot_wireframe()


Note that when registering a body base network, the boundary type is 5.

1.6 Wakes

For this model we need to attach 3 types of wakes:

  1. Wing wake
  2. Body base wake
  3. Body-wing wake

The wing wake can be created in the same way as tutorials 1 and 2.


In [23]:
wake_length = 2.3093 * 50
wing_wake = wing.make_wake(3, wake_length)
wgs.append_network("wingwake", wing_wake, 18)

The body base wake will be attached to edge3 of the Networks aft_body_up and aft_body_low.
The boundary type will be -18 (not 18).


In [24]:
body_base_wake_up = aft_body_up.make_wake(3, wake_length)
wgs.append_network("bbwake_up", body_base_wake_up, -18)
body_base_wake_low = aft_body_low.make_wake(3, wake_length)
wgs.append_network("bbwake_low", body_base_wake_low, -18)

The body-wing wake is a wake that connects the wing wake and body base wake.
The method make_bodywingwake is used to create this special type of wake.
We will attach the body-wing wake to edge4 of aft_body_up.


In [25]:
body_wake = aft_body_up.make_bodywingwake(4, wake_length)
wgs.append_network("bodywake", body_wake, 20)

In [26]:
body_wake.plot_wireframe()


Note that corner 1 of the body base wake must lie at the intersection of the wing TE and the inboard edge of the wing wake. Also, the boundary type for a body wake is 20.

(Read 3.3.2.5 of Vol.2 of the user manual2 for more information on body wakes.)

Next, we create a stl file and check for errors.


In [27]:
wgs.create_stl()

Finally, we create input files for panin.


In [28]:
wgs.create_wgs()
wgs.create_aux(alpha=5., mach=1.4, cbar=2.3093, span=4., sref=6.928, xref=5.943, zref=0.)

2. Analysis

After creating the geometry of the model, run the analysis.
(If you don't know how to do so, check tutorial 1 or 2.)

3. Validation and visualization

First, we'll validate the results by comparing the aerodynamic coefficients to reference 1.


In [29]:
from pyPanair.postprocess import read_ffmf
read_ffmf()


Out[29]:
sol-no alpha beta cl cdi cy fx fy fz mx my mz area
0 1 5.0000 0.0000 0.24304 0.04082 0.00000 0.01948 0.00000 0.24567 0.00000 0.02639 0.00000 31.09798

According to reference 1, the aerodynamic coefficients for an AGARD-B model at an angle of attack of 5.0 degrees, mach number of 5.0, and a Reynolds number of around 10.4 million is,
$$\begin{align}C_L&=0.224 \\ C_D&=0.0435 \\ C_M&=-0.019\end{align}$$
The lift and drag coefficients are fairly close to the results from reference 1, considering the fact that we've slightly tampered with the shape of aft-body.
On the other hand, the pitching moment coefficient does not match with reference 1. This may be due to multiple factors (e.g. modified aft-body, coarse paneling, etc.).

Visualization can be done in the same way as tutorials 1 and 2.


In [30]:
from pyPanair.postprocess import write_vtk, write_tec
write_vtk(n_wake=4)
write_tec()


n_wake =  4
network 1 shape:  (30, 61, 5)
network 2 shape:  (7, 11, 5)
network 3 shape:  (7, 11, 5)
network 4 shape:  (7, 31, 5)
network 5 shape:  (7, 31, 5)
network 6 shape:  (7, 6, 5)
network 7 shape:  (7, 6, 5)
network 8 shape:  (3, 13, 5)
n_wake =  0
network 1 shape:  (30, 61, 5)
network 2 shape:  (7, 11, 5)
network 3 shape:  (7, 11, 5)
network 4 shape:  (7, 31, 5)
network 5 shape:  (7, 31, 5)
network 6 shape:  (7, 6, 5)
network 7 shape:  (7, 6, 5)
network 8 shape:  (3, 13, 5)
network 9 shape:  (30, 2, 5)
network 10 shape:  (7, 2, 5)
network 11 shape:  (7, 2, 5)
network 12 shape:  (2, 7, 5)

In [31]:
import pandas as pd
from pyPanair.postprocess import calc_section_force
calc_section_force(aoa=5., mac=2.3093, rot_center=(5.943,0,0), casenum=1, networknum=1)

section_force = pd.read_csv("section_force.csv")
section_force


C:\Users\Satoshi\Anaconda3\lib\site-packages\pyPanair\postprocess\section_force.py:48: RuntimeWarning: invalid value encountered in true_divide
  norm = np.fliplr(diff_coord / length[:, np.newaxis])  # norm vector between each vertex
Out[31]:
pos chord cd cl cm
0 0.5000 2.5980 0.021581 0.192287 0.030869
1 0.5517 2.5084 0.021051 0.189387 0.024429
2 0.6034 2.4188 0.021226 0.193268 0.021537
3 0.6552 2.3292 0.021615 0.197942 0.018345
4 0.7069 2.2397 0.022067 0.203079 0.014907
5 0.7586 2.1501 0.022584 0.208542 0.011101
6 0.8103 2.0605 0.023062 0.214417 0.006991
7 0.8621 1.9709 0.023675 0.220937 0.002608
8 0.9138 1.8813 0.024343 0.227773 -0.002179
9 0.9655 1.7917 0.024980 0.234842 -0.007376
10 1.0172 1.7021 0.025545 0.242600 -0.012971
11 1.0690 1.6126 0.026340 0.251225 -0.018927
12 1.1207 1.5230 0.027238 0.260195 -0.025410
13 1.1724 1.4334 0.028041 0.269793 -0.032510
14 1.2241 1.3438 0.028790 0.280585 -0.040245
15 1.2759 1.2542 0.029867 0.292796 -0.048703
16 1.3276 1.1646 0.031040 0.305886 -0.057990
17 1.3793 1.0750 0.032120 0.320299 -0.068239
18 1.4310 0.9854 0.033159 0.336954 -0.079728
19 1.4828 0.8959 0.034686 0.356348 -0.092736
20 1.5345 0.8063 0.036447 0.378042 -0.107413
21 1.5862 0.7167 0.038098 0.403418 -0.124325
22 1.6379 0.6271 0.039794 0.434638 -0.144386
23 1.6897 0.5375 0.042376 0.473100 -0.168576
24 1.7414 0.4479 0.045353 0.521272 -0.198391
25 1.7931 0.3583 0.048937 0.585990 -0.237224
26 1.8448 0.2688 0.053278 0.676350 -0.290197
27 1.8966 0.1792 0.061418 0.799434 -0.362585
28 1.9483 0.0896 0.071547 1.283968 -0.611786
29 2.0000 0.0000 NaN NaN NaN

In [32]:
section_force = section_force.dropna()
plt.plot(section_force.pos / 2, section_force.cl * section_force.chord, "s", mfc="None", mec="b")
plt.xlabel("spanwise position [normalized]")
plt.ylabel("$C_l$ * chord")
plt.grid()
plt.show()


This is the end of tutorial 3.

Reference

  1. Damljanović, D., Vitić, A., Vuković, Đ., and Isaković, J. "Testing of AGARD-B calibration model in the T-38 trisonic wind tunnel," Scientific Technical Review, Vol. 56, No. 2, 2006, pp. 52-62.
  2. Sidwell, K. W., Baruah, P. K., Bussoletti, J. E., Medan, R. T., Conner, R. S., and Purdon, D. J., "PAN AIR: A computer program for predicting subsonic or supersonic linear potential flows about arbitrary configurations using a higher order panel method. Volume 2: User's manual (version 3.0)," NASA CR 3252, 1990.

In [ ]: