Why does this simulation "explode"?

Question

I am attempting to simulate a steerable vehicle (plane, hovercraft, boat) in a fluid (gas or liquid, e.g. air or water).

My JavaScript implementation is based on the 3D rigid body airplane simulator from Physics for Game Developers ported to Three.js's Vector3, Matrix3 and Quaternion classes.

Applying a thrust force to the vehicle causes it to accelerate and move through the medium. It's speed causes drag and lift on a control surface / "rudder", depending on its angle.
The resulting force is offset from the common center of gravity of the vehicle parts ("main" and "rudder") and as such also causes a moment (torque) that causes angular motion of the vehicle (steering).

Drag and lift are quadratic to the speed.
Thus drag will limit the maximum speed given a thrust force.

In steady state (no acceleration) the sum of all forces is zero, and that means, the sum of drag and lift are equal and opposing the thrust.

This is where the program can start to oscillate (either stable or it "explodes"), depending on the parameters.

Oscillation seems to happen when the control surface's normal vector and the drag vector become orthogonal, as then the dot product becomes zero, which is used to determine the angle of attack (control surface in medium), which then is also approaching zero.

This happens later rather than sooner with:

smaller RHO (affects both, lift and drag)
smaller lift coefficient lc
larger drag coefficient dc (larger drag means less acceleration and speed)

This JSFiddle shows the problem with a top view (Z is up, vehicle moves along X on XY-plane) on a vector representation of some of the forces at play:

https://jsfiddle.net/ghcse56f/ (easier to play with than snippets here on SO)

It includes a crude plot of some values (different scales). Calculation of lift and drag start around line 215.

Why is this problem happening and how can it be fixed?

I've updated the demo to use a lift coefficient that is proportional to the angle of attack, reduced the drag coefficient and increased RHO - the simulation falls apart in step 227.

//import * as THREE from "three";
//import { Euler, Matrix3, Vector3 } from "three";

// colors
const cred = 0xff0000;
const cgreen = 0x00ff00;
const cblue = 0x0000ff;
const cyellow = 0xffff00;
const ccyan = 0x00ffff;
const cmagenta = 0xff00ff;
const cblack = 0x000000;
const cwhite = 0xffffff;

/** degree to radian */
function rad(deg) {
  return (deg / 180) * Math.PI;
}
/** radian to degree */
function deg(rad) {
  return (rad / Math.PI) * 180;
}

// string of vector components
function vstr(v) {
  return `${v.x}, ${v.y}, ${v.z}, `
}


class Vehicle {
  constructor(scene) {

    this.rudder = 10; // [deg]
    this.thrust = 1; // [N]

    // this.length = 0.3;
    // this.width = 0.15;

    this.fMass = 0.0; // total mass
    this.mInertia = new THREE.Matrix3(); // mass moment of inertia in body coordinates
    this.mInertiaInverse = new THREE.Matrix3(); // inverse of mass moment of inertia
    this.vPosition = new THREE.Vector3(0, 0, 0.001); // position in earth coordinates
    this.vVelocity = new THREE.Vector3(); // velocity in earth coordinates
    this.vVelocityBody = new THREE.Vector3(.0, 0, 0); // velocity in body coordinates
    this.vAngularVelocity = new THREE.Vector3(); // angular velocity in body coordinates
    this.vEulerAngles = new THREE.Vector3(0, 0, 0); // Euler angles in body coordinates
    this.qOrientation = new THREE.Quaternion(); // orientation in earth coordinates
    this.qOrientation.setFromEuler(new THREE.Euler().setFromVector3(this.vEulerAngles));
    this.vForces = new THREE.Vector3(); // total force on body
    this.vMoments = new THREE.Vector3(); // total moment (torque) on body

    this.elements = [{
        type: "main",
        fMass: 1.0,
        vDCoords: new THREE.Vector3(0, 0.0, 0.0),
        vCGCoords: null,
        vLocalInertia: new THREE.Vector3(.001, .001, .001), // TODO: 
        rotates: false,
        fIncidence: 0.0,
        fDihedral: 90.0,
        fArea: 0.5,
      },
      {
        type: "rudder",
        fMass: .01,
        vDCoords: new THREE.Vector3(-.1, 0.0, 0.0),
        vCGCoords: null,
        vLocalInertia: new THREE.Vector3(.001, .001, .001), // TODO: 
        rotates: true,
        fIncidence: this.rudder,
        fDihedral: 90.0,
        fArea: 0.1,
      },
    ];

    // Calculate the vector normal (perpendicular) to each lifting surface.
    // This is required when you are calculating the relative air velocity for lift and drag calculations.
    for (let e of this.elements) {
      let inc = rad(e.fIncidence);
      let dih = rad(e.fDihedral);
      e.vNormal = new THREE.Vector3(
        Math.sin(inc), // 
        Math.cos(inc) * Math.sin(dih), // 
        Math.cos(inc) * Math.cos(dih) // 
      );
      e.vNormal.normalize();
    }

    // Calculate total mass
    for (let e of this.elements) {
      this.fMass += e.fMass;
    }

    // Calculate combined center of gravity location
    let vMoment = new THREE.Vector3();
    for (let e of this.elements) {
      vMoment.add(e.vDCoords.clone().multiplyScalar(e.fMass))
    }
    let vCG = vMoment.divideScalar(this.fMass);

    // Calculate coordinates of each element with respect to the combined CG
    for (let e of this.elements) {
      e.vCGCoords = e.vDCoords.clone().sub(vCG);
    }

    // (This inertia matrix (tensor) is in body coordinates)
    // Here we’re using a vector to represent the three local moment of inertia terms 
    // and we’re also assuming that the local products of inertia are zero for each element.
    let Ixx = 0,
      Iyy = 0,
      Izz = 0,
      Ixy = 0,
      Ixz = 0,
      Iyz = 0;
    for (let e of this.elements) {
      Ixx += e.vLocalInertia.x + e.fMass * (e.vCGCoords.y ** 2 + e.vCGCoords.z ** 2);
      Iyy += e.vLocalInertia.y + e.fMass * (e.vCGCoords.x ** 2 + e.vCGCoords.z ** 2);
      Izz += e.vLocalInertia.z + e.fMass * (e.vCGCoords.x ** 2 + e.vCGCoords.y ** 2);
      Ixy += e.fMass * (e.vCGCoords.x * e.vCGCoords.y);
      Ixz += e.fMass * (e.vCGCoords.x * e.vCGCoords.z);
      Iyz += e.fMass * (e.vCGCoords.y * e.vCGCoords.z);
    }
    // mass, inertia matrix and the inverse
    this.mInertia.e11 = Ixx;
    this.mInertia.e12 = -Ixy;
    this.mInertia.e13 = -Ixz;

    this.mInertia.e21 = -Ixy;
    this.mInertia.e22 = Iyy;
    this.mInertia.e23 = -Iyz;

    this.mInertia.e31 = -Ixz;
    this.mInertia.e32 = -Iyz;
    this.mInertia.e33 = Izz;

    this.mInertiaInverse = this.mInertia.clone().invert();


    // three.js
    this.object = new THREE.Group();

    const axes = new THREE.AxesHelper(0.1);
    axes.position.copy(vCG);
    axes.setColors(cblack, new THREE.Color('rgb(50%,50%,50%)'), new THREE.Color('rgb(50%,50%,50%)'));
    this.object.add(axes);

    this.Athrust = new THREE.ArrowHelper(new THREE.Vector3(1, 0, 0), vCG, 1, cgreen);
    this.object.add(this.Athrust);

    this.Anormal = new THREE.ArrowHelper(this.elements[1].vNormal.clone(), this.elements[1].vCGCoords, 1, ccyan);
    this.object.add(this.Anormal);

    this.Alift = new THREE.ArrowHelper(new THREE.Vector3(0, 1, 0), this.elements[1].vCGCoords, 1, cblue);
    this.object.add(this.Alift);
    this.Adrag = new THREE.ArrowHelper(new THREE.Vector3(-1, 0, 0), this.elements[1].vCGCoords, 1, cred);
    this.object.add(this.Adrag);
    this.Aresultant = new THREE.ArrowHelper(new THREE.Vector3(-1, 0, 0), this.elements[1].vCGCoords, 1, cyellow);
    this.object.add(this.Aresultant);
    this.Aforce = new THREE.ArrowHelper(new THREE.Vector3(1, 0, 0), vCG, 1, cwhite);
    this.object.add(this.Aforce);
    this.Amoment = new THREE.ArrowHelper(new THREE.Vector3(0, 0, 1), vCG, 1, cmagenta);
    this.object.add(this.Amoment);

    // add self to scene
    scene.add(this.object);

  } // ctor


  physics(dt) {

    // reset forces and moments:
    let vFb = new THREE.Vector3();
    let vMb = new THREE.Vector3();
    this.vForces.set(0, 0, 0);
    this.vMoments.set(0, 0, 0);

    // thrust vector acts through the CG
    let Thrust = new THREE.Vector3(1, 0, 0).multiplyScalar(this.thrust);

    // forces and moments in body space:
    let vLocalVelocity = new THREE.Vector3();
    let fLocalSpeed = 0;
    let vDragVector = new THREE.Vector3();
    let fAttackAngle = 0;
    let tmp = 0.0;
    let vResultant = new THREE.Vector3();

    let vtmp = new THREE.Vector3();

    // "rudder" forces only
    for (let e of this.elements) {
      if (e.type == "main")
        continue;

      // Calculate local velocity at element
      // The local velocity includes the velocity due to linear motion of the airplane,
      // plus the velocity at each element due to the rotation of the airplane.
      // Here's the rotational part
      vtmp.crossVectors(this.vAngularVelocity, e.vCGCoords);
      vLocalVelocity.addVectors(this.vVelocityBody, vtmp);

      // Calculate local air speed
      fLocalSpeed = vLocalVelocity.length();
      pixel(fLocalSpeed * 150, 'green');

      // Find the direction in which drag will act.
      // Drag always acts inline with the relative velocity but in the opposing direction
      if (fLocalSpeed > 0) {
        vDragVector = vLocalVelocity.clone().divideScalar(-fLocalSpeed);
      }

      // Find the direction in which lift will act.
      // Lift is always perpendicular to the drag vector
      const vDN = new THREE.Vector3().crossVectors(vDragVector, e.vNormal)
      let vLiftVector = new THREE.Vector3().crossVectors(vDN, vDragVector);
      vLiftVector.normalize();
      console.log("vliftVector", vstr(vLiftVector));

      // Find the angle of attack.
      // The attack angle is the angle between the lift vector and the element normal vector.
      // Note, the sine of the attack angle is equal to the cosine of the angle between the drag vector and the normal vector.
      tmp = vDragVector.dot(e.vNormal);
      console.log("vDragVector", vstr(vDragVector));
      console.log("e.vNormal", vstr(e.vNormal));
      console.log("tmp", tmp);

      if (tmp > 1) tmp = 1;
      else if (tmp < -1) tmp = -1;
      fAttackAngle = deg(Math.asin(tmp));
      console.log("fAttackAngle", fAttackAngle);
      pixel(fAttackAngle * 50, 'magenta');


      // Determine the resultant force (lift and drag) on the element.

      // fluid density
      // const RHO = 1.225;       // kg/m^3 for AIR
      const RHO = 100; // kg/m^3 for WATER
      tmp = 0.5 * RHO * fLocalSpeed ** 2 * e.fArea;

      // lift coefficient
      // sign fixes vector direction
      // depends on angle of incidence - constant to simplify things
      let lc = 0;
      //if (fAttackAngle < 0)
      //  lc = -.1;
      //if (fAttackAngle > 0)
      //  lc = .1;
      lc = fAttackAngle / 2;

      // drag coefficient
      // depends on angle of incidence - constant to simplify things
      const dc = .01;

      let vLift, vDrag;
      vLift = vLiftVector.multiplyScalar(lc);
      vDrag = vDragVector.multiplyScalar(dc);

      this.Alift.setDirection(vLift.clone().normalize());
      this.Alift.setLength(vLift.length() * tmp);
      this.Adrag.setDirection(vDrag.clone().normalize());
      this.Adrag.setLength(vDrag.length() * tmp);

      pixel(vLift.length() * tmp * 150, 'blue');
      pixel(vDrag.length() * tmp * 150, 'red');

      const vLiftDrag = vDrag.add(vLift);
      vResultant = vLiftDrag.multiplyScalar(tmp);

      this.Aresultant.setDirection(vResultant.clone().normalize());
      this.Aresultant.setLength(vResultant.length());

      // Keep a running total of these resultant forces (total force)
      vFb.add(vResultant);

      // Calculate the moment about the CG of this element's force
      // and keep a running total of these moments (total moment)
      vtmp.crossVectors(e.vCGCoords, vResultant);
      vMb.add(vtmp);

    } // elements

    // force
    pixel(vFb.length() * 150, 'yellow'); // drag and lift

    // Now add the thrust
    vFb.add(Thrust);

    pixel(vFb.length() * 150, 'white'); // total
    this.Aforce.setDirection(vFb.clone().normalize());
    this.Aforce.setLength(vFb.length());

    // stop if sim "explodes"
    if (isNaN(vFb.length()))
      return -1;

    // Convert forces from model space to earth space
    this.vForces = vFb.clone().applyQuaternion(this.qOrientation);


    // moments
    this.vMoments.add(vMb);

    this.Amoment.setDirection(vMb.clone().normalize());
    this.Amoment.setLength(vMb.length());


    // INTEGRATOR (EULER)

    // calculate the acceleration in earth space:
    let vAe = this.vForces.clone().divideScalar(this.fMass);

    // calculate the velocity in earth space:
    vAe.multiplyScalar(dt);
    this.vVelocity.add(vAe);

    // calculate the position in earth space:
    let vdp = this.vVelocity.clone().multiplyScalar(dt)
    this.vPosition.add(vdp);


    // calculate the angular velocity in body space:
    let mv = this.vAngularVelocity.clone().applyMatrix3(this.mInertia);
    let cp = new THREE.Vector3().crossVectors(this.vAngularVelocity, mv);
    let sub = this.vMoments.clone().sub(cp);
    sub.multiplyScalar(dt);
    this.vAngularVelocity.add(sub.applyMatrix3(this.mInertiaInverse));

    // calculate the new rotation quaternion:
    // Quaternion * Vector
    let q = this.qOrientation.clone();
    let v = this.vAngularVelocity.clone();
    let vq = new THREE.Quaternion(v.x, v.y, v.z, 0);
    let qv = q.multiply(vq);

    // Q * Scalar
    const scalar = 0.5 * dt;
    qv.x *= scalar;
    qv.y *= scalar;
    qv.z *= scalar;
    qv.w *= scalar;

    // Q+Q
    this.qOrientation.x += qv.x;
    this.qOrientation.y += qv.y;
    this.qOrientation.z += qv.z;
    this.qOrientation.w += qv.w;

    this.qOrientation.normalize();

    // calculate the velocity in body space:
    // (we'll need this to calculate lift and drag forces)
    this.vVelocityBody = this.vVelocity.clone().applyQuaternion(this.qOrientation.clone().conjugate());

    // get the Euler angles for our information
    let u = new THREE.Vector3().setFromEuler(new THREE.Euler().setFromQuaternion(this.qOrientation));

    this.vEulerAngles.x = u.x; // roll
    this.vEulerAngles.y = u.y; // pitch
    this.vEulerAngles.z = u.z; // yaw

    pixel(deg(u.z), 'orange');
    pixel(this.rudder, 'black');

    // threedee
    this.object.position.copy(this.vPosition);
    this.object.rotation.z = this.vEulerAngles.z;
    camera.position.set(this.vPosition.x, this.vPosition.y - 0.0000001, 6);
    camera.lookAt(this.vPosition.x, this.vPosition.y, 0);

    return 0;
  } // physics


} // Vehicle



// minimal graphing
let plot = document.querySelector("#plot");
plot.width = 500;
plot.height = 300;
let ctx = plot.getContext('2d');
ctx.fillStyle = '#333';
ctx.fillRect(0, 0, plot.width, plot.height);

ctx.strokeStyle = '#555';
ctx.beginPath();
ctx.moveTo(0, plot.height / 2);
ctx.lineTo(plot.width, plot.height / 2);
ctx.stroke();

let n = 0;

function pixel(y, color) {
  ctx.fillStyle = color;
  ctx.fillRect(n, plot.height / 2 - y, 1, 1);
}


// three.js
let renderer = new THREE.WebGLRenderer({
  antialias: true
});
renderer.setPixelRatio(window.devicePixelRatio);
renderer.setSize(window.innerWidth / 2, window.innerHeight / 2);
document.body.appendChild(renderer.domElement);

let scene = new THREE.Scene();
scene.background = new THREE.Color('rgb(25%,25%,25%)');

const axes = new THREE.AxesHelper();
scene.add(axes);

const grid = new THREE.GridHelper(20, 20);
grid.rotation.x = Math.PI / 2; // Z is up 
scene.add(grid);

let camera = new THREE.PerspectiveCamera(50, window.innerWidth / window.innerHeight, 0.1, 50);
camera.up.set(0, 0, 1); // Z is up
camera.position.set(0, -0.0000001, 6);

// vehicle
let craft = new Vehicle(scene);

// simulation
let timestart = performance.now();
const dt = 0.01;
const duration = plot.width / 60 * 1000;

function step(time) {
  if (time - timestart > duration) return;

  console.log(n);
  n++;

  const r = craft.physics(dt);
  if (r == -1) {
    console.log("ERROR");
    return;
  }

  renderer.render(scene, camera);
  requestAnimationFrame(step);
}
requestAnimationFrame(step);

window.addEventListener("resize", e => {
  camera.aspect = window.innerWidth / window.innerHeight;
  camera.updateProjectionMatrix();
  renderer.setSize(window.innerWidth, window.innerHeight);
});

body {
  margin: 0px;
}

Plot:
<span style="color:red">red: drag</span>,
<span style="color:green">green: localSpeed</span>,
<span style="color:blue">blue: lift</span>,
<span style="color:gold">yellow: drag+lift force</span>,
<span style="color:magenta">magenta: attackangle</span>,
<span style="color:gray">white: total</span>,
<span style="color:black">black: rudder angle</span>,
<span style="color:orange">orange: vehicle angle</span>
<br>
<canvas id="plot"></canvas>
<br> Vectors:
<span style="color:red">red: drag</span>,
<span style="color:green">green: thrust</span>,
<span style="color:blue">blue: lift</span>,
<span style="color:gold">yellow: drag+lift force</span>,
<span style="color:cyan">cyan: rudder normal</span>,
<span style="color:gray">white: total</span>,

<script src="https://cdnjs.cloudflare.com/ajax/libs/three.js/0.150.0/three.min.js"></script>

The coefficient needs to be proportional to attackAngle (currently either +/-0.1 or 0). — handle, Feb 23 '23 at 16:18
You're outputting the step-by-step data, which is good. It may be worthwhile to collect and export it into a CSV or similar, so you can analyze it in Excel/Sheets. Sometimes it's easier to spot problems if you can observe several previous steps at the same time. — TheJim01, Feb 23 '23 at 19:06
@TheJim01 Yes, plotting the data helps a lot (run the snippet fullscreen or use the JSFiddle to see it). — handle, Feb 24 '23 at 08:53

score 0 · Answer 1 · answered Feb 23 '23 at 21:09

0

I'm not sure how you're calculating the vtmp vector, but I'm getting insanely large values right before crashing:

x: 0,
y: -6e+215,
z: -3e+199

When calculating the .length() of such a huge vector, you're surpassing JavaScript's Number.MAX_VALUE limit, which is 1.7 with 308 zeroes behind it. Instead, it returns Infinity and you're getting fLocalSpeed = Infinity, which immediately causes divisions by 0 and all kinds of other problems with your dot products and crossVectors.

Like I said, I'm not sure how you're calculating vtmp, but with your approach it seems like vAngularVelocity reaches incredibly large values... maybe it's accurate that an airplane would explode when subjected to such large forces! :D

Maybe you need to protect against divisions by 0 or some other undesired result stemming from cross vectors or dot products?

answered Feb 23 '23 at 21:09

M -

26,908
11
49
81

Thanks for having a look! `vtmp` is the cross product of the combined drag and lift force vector and the distance of the "rudder" from the combined center of gravity of the vehicle, i.e. it's a torque (moment). Yeah, that would explode as well with the oscillation of attackangle. It's hard to tell where/when the issue starts just by staring at the numbers... I think the root cause may be the integrator "overshooting" zero, leading to a change of sign and thus the runaway oscillation (ending in Infinity or NaN), depending on the gain (basically the whole lift/drag calculation). – handle Feb 24 '23 at 08:37
In support of that, reducing `dt` increases the progress before the simulation fails considerably. – handle Feb 24 '23 at 08:37

Why does this simulation "explode"?

1 Answers1