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:
Mass: Aircraft mass.
CG: Center of Gravity position.
N1: Rotation rate of the big fan of the turbo fan engine as % of maximum N1.
Trim/El: Elevator Trim and Elevator Flap positions.
WingA: Wing Area.
MInert: Moment of Inertia. This is the tendency of the airplane to resits changes in rotation around the wing axes. Moment of Inertia is the rotation equivalent of Mass for translation.
PFD: Primary Flight Display. It is composed left to right of: Speed Tape, Attitude Indicator, Altimeter, Vertical Speed Indicator.
TAS: True Air Speed. Speed of the aircraft with respect to the air. This is displayed at the top right of the Speed Tape.
CAS: Calibrated Airspeed. This is the speed the aircraft 'feels'. Due to less density at heigher altitudes this value is always less than TAS. Because all speed limits are expressed in CAS this is the main speed displayed on the Speed Tape on the PFD.
Mach: Mach Speed. Above the Crossover altitude of about FL280 ≈ 28,000 ft (depending on speed) the selected speed is expressed in Mach, i.e. % of the speed of sound, instead of CAS. Mach speed is displayed below the Speed Tape.
Alt: Pressure Altitude in feet. This altitude is displayed on the Altimeter on the PFD.
VS: Vertical Speed in feet per minutes. This speed is indicated by a pointer needle on the Vertical Speed Indicator on the PFD. The exact value is displayed additionally below the Altimeter.
G Disp: GForce Display. This is the force the passangers feel as weight relative to standard gravity. A value less than 1 means the passanger feel lighter than on earth. The GForce is displayed below the Attitude Indicator.
AoA: Angle of Attack is the angle the wing profile builds with the wind direction. The AoA is one of the major components of the lift. If the AoA exceeds a certain limit (about 15°), the airplane stalls.
Pitch: is the angle the airplanes roll axes build with the horizon. Only on level flight (VS = 0) this is the same angle as AoA.
Dens, P, T, a: Air Denisity in kg/m^{3}, Pressure in hPa, Temperature in °C and speed of sound in m/s. T and P are displayed below the Attitude Indicator. Air Density and Speed of Sound are displayed above the PFD. The Speed of Sound (a) depends on the Air Density and temperature and hence the Altitude of the aircraft.
Lift, Drag, FTrim: Lift, Drag and Elevator Trim Forces.
Weight: Aircraft Weight Force = Mass times Gravitational Acceleration (g).
Thrust: Thrust Force delivered by the engines. Thrust decreases with increasing TAS and increasing altitude due to decreasing air density.
Forces: Vector Sum of all forces. The simulator uses Earth Centered Earth Fixed coordinates.
Accel: Acceleration of the Aircraft in all directions. This is calculated using Newton's law F = m · a. Integrating the Acceleration over time gives the Velocity Vector Veloc.
MLift: Lift Moment. Depending on the amount of Lift a wing creates, it also creates a Lift Moment around a wing axis (Pitching Moment). Choosing a special position for the center of lift the lift moment can be canceled.
MDmp: Pitch Damping Moment. Depending on how fast the aircraft rotates around the wing axes (Pitch Rate), a moment is created that acts against this rotation and damps it.
Mtot: Total of all Moments acting in the wing axis (Pitching Moments). This includes Moments from the Thrust, Weight, and Trim Forces wich do not act through the Center of Lift.
PchAcc: Pitch Acceleration. This is the rate at which the aircrafts rotation arount the wing axes accelerates. Integrating over time gives the Pitch Rate PchRt.
Veloc: Velocity Vector in ECEF coordinates. The component vertical to the surface of the earth gives the Vertical Speed VS. The component in the direction of motion gives the true air speed TAS. The angle between the Velocity Vector and the longitudinal axes is the Angle of Attack AoA. Integrating over time gives the Position of the aircraft Pos.
Pos: Aircraft Position in Earth Centered Earth Fixed ECEF coordinates. The distance from the surface of the earth (distance from the center of the earth minus radius of the earth) is the Elevation of the aircraft. On standard pressure (ΔP = 0) Altitude is the same as Elevation. The Position also defines the direction of the Gravity Vector g.
PchRt: Pitch Rate indicates how fast the aircraft is rotating around the wing axis. Integrating over time gives the Pitch Angle Pitch.
Pitch: Pitch Angle is the angle between the longitudinal axis and the surface of the earth. It is displayed in the Attitude Indicator as the artificial horizon scale. Additionally it is displayed above the PFD.
To calculate the translational and rotational acceleration of the aircraft, we have to calculate all forces and moments acting on the aircraft.
(1) 
 
with  
where^{'} 

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 ); : }
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.
(2) 
 
with  
where^{'} 

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 ); : },
(3) 
 
with  
and  
where^{'} 

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 ); : }
Because the simulation treats the earth as a globe, the gravitational acceleration
(4) 
 
where^{'} 

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:
(5) 
 
where^{'} 

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) ); : }
(6) 
 
with  
where^{'} 

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
The torque is the vector product
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 forceOriginCG 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; },
For Newton's equation of motion we need the lateral and angular accelerations of the aircraft:
(7) 

 
where^{'} 

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; },
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; },
The stall speed is calculated for the display of the minimum selectable and stallspeed limit as a red bar on the speed tape.
(8) 
 
where^{'} 

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:
Altitude h (→ ρ(h))  < 3000 m, < 10,000 ft 
Speed v  < 130 m/s, < 470 km/h, < 250 kt 
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'; } : }
The GForce Display shows the acceleration the passangers feel due to gravity and airplane maneuvers, expressed as a multiple of the magnitude of gravity
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 GForce. Because the gravitational acceleration decreases with altitude the GForce also decreases with altitude.
(9) 
 
width  
and  
where^{'} 

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 updirection of the airplane and the total acceleration:
Gravity: function( h ) { return this.g * Math.pow( this.EarthRadius, 2 ) / Math.pow( this.EarthRadius + h, 2 ); }, : DrawPfd: function(...) { : // GForce 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; : }
Newton's equations of motion give the connection between the acceleration, velocity and position of a moving object.
(10)  
(11)  
(12)  
where^{'} 

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 );