OpenVOGEL/User's guide/Guide 3 Part 2

Batch scan using the console

Introduction

While the classic approach in flight mechanics is to decompose the model and the aerodynamic actions (like for instance "the change of lift on the empennage due to down-wash") and then write complicated equilibrium equations, we will here take advantage of the aerodynamic method we have and propose a few numerical simulations that will only use the resultant aerodynamic forces and moments. This will let us focus on the final behavior of the model, without having to think about the interactions.

These examples will serve as an illustration of what the software can do in real life problem solving. Nothing prevents you from creating your own customized analysis for other situations different than the ones we will expose here. It is hard (or almost impossible) to build software that fits to everyone's needs, so we take for granted that you will be preferring to do so. This is the maximum freedom that we can provide, and it relays on the power and flexibility provided by the Console and the scripts.

Before starting with the standard batch analysis that will help us solve different design aspects of the aircraft, we will see how to run the simplest simulation on the console. The raw data that results from that simulation can be used in customized post-processing, so it is maybe interesting to see.

Steady analysis for a single simulation

In the simplest form of automatic analysis through the Console, you have to create yourself (in Tucan) the set of *.vog files you need. Then you open the Console and enter, for each file:

load
enter file name:
<your file path>
[there comes a message describing what has been loaded]
steady
[the calculation starts]

Scanning the incidence angle

If you want to quickly analyse the aerodynamic loads and the previous method seems too tedious, there are faster alternatives. The first automation procedure we will see is called "alfa scan" and will let you run a series of simulations in one go for a predefined range of incidence angle.

The first thing you need to do is to create in Tucan the model at what you consider is a zero incidence angle. When you are ready, open the Console and load the file:

load
enter file name:
<your file path>
[there comes a message describing what has been loaded]

Next you need to specify that you want to run an alfa scan:

alfa_scan;<alfa_1>;<alfa_2>;<alfa_step>
[the calculation starts for each incidence in the given set]

The starting incidence alfa_1 and the final incidence alfa_2 will always be included, and the alfa_step is the increment that is maintained until alfa_2 is reached.

As result of the scan, an output file is written in the same directory as the input file named <input_file_name>_batch.dat. Here is an example of the output file:

OpenVOGEL alfa scan

A = 1.485717E+000m²
L = 1.485717E+000m
q = 1.159693E+004Pa

# Force coefficients
Alfa   CL             CDi            CDp           
 0.000  7.443804E-001  4.707133E-002  1.146459E-002
 2.000  8.550402E-001  6.360122E-002  1.446163E-002
 ...

# Forces and moments ([N] and  [Nm])
Alfa   Fx             Fy             Fz             Mx             My             Mz            
 0.000  1.004243E+003  1.317944E-014  1.278442E+004  1.558601E+004 -1.781820E+003 -1.224195E+003
 2.000  8.247472E+002  1.538323E-014  1.470960E+004  1.793306E+004 -2.047669E+003 -1.005339E+003
 ...

Longitudinal control performance in straight flight

While the alfa_scan is intended for a simple performance analysis, if the airplane is equipped with control surfaces, it does not tell much about the static equilibrium. In practice, when we compute the polar curve of an airplane with the goal of studying its performance, we need to refer to an equilibrium state. Therefore, the equilibrium and the performance must be analyzed together.

When we design an airplane, one of the many basic questions that arise is: where will the center of gravity be allowed to be so that the aircraft is stable and in static equilibrium throughout the desired flight envelope Some aircraft are designed in such a way that the equilibrium is obtained by displacing the center of gravity (the most common example being hang gliders). However, all other aircraft normally include either a flap or an extra surface (called empennage) equipped with a flap (called elevator) to provide longitudinal control. In this later case, the empennage itself increases the stability of the model while the deflection of the flap provides for control. The size and position of the empennage and elevator play a key role, so the previous question translates into an optimization problem with many variables involved.

One of the things that only a few books about flight mechanics make readers aware of is that the deflection of the elevator not only affects the pitching moment, but also the lift and the drag. For wing-empennage configurations the influence on the lift is less accentuated than for a flying wing, were the effect is very strong. So even for the simplest analysis of longitudinal stability (finding the equilibrium states for straight and leveled flight), we require dealing with many variables. Because the lift and the moment are coupled through the incidence and control deflection angles, finding the equilibrium states is not so straightforward as solving a simple set of uncoupled equations. If only the lift and moment are included in the equilibrium equations, and for both actions the associated coefficients are considered as a linear function of alpha and delta, then at most we can reduce the problem to a system of two linear equations and run only four simulations. You can do that (and save calculation time), but we don't enforce it here to preserve generality. If you want to add propulsive forces, then the problem gets even more complex and it is probably easier to move to a multi-variables graph and brute calculation force.

OpenVOGEL offers you a basic tool to find the equilibrium states in straight flight, and is called alfa_delta_scan. The approach is quite simple: for a range of incidence angles, simulate flight in a range of elevator deflections. Since both angles are limited by nature (the stall and the flap hing constraints respectively), by providing the right ranges the analysis can cover the whole flight envelope. The resulting air loads are then transferred to a scrip in Scilab to make a plot. There, the airloads (lift and pitching moment) are interpolated within the given ranges and contour lines are drawn for the center of gravity and lift coefficient. The resulting graph (in the alpha-delta space) provides the answer to three basic questions:

What is the forwardmost CG position that can still be held at maximum lift?

This state is represented by the CG curve that crosses the maximum alpha and delta values (upper-right corner of the alpha-delta plot for a trailing empennage, and lower-right corner for a forward empennage). If the CG position is limited to that coordinate, the aircraft will still be able to maintain leveled flight at the minimum [stall] speed (or said differently, if the CG exceeds that limit, the aircraft will not be able to fly at the minimum [stall] speed -and imagine how cruel that situation would be for the pilot that is trying to land-). Positioning the CG after that coordinate is a typical mistake for amateur air modelers. Be careful, because you normally will not realize of the issue until you are about to land! (because that is when you want to trade some air speed for extra circulation around the wing -call it CL-). Actually, the most critical condition for this state is when the pitching moment is at its maximum, and this typically occurs when recovering from a dive at maximum lift with all [main wing] flaps down. In OpenVOGEL you can simulate that condition using flow rotation. However we will leave that case for another section.

What is the rearmost CG position that makes the control too sensitive or totally unstable?

This state is represented by the CG curve that is (almost) horizontal. If the CG is located too close to that coordinate, the aircraft will be completely unstable because even the smallest elevator deflection will cause huge variations in the [equilibrium] incidence angle. Now even that is theoretical, because flight dynamics will make the equilibrium simply impossible, so the aircraft will end up following an erratic path (unlike the previous issue, this you would notice from the moment you try to take off). A correct design (for an aircraft that is free of any kind of computer-based control) will limit the CG to a point that is sufficiently far from that coordinate. This CG coordinate is commonly referred to as neutral point or aerodynamic center. Many accidents in civil aviation (and a lot more in air modelling!) have been sadly attributed to the CG exceeding this point, so it is worth understanding it well. Normally, certification norms for civil aircraft will put an upper limit to the g-factor to stick force gradient for the different aircraft categories. This magnitude can be obtained from the slope of the CG curves (but also requires knowing the elevator hing moment and control mechanism properties). For air models or drones, this last value is useless (because the pilot doesn't feel the onboard rigidity of the control) and the limit can be set directly from the slope.

What are the required elevator deflections for a given lift coefficient and center of gravity?

These states are given by the points where a given CG curve crosses the CL curves. Since the CG normally only changes a bit during flight, these points are important to understand the aircraft performance regarding range and autonomy.

Note that because the equilibrium is considered using the Z force component and the Y moment component, the center of gravity must lay along the X axis for the solution to be valid. Also, because the throttle is not included, the throttle force should also be acting on the X axis. For the rest, the method makes no assumptions about how the actions behave and it even includes the contribution of the drag to the moment, so it can be used for any configuration. If you need to add throttle or any other action, you could include it manually on the Scilab script.

To execute an alfa delta scan, you would enter the next command:

alfa_delta_scan;<alfa_1>;<alfa_2>;<alfa_step>;<surface_name>;<region_index>;<delta_1>;<delta_2>;<delta_step>

The first three arguments are the same as for the alfa scan. The following two arguments are the identification of the flap that will be deflected: the name of the lifting surface (case insensitive) and the index of the region (1-based). The last three arguments are similar to the first three ones, but they control the flap deflection range.

If after entering the command the surface or the flap are not valid (not found), the procedure will be aborted.

The next data is an example of an automatically generated Scilab script created by OpenVOGEL alfa delta scan. You can try running it on Scilab to have an idea of how the results look like.

// OpenVOGEL automatic script for alfa-delta scan
// Kernel version: 2.1-2020.05
// Original model: C:\...\scan.vog

X = [    2.00    4.00    6.00    8.00   10.00   12.00]
Y = [  -20.00  -15.00  -10.00   -5.00    0.00    5.00   10.00   15.00   20.00]
CL = [
  4.608154E-001  4.204594E-001  3.785778E-001  3.356219E-001  2.920815E-001  2.484650E-001  2.052767E-001  1.630145E-001  1.221453E-001
  6.028028E-001  5.630714E-001  5.216956E-001  4.791117E-001  4.357995E-001  3.922644E-001  3.490148E-001  3.065488E-001  2.653438E-001
  7.425317E-001  7.035352E-001  6.627821E-001  6.206924E-001  5.777329E-001  5.344030E-001  4.912140E-001  4.486656E-001  4.072419E-001
  8.789988E-001  8.408440E-001  8.008280E-001  7.593512E-001  7.168683E-001  6.738680E-001  6.308602E-001  5.883504E-001  5.468256E-001
  1.011228E+000  9.740175E-001  9.348486E-001  8.941022E-001  8.522159E-001  8.096694E-001  7.669651E-001  7.246132E-001  6.831051E-001
  1.138275E+000  1.102106E+000  1.063892E+000  1.023990E+000  9.828203E-001  9.408495E-001  8.985733E-001  8.564985E-001  8.151244E-001
]
CFz = [
  4.608527E-001  4.203977E-001  3.784443E-001  3.354460E-001  2.918928E-001  2.482961E-001  2.051587E-001  1.629757E-001  1.222121E-001
  6.026229E-001  5.627078E-001  5.212022E-001  4.785446E-001  4.352170E-001  3.917455E-001  3.485989E-001  3.063077E-001  2.653303E-001
  7.414089E-001  7.021558E-001  6.612238E-001  6.190391E-001  5.760721E-001  5.328664E-001  4.898604E-001  4.475741E-001  4.065132E-001
  8.759254E-001  8.374476E-001  7.972121E-001  7.556280E-001  7.131549E-001  6.703403E-001  6.276564E-001  5.854794E-001  5.444539E-001
  1.005720E+000  9.681273E-001  9.287057E-001  8.878469E-001  8.459948E-001  8.036856E-001  7.614793E-001  7.196228E-001  6.786699E-001
  1.128576E+000  1.091976E+000  1.053485E+000  1.013476E+000  9.723749E-001  9.306980E-001  8.890987E-001  8.477344E-001  8.070631E-001
]
CMy = [
 -3.864378E+000 -3.398209E+000 -2.912021E+000 -2.411109E+000 -1.901380E+000 -1.388781E+000 -8.795103E-001 -3.797225E-001  1.048313E-001
 -4.860076E+000 -4.398087E+000 -3.914910E+000 -3.415767E+000 -2.906407E+000 -2.392466E+000 -1.881052E+000 -1.377425E+000 -8.878010E-001
 -5.835478E+000 -5.379238E+000 -4.900740E+000 -4.405003E+000 -3.897652E+000 -3.384237E+000 -2.871843E+000 -2.366306E+000 -1.873239E+000
 -6.782238E+000 -6.333272E+000 -5.861016E+000 -5.370292E+000 -4.866604E+000 -4.355576E+000 -3.843656E+000 -3.337756E+000 -2.842890E+000
 -7.689027E+000 -7.248774E+000 -6.784290E+000 -6.300191E+000 -5.801808E+000 -5.294886E+000 -4.785464E+000 -4.280343E+000 -3.785410E+000
 -8.549503E+000 -8.119354E+000 -7.664151E+000 -7.188262E+000 -6.696827E+000 -6.195587E+000 -5.690682E+000 -5.188234E+000 -4.694368E+000
]
clf
title("Stability plot", "fontsize", 4)
xlabel("alpha [degrees]", "fontsize", 3)
ylabel("delta [degrees]", "fontsize", 3)
xgrid(3)
legends(["iso-CL", "iso-Xcg"], [2,5], "lr")
// CL countour lines
N_CL = 30
Stl_CL = 2*ones(1,N_CL)
contour2d(X, Y, CL, N_CL, Stl_CL)
// Expand CL And CM to refine Xcg
CFz_Spline = splin2d(X, Y, CFz)
CMy_Spline = splin2d(X, Y, CMy)
X_Int = linspace(  2.00,  12.00, 100)
Y_Int = linspace(-20.00,  20.00, 100)
[X_Grid,Y_Grid] = ndgrid(X_Int, Y_Int)
CFz_Int = interp2d(X_Grid, Y_Grid, X, Y, CFz_Spline)
CMy_Int = interp2d(X_Grid, Y_Grid, X, Y, CMy_Spline)
Xcg_Int = CMy_Int./ CFz_Int
N_Xcg = 45
Stl_Xcg = 5*ones(1,N_Xcg)
contour2d(X_Int, Y_Int, Xcg_Int, N_Xcg, Stl_Xcg)

Circular path recovery performance

In the alpha-delta analysis we focused our attention in the longitudinal control characteristics during steady straight flight. However, we mentioned in that section that there is a more restrictive forward-most limit for the CG, and that this happens when the airplane rotates following a circular path (in practice: when it recovers from a dive) at low speed. Imagine the airplane is flying on approach close to the stall airspeed, and therefore the lift coefficient and incidence angle are at the maximum allowed. Now imagine that the airplane needs to go around as soon as possible in an emergency situation. The pilot will push the gas lever to full throttle and gradually pull the elevator up at full [permitted] deflection, trying to bend the trajectory upwards as much as possible. How will the aircraft behave for a given mass configuration?

Simulation approach

OpenVOGEL lets you simulate the rotation of the airplane by adjusting the stream angular velocity (the vector called Omega). Now lets figure out a simple mathematical description of this situation. The aircraft moves in the XZ plane, so some components of the velocity and angular velocity are null:

$\omega _{x}=\omega _{z}=0$

$V_{y}=0$

Then, the curvature of the trajectory depends on the angular velocity to speed ratio:

$\kappa ={\frac {\omega _{y}}{V}}$

with the airflow velocity as:

$V={\sqrt {V_{x}^{2}+V_{z}^{2}}}$

Since the kinematics of the flow will be the same for a same curvature, we can fix a reference velocity and adjust only the angular velocity to obtain a range of curvatures (different circular paths). This is the basis of what we call omega_scan.

With the aid of a free body diagram (see picture) we can write the equilibrium equations of the aircraft at the lowest point of the trajectory. In the vertical axis aligned with the lift force (not to be confused with the Z axis!) we have:

$L=m(g+\kappa V^{2})$

Note that here we did not include the propulsive forces. Around the Y axis we have the moment equilibrium equation:

$M_{y}+x_{cg}F_{z}=0$

Note that we have simplified the dynamics by introducing a steady circular path, however the situation will still resemble the original problem.

From the moment equation we can directly obtain the CG coordinate for each circular path. For this simplified approach, this coordinate will not depend on the velocity nor on the mass, so to make it clear, there is a functional relationship between the curvature and the CG coordinate.

From the force equilibrium equation we can obtain the necessary aircraft velocity for a given mass. Here we have to take into account that the lift depends on the dynamic pressure, and so it is also a function of the square of the velocity. By doing some arithmetic hocus-pocus, we obtain:

$V={\sqrt {\frac {2mg}{\rho SC_{L}-2m\kappa }}}$

and the associated g-factor (the total acceleration as fraction of the gravitational acceleration) will be:

$n=1+{\frac {\kappa V^{2}}{g}}$

What does this all mean? If you got lost in the math (it can happen), this means that by taking the results of our omega_scan set of simulations at maximum lift and maximum elevator deflection, we can draw a diagram of contour lines where the maximum g-factor is plotted for a range of mass and a range of CG coordinates. If we then pick a minimum target g-factor for high incidence maneuvering, the plot will give us the forward-most allowed CG position for a given mass. Extremely useful, isn't it? This also means that by running an alfa_delta_scan and an omega_scan, you will be able to understand most of the longitudinal control of your model at [quasi] steady flight (not the dynamic response! we have not been there yet!).

We can bring an example in numbers to make it clearer. In the model of the depicted diagram, if the landing mass is 1200kg and we still want to be able to recover at 1.3g, the CG should not be placed behind X=0.94. The more gees we want, the closer the CG must be placed to the aerodynamic center.

Command

The omega scan analysis requires calling the omega_scan command (no surprise). A typical script for this simulation would be as follows:

readback
load;C:\...\omega_model.vog
set_velocity;80
set_alfa;12
set_delta;empenage;1;20
set_altitude;1000
pause
omega_scan;0.4;10;2000;4000;50

As you see, the script allows you to easily adapt the incidence and deflection angles, which are not part of the arguments of the omega_scan function. The signature of this functions is as follows:

omega_scan;<Maximum Omega>;<Number of simulations>;<Minimum mass>;<Maximum mass>;<Mass range refinement>

The first argument is the maximum Y component of the angular velocity and the second one the number of aerodynamic simulations (which equals the number of evaluated CG coordinates). The scan will cover the range from 0 rad/s (straight path) to the provided value. The next two following arguments are the minimum and maximum mass respectively, and the last one is the number of control points in the mass range. As result, the command will return the Scilab script that can be executed in the Scilab environment to obtain a design diagram. The name of the script file is equal to the name of the original OpenVOGEL project file, with the string _script appended and a new .sce extension.

Ultimate state (sharpest recovery)

Maybe you have already noticed that the velocity equation is a bit tricky, because it has a singular point and might adopt complex solutions above a certain path curvature or mass. In fact, this means that the aircraft will not be able to exceed a certain path curvature for a given mass (or put the other way around, it will not be able to exceed a certain mass if you want it to be able to pull up with a certain curvature). This is logic, since the lift is not unlimited, and the $C_{L}$ value that we provide is the maximum for any possible air speed. So you will have to be careful with the range of mass and curvature (in the form of angular velocity and airspeed) you provide.

If you provide an unrealistic mass range, the Scilab script will warn you that the maximum range has been exceeded and it will reduce the upper limit to 90% of the critical mass. That way you don't get to see the mess on the diagram.

Shortcomings

Lets now be critic about the shortcomings of this approach:

The center of rotation is not exactly at the CG, it is at the origin, so this can affect the pitching moment. So when you run this simulation, try to locate the origin close to where you think the CG will be located (you can loop to refine if you like).
The real path in recovery is not perfectly circular and the speed is not constant. A hard recovery involves an unsteady transit, during which the aerodynamic forces evolve differently.
The aerodynamic forces are obtained for a single reference velocity (the Reynolds number is constrained).
We didn't include the propulsive forces (although for most situations the effect might be negligible).

These negative points do not mean that this analysis is useless. For most practical applications, it will give you a fast an relatively accurate description of the aircraft performance. In principle, it should be better than most traditional mathematical descriptions.