We write a simple * game physics of the aircraft
* - in three dimensions.
Warning: further reasoning may well be wrong, my experience is limited to playing simulators and the course of theoretical mechanics. Knowledge of aerodynamics and game physics engines is very scarce. A picture for attracting attention is a photograph, not a screenshot.
“What could be simpler than a plane? The lift is proportional to the square of speed, the engine pulls forward, everything is simple ”- such a thought came to my mind in the summer, and I sat down to write a game. The summer passed, several rakes were collected, and the list of what I planned to add to the project grew very much.
')
In this article I will try to write about the physical component, without being distracted by the graphics and other things. I hope this will help someone - there is not a lot of information on this topic on the Internet.
Examples of working and not-so-code will be on Scala (to understand the essence, the language is not necessary to know). In addition, I use classes for vectors, matrices, and quaternions from libgdx. If someone has questions on the features of their implementation - the engine code is open. For convenience, vectors have been added methods of the type + =, * =, as well as +, -, the benefit in scala is so possible. It also turned out to be convenient to add methods: = for assignment by value.
Code:
Equivalent to
With the exception of quaternions:
The second option is much more familiar to me, and I never use set. I will not lay out my code yet, there is a little trash, which from time to time is being overwritten beyond recognition.
Write the class for the model position:
As you might guess, the class will have two final fields describing the position of the center of mass of the model and its orientation relative to the world.
You can get the matrix in libgdx like this:
We introduce a class for the derivative:
This class will be convenient not only to describe the speed (first derivative), but also for accelerations and forces.
Stop. How did the quaternion become a vector? Indeed, the angular velocity, acceleration and momentum are usually described by vectors. The essential difference between the angular velocity and orientation is that the orientations are “looped around”, turning two pi is equivalent to “zero” turning. On the contrary, the angular velocity is in principle unlimited.
You can enter the operation of the logarithm of the quaternion q- vector v, which is directed along the direction of the axis of rotation, and its length is equal to the angle of rotation in radians. The exhibitor is the inverse operation. q == exp (v * 2) == exp (v) * exp (v)
The mapping of vector-> quaternion is unambiguous, but the opposite is not. The rotation angle alpha in time dt can correspond to the angular velocity (alpha + 2 * pi * n) / dt, where n is any integer.
Obviously, during the time dt at the angular velocity w, the rotation is q = exp (w * dt).
Code:
What is an airplane flight from an abstract point of view? By solving a system of second-order differential equations! Set the state of the aircraft:
The class for the real state of the aircraft depends on the features of the model and the degree of its physical development, but it will implement this interface, necessary for drawing the aircraft on the screen.
1) the calculation of the forces acting on the aircraft in a given state
2) numerical solution diffura
The first part, as the most interesting, we leave at the end of the article.
For the graphic part, we write the interface:
I will not give the implementation, but the essence is simple - a sequence of states is stored inside with times t1, t2, t3, etc., t (n + 1)> t (n). States are completed if necessary, in the getState method, the next two are interpolated. Thus, it is possible, for example, to read physics 10 times per second and at the same time observe smooth motion at 60 fps.
We write the following interface:
The angular momentum is L = I * w, and L and w (angular velocity) are transformed as vectors. Thus, in the transformation L '= qL, w' = qw it turns out:
L '= I' * w '
qL = I '* qw
L = q ^ (- 1) * I '* q * w
We get I = q ^ (- 1) * I '* q, or I' = q * I * q ^ (- 1).
The transformation w '= position.rot * w translates the angular velocity from the local coordinate system to the global one.
The getForceAndMoment method will be discussed later, it calculates the forces and torque acting on the aircraft.
I do not very well imagine how to accurately calculate the movement of a model that moves and rotates with accelerations, so the simplest method was chosen with a fixed step of 20ms.
More information about the rotation can be found in this article: habrahabr.ru/post/264099 . Honestly, I am not a fan of tensors, I just took a formula from there in vector form in order to obtain angular acceleration. Calculations are made in the coordinate system of the world. By the way, when external forces were disconnected, I was able to observe a movement quite similar to the Janibekov effect.
The plane must be managed:
Values ​​are cut to the interval [-1, 1], engine thrust from 0 to 1.
Well, let's move on to the most important part - we will find strength. Here is a little Terra incognita, my knowledge of aerodynamics is very superficial. Errors and inaccuracies are possible.
The first thing that comes to mind is lifting force. In order not to create garbage, on the Internet, a reference book of aviation profiles was found with charts of lift coefficient depending on the angle of attack. The essence turned out to be quite simple - Cy (lift coefficient) grows quite linearly to critical angles, reaches about one, and then the flow from the wing breaks down, and the lift begins to decrease. Also, the graph of the coefficient for the abstract wing can be viewed in the English Wikipedia :
Then I was trapped by a rake - if you read another article on Wikipedia on frontal resistance, you can see that there is some kind of inductive resistance. The surprise is that the lifting force is considered to be in the direction perpendicular to the direction of speed, and not perpendicular to the surface of the wing (as I thought). Since the difference in air pressure above and below the wing still leads to a force perpendicular to the surface of the wing, the projection of this force on the direction opposite to the movement is non-zero.If I understand correctly, this is inductive force.
And no. See comment below. Further text and code leave unchanged.
If we assume that the lift is directed upwards in the reference system of the aircraft, then the inductive force does not seem to be necessary - it has already been taken into account. The orientation of the axes is the same as in openGL:
Ox - right
Oy - up
Oz - back
In addition to lift, you will need air resistance force:
I do not have sufficient knowledge in aerodynamics. The main goal is simplicity of formulas and, if possible, adequate aircraft behavior for flight angles of attack and glide. 0.5f - the consequences of divider 2 in the formulas. 0.1f - the consequences of adjusting the coefficients.
The model is as simple as possible: no propeller steps, let the engine spend all its power to accelerate the aircraft. No bonuses to the moment of inertia, no torque when changing the number of revolutions. However, the revolutions either. Power = force * speed. To prevent the plane from flying up like a rocket, we will limit the maximum force (using the minimum speed limit).
There is an interesting point - the air accelerated by a screw goes directly to the control planes of the tail, and the plane, in principle, is controlled by the tail a little even on the runway. In addition, the air resistance is trying to twist the plane in the opposite direction of the screw rotation. And up to the heap - the engine has a moment of inertia, when increasing / decreasing the speed of rotation, the plane will also “twist” a little. I will neglect all this ...
As with the wing, a familiar multiplier appears with a square of velocity and air density:
There is no symmetry in pitch, the aircraft (and the pilot) are much better tolerated by positive overloads than negative ones. In addition, the aircraft itself is stable (if it is not the Su-47) and tends to return to the position “nose forward”:
There is another force with which behavior becomes more interesting. When looking at the plane from the front or from the back, you can see that the wing is slightly bent up by the Latin letter V - this is dictated by concern for the stability of the flight. If the plane does not fly straight ahead, but slightly moves sideways, the lifting forces on the left and right will be different, and it will begin to rotate.
forceLocal.y - lifting force
Something happened against which my sense of beauty protested, but otherwise I would have to greatly complicate the model. Before adding strength, I still try to justify it. If a straight flying plane rotates, for example, with a roll to the left, then the angle of attack of the left wing increases, and the right wing - on the contrary, and this effect slows down the rotation. On other axes, there is probably something similar (in the StateGenerator class, very slight friction during rotation is made for the stability of the computational scheme, and here it is simply so that the plane does not become like a pendulum):
We translate into the global coordinate system:
The system of units - meters, kilograms, seconds. The “fitting coefficients” are given for a reason — I tried to match them to the I-16 parameters. Mass 1400, power 750hp, or (750 * 735.5) Watts.The moment of inertia (estimated) - 5000 along the OX, OY and much less along the OZ (such as the main mass is concentrated in the fuselage of the aircraft, and it is quite long).
Imp5 reported more accurate data: the main moments of inertia are 2440, 5520, 3080 along the “forward”, “up” and “right” axes, respectively.
This physical model does not take into account the rotation of the aircraft, and fall into a spin does not work. In the future, I plan to take several points on each wing and individually for each point calculate the angles of attack and speed of movement relative to the air. The control of the tail and ailerons is implemented as a change in the parameters of the wing pieces. Perhaps then the rotation of the aircraft will honestly fade due to air resistance.
The code that calculates the forces and movement of the aircraft can be replaced at any time with something more serious.
PS I take off my hat to domestic developers of a simulator about airplanes of the Second World War, from which my interest in aviation once began long ago.
It was worth trying to write the physics himself to understand what a titanic work they had done. For example, the rotation of a longitudinally located engine leads to the fact that with a bend in one direction, the nose of the aircraft leads upwards and in the other - downwards. I do not have this effect, like many others. On the one hand, a trifle, but from such trifles a unique behavior for each model is formed.
Warning: further reasoning may well be wrong, my experience is limited to playing simulators and the course of theoretical mechanics. Knowledge of aerodynamics and game physics engines is very scarce. A picture for attracting attention is a photograph, not a screenshot.
“What could be simpler than a plane? The lift is proportional to the square of speed, the engine pulls forward, everything is simple ”- such a thought came to my mind in the summer, and I sat down to write a game. The summer passed, several rakes were collected, and the list of what I planned to add to the project grew very much.
')
In this article I will try to write about the physical component, without being distracted by the graphics and other things. I hope this will help someone - there is not a lot of information on this topic on the Internet.
Examples of working and not-so-code will be on Scala (to understand the essence, the language is not necessary to know). In addition, I use classes for vectors, matrices, and quaternions from libgdx. If someone has questions on the features of their implementation - the engine code is open. For convenience, vectors have been added methods of the type + =, * =, as well as +, -, the benefit in scala is so possible. It also turned out to be convenient to add methods: = for assignment by value.
Code:
v:=(0, 1, 2) v2:=v Equivalent to
v.set(0,1,2) v2.set(v), With the exception of quaternions:
set(x,y,z,w) , but :=(w,x,y,z)The second option is much more familiar to me, and I never use set. I will not lay out my code yet, there is a little trash, which from time to time is being overwritten beyond recognition.
So, let's begin.
Write the class for the model position:
class Position(val pos:Vector3f, val rot:Quaternion) As you might guess, the class will have two final fields describing the position of the center of mass of the model and its orientation relative to the world.
You can get the matrix in libgdx like this:
matrix.idt().translate(position.pos).rotate(position.rot)We introduce a class for the derivative:
class Derivative(val linear: Vector3f, val angular: Vector3f) This class will be convenient not only to describe the speed (first derivative), but also for accelerations and forces.
Stop. How did the quaternion become a vector? Indeed, the angular velocity, acceleration and momentum are usually described by vectors. The essential difference between the angular velocity and orientation is that the orientations are “looped around”, turning two pi is equivalent to “zero” turning. On the contrary, the angular velocity is in principle unlimited.
You can enter the operation of the logarithm of the quaternion q- vector v, which is directed along the direction of the axis of rotation, and its length is equal to the angle of rotation in radians. The exhibitor is the inverse operation. q == exp (v * 2) == exp (v) * exp (v)
The mapping of vector-> quaternion is unambiguous, but the opposite is not. The rotation angle alpha in time dt can correspond to the angular velocity (alpha + 2 * pi * n) / dt, where n is any integer.
Obviously, during the time dt at the angular velocity w, the rotation is q = exp (w * dt).
Code:
class QuaternionS extends Quaternion{ ... def :=(log: Vector3): Unit = { val len = log.len if (len > 0.0001f) { val a = len / 2 val m = Math.sin(a).toFloat / len this :=(Math.cos(a).toFloat, log.x * m, log.y * m, log.z * m) } else this :=(1, 0, 0, 0) } def log(result: Vector3): Unit = { val xyz = Math.sqrt(x*x + y*y + z*z) // if (xyz > 0.0001f) { val m = (Math.acos(w).toFloat * 2f) / xyz result.set(m * x, m * y, m * z) } else { result.set(0f, 0f, 0f) } } } What is an airplane flight from an abstract point of view? By solving a system of second-order differential equations! Set the state of the aircraft:
trait TimeStamp{ def time: Long // } trait VisibleState extends TimeStamp{ def position: Position def speed: Derivative } The class for the real state of the aircraft depends on the features of the model and the degree of its physical development, but it will implement this interface, necessary for drawing the aircraft on the screen.
The task can be divided into two independent parts.
1) the calculation of the forces acting on the aircraft in a given state
2) numerical solution diffura
The first part, as the most interesting, we leave at the end of the article.
For the graphic part, we write the interface:
abstract class Avatar(val player: Player, val model: Airplane) { def getState(time: Long): VisibleState } I will not give the implementation, but the essence is simple - a sequence of states is stored inside with times t1, t2, t3, etc., t (n + 1)> t (n). States are completed if necessary, in the getState method, the next two are interpolated. Thus, it is possible, for example, to read physics 10 times per second and at the same time observe smooth motion at 60 fps.
We write the following interface:
trait Airplane{ def mass: Float def I: Matrix4 private val tempM = new Matrix4() private val tempQ = new QuaternionS() def getInertiaTensor(orientation: Quaternion): Matrix4 = { tempM.set(orientation) tempM.mul(I) tempQ := orientation tempQ.conjugate() tempM.rotate(tempQ) tempM } def getForceAndMoment(state: State, input: PlaneInput): Derivative } The angular momentum is L = I * w, and L and w (angular velocity) are transformed as vectors. Thus, in the transformation L '= qL, w' = qw it turns out:
L '= I' * w '
qL = I '* qw
L = q ^ (- 1) * I '* q * w
We get I = q ^ (- 1) * I '* q, or I' = q * I * q ^ (- 1).
The transformation w '= position.rot * w translates the angular velocity from the local coordinate system to the global one.
The getForceAndMoment method will be discussed later, it calculates the forces and torque acting on the aircraft.
I do not very well imagine how to accurately calculate the movement of a model that moves and rotates with accelerations, so the simplest method was chosen with a fixed step of 20ms.
class StateGenerator { // , // , , private val inversed = new Matrix4() private val omega = new Vector3f() private val iOm = new Vector3f() private val angularAcceleration = new Vector3f() private val inv = new QuaternionS() def nextState(now: State, timeStep: Long, model: Airplane, planeInput: PlaneInput): State = { assert(timeStep > 0) val dt = timeStep * 0.001f // val next = new State(now.time + timeStep) next.position := now.position next.speed := now.speed // val forces = model.getForceAndMoment(now, planeInput) //linear next.speed.linear += forces.linear * (dt / model.mass) next.position.pos += (now.speed.linear + next.speed.linear) * (0.5f * dt) //angular val I = model.getInertiaTensor(now.position.rot) inversed.set(I).inv() omega := now.speed.angular iOm := omega iOm.mul(I) angularAcceleration.setCross(iOm, omega) angularAcceleration += forces.angular angularAcceleration.mul(inversed) next.speed.angular.madd(angularAcceleration, dt) next.speed.angular *= 1f - dt * 0.1f //. . - 10 . val angSp = (next.speed.angular + now.speed.angular) / 2 next.position.rot := new QuaternionS(angSp * dt) * now.position.rot next // Scala return } } More information about the rotation can be found in this article: habrahabr.ru/post/264099 . Honestly, I am not a fan of tensors, I just took a formula from there in vector form in order to obtain angular acceleration. Calculations are made in the coordinate system of the world. By the way, when external forces were disconnected, I was able to observe a movement quite similar to the Janibekov effect.
Forces acting on the aircraft
The plane must be managed:
trait PlaneInput extends TimeStamp { def yaw: Float // right is positive def pitch: Float // up is positive def roll: Float // clockwise is positive def engineTrust: Float } Values ​​are cut to the interval [-1, 1], engine thrust from 0 to 1.
Well, let's move on to the most important part - we will find strength. Here is a little Terra incognita, my knowledge of aerodynamics is very superficial. Errors and inaccuracies are possible.
The first thing that comes to mind is lifting force. In order not to create garbage, on the Internet, a reference book of aviation profiles was found with charts of lift coefficient depending on the angle of attack. The essence turned out to be quite simple - Cy (lift coefficient) grows quite linearly to critical angles, reaches about one, and then the flow from the wing breaks down, and the lift begins to decrease. Also, the graph of the coefficient for the abstract wing can be viewed in the English Wikipedia :
Then I was trapped by a rake - if you read another article on Wikipedia on frontal resistance, you can see that there is some kind of inductive resistance. The surprise is that the lifting force is considered to be in the direction perpendicular to the direction of speed, and not perpendicular to the surface of the wing (as I thought). Since the difference in air pressure above and below the wing still leads to a force perpendicular to the surface of the wing, the projection of this force on the direction opposite to the movement is non-zero.
And no. See comment below. Further text and code leave unchanged.
If we assume that the lift is directed upwards in the reference system of the aircraft, then the inductive force does not seem to be necessary - it has already been taken into account. The orientation of the axes is the same as in openGL:
Ox - right
Oy - up
Oz - back
// def cy(angle: Float) = { val pi = 3.1415928535f if (angle > maxAngle) (pi - angle) / (pi - maxAngle) else if (angle < -maxAngle) -(pi + angle) / (pi - maxAngle) else angle / maxAngle } def liftingForce(speed2: Float, angle: Float, airDensity: Float): Float = cy(angle) * airDensity * speed2 * 0.5f * planesS // . val forceLocal = new Vector3f val speedLocal = state.position.rot^(-1) * state.speed.linear val sp2 = state.speed.linear.length2 val airDensity = WorldParams.atmosphere.aitDensity(state.position.pos) val attackAngle = -Math.asin(speedLocal.y / (speedLocal.length + 0.1f)).toFloat val steeringAngle = Math.asin(speedLocal.x / (speedLocal.length + 0.1f)).toFloat // , ( 0.1 /) , "" . - . forceLocal.y += liftingForce(sp2, attackAngle, airDensity) forceLocal.x -= steeringForce(sp2, steeringAngle, airDensity) forceLocal += speedLocal/(0.1f + speedLocal.length) * -dragForce(sp2, attackAngle, steeringAngle, airDensity) In addition to lift, you will need air resistance force:
dragForce and the force that occurs when the plane flies a little sideways: steeringForce .I do not have sufficient knowledge in aerodynamics. The main goal is simplicity of formulas and, if possible, adequate aircraft behavior for flight angles of attack and glide. 0.5f - the consequences of divider 2 in the formulas. 0.1f - the consequences of adjusting the coefficients.
private def steeringForce(speed2: Float, angle: Float, airDensity: Float): Float = angle * airDensity * speed2 * 0.5f * planeSVert private def dragForce(speed2: Float, attackAngle: Float, steeringAngle: Float, airDensity: Float): Float = { speed2 * (0.5f * airDensity * crossSectionalArea + Math.abs(attackAngle) * 0.1f * planesS + Math.abs(steeringAngle) * 0.1f * planeSVert) } Add traction motor
The model is as simple as possible: no propeller steps, let the engine spend all its power to accelerate the aircraft. No bonuses to the moment of inertia, no torque when changing the number of revolutions. However, the revolutions either. Power = force * speed. To prevent the plane from flying up like a rocket, we will limit the maximum force (using the minimum speed limit).
val speed = Math.max(minSpeed, -speedLocal.z) forceLocal.z -= input.engineTrust * maxPower / speed Control
There is an interesting point - the air accelerated by a screw goes directly to the control planes of the tail, and the plane, in principle, is controlled by the tail a little even on the runway. In addition, the air resistance is trying to twist the plane in the opposite direction of the screw rotation. And up to the heap - the engine has a moment of inertia, when increasing / decreasing the speed of rotation, the plane will also “twist” a little. I will neglect all this ...
As with the wing, a familiar multiplier appears with a square of velocity and air density:
val mul = spLocal.z * spLocal.z * airDensity result.angular.x = (input.pitch + 0.3f) * mul * 2f result.angular.y = -input.yaw * mul * 1f result.angular.z = -input.roll * mul * 5 There is no symmetry in pitch, the aircraft (and the pilot) are much better tolerated by positive overloads than negative ones. In addition, the aircraft itself is stable (if it is not the Su-47) and tends to return to the position “nose forward”:
result.angular.x -= attackAngle * airDensity * mul * 5f result.angular.y -= steeringAngle * airDensity * mul * 5f Do not forget anything?
There is another force with which behavior becomes more interesting. When looking at the plane from the front or from the back, you can see that the wing is slightly bent up by the Latin letter V - this is dictated by concern for the stability of the flight. If the plane does not fly straight ahead, but slightly moves sideways, the lifting forces on the left and right will be different, and it will begin to rotate.
result.angular.z += forceLocal.y * steeringAngle * 1f forceLocal.y - lifting force
Add "friction" to the rotation
Something happened against which my sense of beauty protested, but otherwise I would have to greatly complicate the model. Before adding strength, I still try to justify it. If a straight flying plane rotates, for example, with a roll to the left, then the angle of attack of the left wing increases, and the right wing - on the contrary, and this effect slows down the rotation. On other axes, there is probably something similar (in the StateGenerator class, very slight friction during rotation is made for the stability of the computational scheme, and here it is simply so that the plane does not become like a pendulum):
result.angular.mulAdd(spAngLocal, -airDensity * 5000f) We translate into the global coordinate system:
result.angular.mul(state.position.rot) forceLocal.mul(state.position.rot) result.linear := forceLocal Note
The system of units - meters, kilograms, seconds. The “fitting coefficients” are given for a reason — I tried to match them to the I-16 parameters. Mass 1400, power 750hp, or (750 * 735.5) Watts.
Imp5 reported more accurate data: the main moments of inertia are 2440, 5520, 3080 along the “forward”, “up” and “right” axes, respectively.
This physical model does not take into account the rotation of the aircraft, and fall into a spin does not work. In the future, I plan to take several points on each wing and individually for each point calculate the angles of attack and speed of movement relative to the air. The control of the tail and ailerons is implemented as a change in the parameters of the wing pieces. Perhaps then the rotation of the aircraft will honestly fade due to air resistance.
The code that calculates the forces and movement of the aircraft can be replaced at any time with something more serious.
PS I take off my hat to domestic developers of a simulator about airplanes of the Second World War, from which my interest in aviation once began long ago.
It was worth trying to write the physics himself to understand what a titanic work they had done. For example, the rotation of a longitudinally located engine leads to the fact that with a bend in one direction, the nose of the aircraft leads upwards and in the other - downwards. I do not have this effect, like many others. On the one hand, a trifle, but from such trifles a unique behavior for each model is formed.
Source: https://habr.com/ru/post/266367/
All Articles