On this page is the Simulation Model explained used on the page How Airplanes follow the Curvature of the Earth, Physics Simulation.

Detailed descriptions of the models not described here (Standard Atmosphere, Turbofan Engine, Speed Conversion), further informations and the source code can be found here:

Schematic Diagram of Aircraft Model

Lift Force

To calculate the translational and rotational acceleration of the aircraft, we have to calculate all forces and moments acting on the aircraft.

This equation is valid for ideal gases. In some limits air can be regarded as an ideal gas. The lift force is a vector of magnitude F_L and the direction is pointing perpendicular upwards to the direction of motion SpeedPerp.

  CL: function( a ) {
    // Compute lift coefficient; a = AoA in rad
    return this.Cl0 + this.ClAoA * a;
  },

  CompSpeedsAndBaro: function() {
    :
    AirSpeedModel.Update( this.Altitude, this.TAS );
    this.DynamicPressure = AirSpeedModel.q;
    :    
  }

  ComputeForces: function() {
    :
    var lift = this.DynamicPressure * this.CL(this.AoA) * this.WingArea;
    this.FLift = V2.Scale( this.SpeedPerp, lift );
    :
  }

Elevator Force

The elevator consists of the big trim area and a smaller elevator flap. Both are symmetrically about the Z axis, so they create no force up or down when the angles to the wind are zero. The big elevator trim part is used to trim the airplane for a certain altitude and speed. The trim elevator can be rotated as a whole around a lateral axis. The elevator flap rotates with it but can be angled separately up or down to increase or decreasy the elevator force.

The elevator force is calculated like the lift force but there are 2 separate forces for the trim and the flap part. The flap angle is added to the trim angle.

The elevator trim force acts perpendicular to the surface of the elevator trim. The drag force is accounted for in the overall drag of the airplane.

  CT: function( trimAngle, elevatorAngle ) {
    // Compute trim lift coefficient for trimAngle and elevator angle
    return this.Ct * trimAngle + this.Ce * elevatorAngle;
  },

  CompForces: function() {
    :
    var windNoseAngle = V2.Angle( this.SpeedDir, this.Xaircraft );
    var effTrimAngle = this.TrimAngle - windNoseAngle;
    var effElevAngle = this.ElevatorAngle + this.TrimAngle - windNoseAngle;
    var trimLift = this.DynamicPressure * this.CT( effTrimAngle, effElevAngle ) * this.ElevatorArea;
    var trimLiftVect = V2.Rotate( this.Xaircraft, -this.TrimAngle - Math.PI/2 );
    this.FTrim = V2.Scale( trimLiftVect, trimLift );
    :
  },

Drag Force

The drag force per definition always acts in the opposite direction of motion -SpeedDir.

  CD: function( a ) {
    // Compute drag coefficient; a = AoA in rad
    var cl = this.CL( a );
    return this.Cd0 + this.Cd1 * cl * cl;
  },

  ComputeForces: function() {
    :
    var drag = this.DynamicPressure * this.CD(this.AoA) * this.WingArea;
    this.FDrag = V2.Scale( this.SpeedDir, -drag );
    :
  }

Weight Force

Because the simulation treats the earth as a globe, the gravitational acceleration $\vec g$ depends on the location of the aircraft on earth; $\vec g$ always points towards the center of the earth, which is the origin of the Earth Centered Earth Fixed ECEF coordinate system used by the simulation.

The gravitational acceleration depends on the distance from the center of the earth, and hence from the altitude of the airplane. Note: in the simulation the radius of the earth can be chosen in a wide range to speed up the flight around the curvature. R is initially set to 1/20 of the real radius of the earth.

The weight force is then:

  Gravity: function( h ) {
    return this.g * Math.pow( this.EarthRadius, 2 ) / Math.pow( this.EarthRadius + h, 2 );
  },

  ComputeForces: function() {
    :
    this.FWeight = V2.Scale( this.Zsurface, -this.Mass * this.Gravity(this.Elevation) );
    :
  }

Aerodynamic Moments

Each force acting on the aircraft that does not act on the center of gravity contributes to the pitch moment. To calculate the torque exerted by such a force we first need the vector from the CG to the point where the force is acting on, the displacement vector $\vec r$. Because the displacement vector is in the coordinate system of the airplane, but the forces are given in ECEF coordinates, we have to transform the force vector $\vec F$ into the airplane coordinate system. This yields the transformed force vector $\vec F'$.

The torque is the vector product $\vec M = \vec r \times \vec F'$. The torque is a vector pointing in the direction of the wing axis. We can discard the X and Z component of $\vec M$, because they are always zero. The Y component is our torque value, see function TorqueAroundCG().

  CM: function( a ) {
    // Compute momentum coeffiecient; a = AoA in rad
    return this.Cm0 + this.CmAoA * a;
  },

  TorqueAroundCG: function( forceOrigin, forceVect ) {
    // computes the torque of a force with origin not at the center of gravity
    // around the center of gravity CG.
    // The torque is the cross product (= vector product) of the vector forceOrigin-CG and
    // the forceVect. This product in 3D only has a component in the wing axis Y. 
    // This Y component is just our torque. I can use the V2.VectProd() function to
    // calculate the Y component.

    // transform forceVect from ECEF into aircraft coordinate system
    var forceVect = [ 
      V2.ScalarProd( this.Xaircraft, forceVect ), 
      V2.ScalarProd( this.Zaircraft, forceVect ) ];
    var torque = V2.VectProd( V2.Sub( forceOrigin, this.CenterOfGravity ), forceVect );
    return torque;
  },

  CompForces: function() {
    :
    // compute the various torques and total PitchTorque

    var liftMoment = this.DynamicPressure * this.CM(this.AoA) * this.WingArea * this.MAC;
    var liftTorque = this.TorqueAroundCG( this.CenterOfLift, this.FLift );
    var trimTorque = this.TorqueAroundCG( this.ElevatorPos, this.FTrim );
    var thrustTorque = this.TorqueAroundCG( this.EnginePos, this.FThrust );

    // force daming up/down oscillations depending on PitchRate
    // acting on the horizontal stabilizer
    var horizStabDistFromCG = this.CenterOfGravity[0] - this.ElevatorPos[0];
    var horizStabVertSpeed = this.PitchRate * horizStabDistFromCG;
    var pitchDampingForce = -0.5 * this.AirDensity * Math.pow(horizStabVertSpeed,2) * this.Cpd * this.ElevatorArea;
    var pitchDampingTorque = horizStabDistFromCG * pitchDampingForce;

    // total pitch torque
    this.PitchTorque = liftMoment + liftTorque + trimTorque + thrustTorque + pitchDampingTorque;
  },

Accelerations

For Newton's equation of motion we need the lateral and angular accelerations of the aircraft:

The moment of inertia I is for rotation what the mass m is for linear motion.

  CompAccel: function() {
    // Compute all forces, moments and the acceleration and PitchRate of the aircraft.

    // save previous state (forces..) for calculation of trend vectors

    this.SaveOldState();

    // compute all forces and moments for current aircraft state

    this.CompForces();

    // Compute acceleration from sum of all forces

    var accel = V2.Scale( this.Ftot, 1 / this.Mass );

    var pitchAccel = this.PitchTorque / this.MomentOfInertia;

    this.State.Accel[0] = accel[0];
    this.State.Accel[1] = accel[1];
    this.State.Accel[2] = pitchAccel;
  },

Pressure Altitude, CAS and Mach Speed

The altitude displayed on the Altmeter is derived from the static air pressure outside the airplane. That's why it's called Pressure Altitude. But the airpressure is not always the same for a certain altitude. In a high pressure area the measured air pressure is higher than expected for this altitude. Such pressure differences can be simulated with the ΔP slider, which adds or subtracts a certain pressure from the pressure of the standard atmosphere. The effect on the altimeter is an increase or decrease of the displayed altitude. The airplane will climp or descend to the trimmed altitude if you change the ΔP.

If ΔP is not zero then the elevation above sea level is not equal to the altitude displayed. Some air speeds like CAS and Mach, the speed of sound and some atmospheric parameters depend on the pressure altitude, not the elevation.

Barometric calculation for the model of the International Standard Atmosphere are done by the BaroModel. The conversion of true air speed TAS into calibrated air speed CAS and mach speed are done by the AirSpeedModel.

  CompSpeedsAndBaro: function() {
    // Computes barometric parameters like air pressure, density, temperature
    // speed of sound for current Elevation.
    // Computes the air parameter dependent different air speeds:
    // TAS, Density(Elevation) -> CAS, Mach, speed of sound (a)

    // Use BaroModel to calculate all air parameters for Elevation

    BaroModel.Update( this.Elevation );

    // Compute Pressure Altitude due to non standard pressure DeltaPressure

    if (this.DeltaPressure == 0) {

      this.Altitude = this.Elevation;

    } else {

      var p = BaroModel.p; // standard pressure for Elevation
      // compute pressure altitude
      this.Altitude = BaroModel.HofPressure( p + this.DeltaPressure );
      // compute parameters from pressure altitude (Elevation+DeltaPressure)
      BaroModel.Update( this.Altitude );

    }

    // Use AirSpeedModel to calculate various air speeds from TAS

    AirSpeedModel.Update( this.Altitude, this.TAS );
    this.CAS  = AirSpeedModel.cas;
    this.Mach = AirSpeedModel.mach;
    this.SpeedOfSound = AirSpeedModel.a;

    // Save air parameters in AircraftModel

    this.AirPressure     = AirSpeedModel.p;
    this.AirTemperature  = AirSpeedModel.T - 273.15; // K to °C
    this.AirDensity      = AirSpeedModel.rho;
    this.DynamicPressure = AirSpeedModel.q;

    // compute stall speed and LowestManeuveringSpeed for speed tape red bar

    var Clmax = this.CL( 0.9 * this.MaxAoA );
    var stallSpeed =
      Math.sqrt( 2 * this.Mass * this.Gravity(this.Elevation) /
      (this.AirDensity * Clmax * this.WingArea) );

    AirSpeedModel.Update( this.Altitude, stallSpeed ); // convert to CAS
    this.StallSpeed = AirSpeedModel.cas;

    this.LowestManeuverSpeed = this.StallSpeed * 1.15;

    // compute max operating speed as the smaller of
    // max CAS limit and altitude dependent Mach limit

    // mach limit for altitude
    var mach = AirSpeedModel.AofH( this.Altitude ) * this.LimitMach;
    AirSpeedModel.Update( this.Altitude, mach ); // convert to CAS
    var vmo = AirSpeedModel.cas;

    if (vmo > this.LimitCAS) vmo = this.LimitCAS;
    this.MaxOperatingSpeed = vmo;
  },

Stall Speed

The stall speed is calculated for the display of the minimum selectable and stallspeed limit as a red bar on the speed tape.

The altitude dependent air density ρ(h) can be calculated using the Model of the Standard Atmosphere. Wing area A and lift coefficients C_L rsp. C_L,max can be found in some datasheets of the aircraft.

The extension of slats and flaps increases the wing area and increases the maximum lift coefficient. This reduces the stall speed, i.e. the aircraft can fly more slowly without getting into the stall than if the flaps are not extended. In this simulation landing flaps are not implemented.

If the TAS stall speed is converted into a CAS stall speed, the influence of the air density is eliminated. So the stall speed in CAS is independent on the altitude, but only on the weight of the aircraft.

See Fluggeschwindigkeiten, IAS, TAS, EAS, CAS, Mach for the equations to convert between different air speeds. An online calculator can be found at Air Speeds Calculator (Calculators and Units Converters).

The equation (8) is derived for incompressible gases. Therefore it is a good approximation for the following conditions:

  CompSpeedsAndBaro: function() {
    :
    // compute stall speed and LowestManeuveringSpeed for speed tape red bar
  
    var Clmax = this.CL( 0.9 * this.MaxAoA );
    var stallSpeed = Math.sqrt( 2 * this.Weight * this.g / (this.AirDensity * Clmax * this.WingArea) );

    AirSpeedModel.Update( this.Altitude, stallSpeed ); // convert to CAS
    this.StallSpeed = AirSpeedModel.cas;

    this.LowestManeuverSpeed = this.StallSpeed * 1.15;

    // indicate SPEED and STALL
    if (this.CAS < this.StallSpeed) {
      this.LimitExceeded = 'STALL';
    } else if (this.CAS < this.LowestManeuverSpeed) {
      this.LimitExceeded = 'SPEED';
    }
    :
  }

G-Force Meter

The G-Force Display shows the acceleration the passangers feel due to gravity and airplane maneuvers, expressed as a multiple of the magnitude of gravity $g_o$ at the surface of the earth. The acceleration is a vector, but only its relative magnitude is displayed. If the vector points in the upper hemisphere of the airplanes coordinate system, the sign is positive, else negative.

A value of 1 means the passanger feel exactly the same weight as on the ground. As the airplane flies over a curved surface there is always a small centrifugal acceleration which reduces the G-Force. Because the gravitational acceleration decreases with altitude the G-Force also decreases with altitude.

The dot product of 2 vectors is positive, if the angle between the vectors is less than 90°, else it is negative. This feature can be used to get the sign of G from the up-direction of the airplane and the total acceleration: $\mathrm{sign}( \vec a \cdot \hat z_\mathrm{ac})$.

  Gravity: function( h ) {
    return this.g * Math.pow( this.EarthRadius, 2 ) / Math.pow( this.EarthRadius + h, 2 );
  },
  :
  DrawPfd: function(...) {
    :
    // G-Force
    var gAtElevation = this.Gravity( this.Elevation );
    var gForce = V2.Add( V2.Scale(this.Zsurface, -gAtElevation), V2.Scale( this.State.Accel, -1 ) );
    var gForceSign = V2.ScalarProd( this.Zaircraft, gForce ) > 0 ? -1 : 1;
    var gForceMagnitude = gForceSign * V2.Length( gForce ) / this.g;
    :
  }

Solving the Equation of Motion

Newton's equations of motion give the connection between the acceleration, velocity and position of a moving object.

In complex cases like this aircraft simulation, this equations can not be solved analytically. But we can solve this numerically, step by step.

We saw that the forces acting on the aircraft depend on the current position, orientation, velocity and pitch rate. The State variable of the aircraft simulation contains all this components plus the accelerations.

From the sum of the forces and moments at any time t we can calculated the acceleration and angular acceleration of the aircraft as shown above. To calculate from this the new State, a time step Δt later, we have to numerically integrate the accelerations with respect to time twice. The first integration yields the next velocity and pitch rate, the second integration yield the next position and pitch after a small time step.

Then we calculate from the new State all the forces, moments and accelerations again, and so forth.

For the schedule of the calculations and drawing I implemented the Sim Module. Here is how it is used:

  OnSimUpdate: function() {
    // is called by simulator to calculate the next AC state
    this.UpdateAircraftState();
  },

  UpdateAircraftState: function() {
    :
    this.CompAccel();
    this.Sim.CompNewStates( this.State );  // double integration of State
  },

  OnSimDraw: function() {
    // is called by simulator to draw currect AC state
    this.Update();
  },

  Update: function( field ) {
    :
    ControlPanels.Update();
    this.Draw();
  },
:
var AircraftSim = new Sim( {
  SimObj:        AircraftModel,
  ResetFuncs:    function(sim) { sim.SimObj.OnSimReset(); },
  TimeStepFuncs: function(sim) { sim.SimObj.OnSimUpdate(); },
  FrameFuncs:    function(sim) { sim.SimObj.OnSimDraw(); },
} );

AircraftModel.Init( AircraftGraph, AircraftSim );

AircraftSim.Run( true );