Initial implementation of a new GroundReactions class which ultimatelly handle surface displacements for the likes of rocky terrain or ocean waves across FDM's.

This first implementation handles static and roling friction coefficients of tyres on dry and wet asphalt, snowy frozen lakes and frozen lakes with a water layer.

Also handles the bumpiness of different landclasses.

TODO:
 - Implement the same functionality for YASIM.
 - Add ocean wave hydrodynamics.
 - handle dry, moist and wet mud.

src/FDM/CMakeLists.txt

@@ -140,6 +140,7 @@ set(SOURCES
 	flightProperties.cxx
 	TankProperties.cxx
 	groundcache.cxx
+        groundreactions.cxx
 	${SP_FDM_SOURCES}
 	ExternalNet/ExternalNet.cxx
 	ExternalPipe/ExternalPipe.cxx
@@ -158,6 +159,7 @@ set(HEADERS
 	flight.hxx
 	flightProperties.hxx
 	groundcache.hxx
+        groundreactions.hxx
 	${SP_FDM_HEADERS}
 	)
 

src/FDM/JSBSim/JSBSim.cxx

@@ -39,6 +39,7 @@
 #include <simgear/bvh/BVHMaterial.hxx>
 
 #include <FDM/flight.hxx>
+#include <FDM/groundreactions.hxx>
 
 #include <Aircraft/controls.hxx>
 #include <Main/globals.hxx>
@@ -1403,39 +1404,15 @@ FGJSBsim::get_agl_ft(double t, const FGColumnVector3& loc, double alt_off,
   if (terrain_active)
   {
     static bool material_valid = false;
-    if (material) {
-      // Traction Coefficient factors for Tires vs Dry Asphalt:
-      // Dry Asphalt   1.0             Wet Asphalt     0.75
-      // Dry Ice/Snow  0.375           Wet Ice         0.125
-      bool rain = (precipitation->getDoubleValue() > 0.0);
-      bool is_lake = (*material).solid_is_prop();
-      bool solid = (*material).get_solid();
-      double friction_fact, pressure;
-      if (is_lake && solid) { // on ice
-        GroundReactions->SetBumpiness(0.1);
-        if (rain) {
-          if (temperature->getDoubleValue() > 0.0) { // Wet Ice
-            GroundReactions->SetRollingFFactor(0.05);
-            friction_fact = 0.125*(*material).get_friction_factor();
-          } else { // Snow
-            GroundReactions->SetRollingFFactor(0.15);
-            friction_fact = 0.375*(*material).get_friction_factor();
-          }
-        } else { // Dry Ice
-          GroundReactions->SetRollingFFactor(0.05);
-          friction_fact = 0.375*(*material).get_friction_factor();
-        }
-        pressure = (*material).get_load_resistance()*1000;
-      } else { // not on ice
-        GroundReactions->SetRollingFFactor((*material).get_rolling_friction()/0.02);
-        GroundReactions->SetBumpiness((*material).get_bumpiness());
-        friction_fact = (*material).get_friction_factor();
-        if (rain) friction_fact *= 0.75;
-        pressure = (*material).get_load_resistance();
-      }
-      GroundReactions->SetStaticFFactor(friction_fact);
-      GroundReactions->SetMaximumForce(pressure*0.00014503773800721815);
-      GroundReactions->SetSolid(solid);
+
+    if (updateGroundReactions(material)) {
+      setPosition(pt);
+      setHeading(Propagate->GetEuler(FGJSBBase::ePsi));
+      GroundReactions->SetBumpHeight(getBumpHeight());
+      GroundReactions->SetRollingFFactor(getRolingFrictionFactor());
+      GroundReactions->SetStaticFFactor(getStaticFrictionFactor());
+      GroundReactions->SetMaximumForce(getPressure()*0.00014503773800721815);
+      GroundReactions->SetSolid(getSolid());
       GroundReactions->SetPosition(pt);
       material_valid = true;
     } else {

src/FDM/JSBSim/models/FGSurface.cpp

@@ -52,7 +52,6 @@ CLASS IMPLEMENTATION
 FGSurface::FGSurface(FGFDMExec* fdmex, int number) :
    contactNumber(number)
 {
-  eSurfaceType = ctBOGEY;
   _PropertyManager = fdmex->GetPropertyManager();
   resetValues();
 }
@@ -116,6 +115,12 @@ void FGSurface::bind(void)
 
 float FGSurface::GetBumpHeight()
 {
+  // bump hight might have been manually set
+  if (bumpHeight < DBL_MAX) {
+    float rv = bumpHeight;
+    bumpHeight = DBL_MAX;
+    return rv;
+  }
   if (bumpiness < 0.001) return 0.0f;
 
   double x = pos[0]*0.1;

src/FDM/JSBSim/models/FGSurface.h

@@ -95,6 +95,9 @@ public:
       pos[0] = pt[0]; pos[1] = pt[1]; pos[2] = pt[2];
   }
 
+  /// Sets the height of the bump at the provided offset
+  void SetBumpHeight(float height) { bumpHeight = height; }
+
 
   /// Gets the static friction factor of the surface area
   double GetStaticFFactor(void) { return staticFFactor; }
@@ -112,25 +115,28 @@ public:
   bool GetSolid(void) { return isSolid; }
 
   /// Returns the height of the bump at the provided offset
-  float  GetBumpHeight();
+  float GetBumpHeight();
 
   std::string GetSurfaceStrings(std::string delimeter) const;
   std::string GetSurfaceValues(std::string delimeter) const;
 
 protected:
-  ContactType eSurfaceType;
-  double staticFFactor, rollingFFactor;
-  double maximumForce;
-  double bumpiness;
-  bool isSolid;
+  ContactType eSurfaceType = ctBOGEY;
+  double staticFFactor = 1.0;
+  double rollingFFactor = 1.0;
+  double maximumForce = DBL_MAX;
+  double bumpiness = 0.0;
+  bool isSolid = true;
 
-  double staticFCoeff, dynamicFCoeff;
+  double staticFCoeff = 1.0;
+  double dynamicFCoeff = 1.0;
 
 private:
-  int contactNumber;
-  double pos[3];
+  int contactNumber = 0;
+  double bumpHeight = DBL_MAX;
+  double pos[3] = { 0.0, 0.0, 0.0 };
 
-  FGPropertyManager* _PropertyManager;
+  FGPropertyManager* _PropertyManager = NULL;
 
   static std::string _CreateIndexedPropertyName(const std::string& Property, int index);
 };

src/FDM/flight.hxx

@@ -81,6 +81,7 @@
 #include <simgear/structure/subsystem_mgr.hxx>
 #include <simgear/props/tiedpropertylist.hxx>
 #include <FDM/groundcache.hxx>
+#include <FDM/groundreactions.hxx>
 #include <FDM/AIWake/AIWakeGroup.hxx>
 
 namespace simgear {
@@ -213,6 +214,7 @@ class FGInterface : public SGSubsystem
     // the ground cache object itself.
     FGGroundCache ground_cache;
 
+    GroundReactions _groundReactions;
     AIWakeGroup wake_group;
 
     void set_A_X_pilot(double x)
@@ -762,6 +764,26 @@ public:
     // the wire end position.
     void release_wire(void);
 
+    // Ground reactions
+    bool updateGroundReactions(simgear::BVHMaterial const*& material) {
+        return _groundReactions.update(material);
+    }
+    inline void setHeading(float h) { _groundReactions.setHeading(h); }
+    inline void setPosition(double pt[3]) {
+        _groundReactions.setPosition(pt);
+    }
+    inline float getPressure() { return _groundReactions.getPressure(); }
+    inline float getBumpiness() { return _groundReactions.getBumpiness(); }
+    inline float getBumpHeight() { return _groundReactions.getBumpHeight();}
+    inline float getStaticFrictionFactor() {
+        return _groundReactions.getStaticFrictionFactor();
+    }
+    inline float getRolingFrictionFactor() {
+        return _groundReactions.getRolingFrictionFactor();
+    }
+    inline bool getSolid() { return _groundReactions.getSolid(); }
+
+
     // Manages the AI wake computations.
     void add_ai_wake(FGAIAircraft* ai) { wake_group.AddAI(ai); }
     void reset_wake_group(void) { wake_group.gc(); }

src/FDM/groundreactions.cxx

@@ -0,0 +1,118 @@
+/*
+ * groundreactions.hxx
+ * caclulate ground reactions based on ground properties and weather scenarios
+ *
+ * SPDX-FileCopyrightText: (C) 2023 Erik Hofman <erik@ehof,am.com>
+ * SPDX-License-Identifier: GPL-2.0-or-later
+ */
+
+
+#include <cmath>
+
+#include <simgear/bvh/BVHMaterial.hxx>
+#include <simgear/math/sg_random.hxx>
+#include <Main/fg_props.hxx>
+
+#include "groundreactions.hxx"
+
+GroundReactions::GroundReactions()
+{
+    precipitation = fgGetNode("/environment/rain-norm", true);
+    temperature = fgGetNode("/environment/temperature-degc",true);
+    dew_point = fgGetNode("/environment/dewpoint-degc", true);
+    ground_wind = fgGetNode("/environment/config/boundary/entry[0]/wind-speed-kt",true);
+    turbulence_gain = fgGetNode("/environment/turbulence/magnitude-norm",true);
+    turbulence_rate = fgGetNode("/environment/turbulence/rate-hz",true);
+    turbulence_model = fgGetNode("/environment/params/jsbsim-turbulence-model",true);
+
+    wind_from_north= fgGetNode("/environment/wind-from-north-fps",true);
+    wind_from_east = fgGetNode("/environment/wind-from-east-fps" ,true);
+    wind_from_down = fgGetNode("/environment/wind-from-down-fps" ,true);
+}
+
+bool
+GroundReactions::update(simgear::BVHMaterial const*& material)
+{
+    if (material)
+    {
+        float friction_factor = (*material).get_friction_factor();
+
+        // TODO: Moist Mud, Wed Mud, Frozen Mud
+        // Traction coefficient factors for tires vs. dry asphalt:
+        // Dry Asphalt   1.0             Wet Asphalt     0.75
+        // Dry Ice/Snow  0.375           Wet Ice         0.125
+        wet = (precipitation->getDoubleValue() > 0.0);
+        snow = (wet && temperature->getDoubleValue() <= 0.0);
+        water_body = (*material).solid_is_prop();
+        solid = (*material).get_solid();
+        frozen = (water_body && solid);
+        if (frozen) // on ice
+        {
+            bumpiness = 0.1f;
+            if (wet)
+            {
+                if (!snow) // Wet Ice
+                {
+                    roling_friction_factor = 0.05f;
+                    static_friction_factor = 0.125f*friction_factor;
+                }
+                else // snow
+                {
+                    roling_friction_factor = 0.15f;
+                    static_friction_factor = 0.375f*friction_factor;
+                }
+            }
+            else // Dry Ice
+            {
+                roling_friction_factor = 0.05f;
+                static_friction_factor = 0.375f*friction_factor;
+            }
+            pressure = (*material).get_load_resistance()*1000;
+        }
+        else // not on ice
+        {
+            roling_friction_factor = (*material).get_rolling_friction()/0.02f;
+            bumpiness = (*material).get_bumpiness();
+            static_friction_factor = friction_factor;
+            if (wet) static_friction_factor *= 0.75f;
+            pressure = (*material).get_load_resistance();
+        }
+        return true;
+    }
+    return false;
+}
+
+void
+GroundReactions::setPosition(const double pt[3])
+{
+    pos[0] = pt[0];
+    pos[1] = pt[1];
+    pos[2] = pt[2];
+}
+
+float
+GroundReactions::getBumpHeight()
+{
+    static const float maxGroundBumpAmplitude = 0.4;
+
+    if (bumpiness < 0.001f) return 0.0f;
+
+    float s = sinf(heading);
+    float c = cosf(heading);
+    float px = pos[0]*0.1f;
+    float py = pos[1]*0.1f;
+    px -= floor(px);
+    py -= floor(py);
+    px *= 1000.0f;
+    py *= 1000.0f;
+
+    int x = c*px - s*py;
+    int y = s*px + c*py;
+
+    // TODO: differentiate between land/water and
+    //       implement hydrodynamcs
+
+//  pc_init(x+y);
+    return pc_map_rand(x, y, 1) * bumpiness*maxGroundBumpAmplitude;
+}
+

src/FDM/groundreactions.hxx

@@ -0,0 +1,71 @@
+/*
+ * groundreactions.hxx
+ * caclulate ground reactions based on ground properties and weather scenarios
+ *
+ * SPDX-FileCopyrightText: (C) 2023 Erik Hofman <erik@ehof,am.com>
+ * SPDX-License-Identifier: GPL-2.0-or-later
+ */
+
+
+#pragma once
+
+# include <float.h>
+
+#include <simgear/math/SGMath.hxx>
+
+class GroundReactions
+{
+public:
+    GroundReactions();
+    ~GroundReactions() = default;
+
+    // call once before processing all gear locations
+    bool update(simgear::BVHMaterial const*& material);
+
+    // first set the gear position in local frame
+    void setPosition(const double pt[3]);
+
+    // and the true heading in radians
+    inline void setHeading(float h) { heading = h; };
+
+    // then get the bum height at that position
+    float getBumpHeight();
+
+
+   inline float getPressure() { return pressure; }
+   inline float getBumpiness() { return bumpiness; }
+   inline float getStaticFrictionFactor() { return static_friction_factor; }
+   inline float getRolingFrictionFactor() { return roling_friction_factor; }
+
+   inline bool getSolid() { return solid; }
+
+private:
+    SGPropertyNode_ptr precipitation;
+    SGPropertyNode_ptr temperature;
+    SGPropertyNode_ptr dew_point;
+    SGPropertyNode_ptr ground_wind;
+    SGPropertyNode_ptr turbulence_gain;
+    SGPropertyNode_ptr turbulence_rate;
+    SGPropertyNode_ptr turbulence_model;
+
+    SGPropertyNode_ptr wind_from_north;
+    SGPropertyNode_ptr wind_from_east;
+    SGPropertyNode_ptr wind_from_down;
+
+    SGVec3d pos;
+    float heading;
+
+    float random = 0.0f;
+    float bumpiness = 0.0f;
+    float static_friction_factor = 1.0f;
+    float roling_friction_factor = 1.0f;
+    float pressure = FLT_MAX;
+
+    bool wet = false;		// is it wet?
+    bool snow = false;		// is it snowed over?
+    bool solid = true;		// is it solid ground or water?
+    bool water_body = false;	// if water, can it ge frozen over?
+    bool frozen = false;	// is it actuall frozen?
+    bool mud = false;		// is it mud?
+};
+