CFD Analysis of MARID

 

CFD Analysis of MARID: Extracting Aerodynamic Stability Derivatives with SimScale

 

MARID is an autonomous long-endurance drone built around a single aft-mounted electric thruster, rotating folding wings, and a liquid hydrogen fuel system. Unlike a conventional aircraft, it has no horizontal tail — pitch stability must come from the wing geometry itself, or from active control. Before flying fast, we needed numbers: specifically, how lift scales with angle of attack (CL_α), and whether the aircraft self-corrects pitch disturbances or amplifies them (Cm_α).

To find out, I ran three RANS CFD simulations on SimScale at α = 0°, 5°, and 15° at V = 60 m/s.

Geometry and Mesh

The simulation geometry is the MARID-alt configuration (B9_TE5.stl) — 77,634 triangles including the V-belly keel fillet. The domain is a large rectangular box with the drone nose pointing into the inlet face. SimScale generated a snappyHexMesh with 6.4 million cells using hex-dominant internal filling and automatic surface layers.

Figure 1 — SimScale mesh domain. The blue bounding box is the outer flow domain; the drone sits at the center with the nose facing the inlet (left face).

 

 

Solver settings:

• Turbulence model: incompressible k-ω SST
• Algorithm: SIMPLE, 1,000 iterations per run
• Reference area S = 1.59 m², chord c = 0.35 m
• ρ = 1.225 kg/m³, V = 60 m/s → q·S = 3,503.6 N, q·S·c = 1,226.3 N·m

Angle of attack was encoded directly in the inlet velocity: U_x = −V cos(α),  U_z = +V sin(α) (nose-up positive). This avoids rotating the geometry between runs.

Flow Visualization

Figure 2 — Surface pressure distribution at low AoA. Pressure loading is spread broadly across the wing planform; the leading-edge suction peak is visible along the wing tips.

 

 

At low AoA the surface pressure is remarkably uniform across the upper wing surface, with the expected low-pressure band along the leading edges. The keel fillet channels pressure cleanly into the belly with no obvious separation.

Figure 3 — Turbulent kinematic viscosity (ν_t) in a longitudinal plane. Blue indicates near-laminar flow at the wing surfaces; the red/orange wake region behind the aft body extends several body lengths downstream.

 

 

The ν_t slice is the most physically revealing image. A thin laminar boundary layer (blue) clings to both wing surfaces, transitioning into a pronounced turbulent wake behind the aft body. The high-ν_t region (red/orange) marks where turbulent mixing dominates. This wake is structurally significant: the aft body is not generating the tail-like restoring pitch moment a conventional aircraft relies on — it is simply shedding momentum into a wide wake. That omission shows up directly in the pitching moment data below.

Run 1 — AoA = 0° (Baseline)

Figure 4a — Force convergence at α = 0°. Total force z (orange, lift) converges to +714 N; total force x (cyan, drag) converges to −394 N.

Figure 4b — Moment convergence at α = 0°. Total moment y (red, pitching moment about CFD origin) converges to +1,018 N·m after ~400 iterations.

 

 

The solver converges cleanly after roughly 300 iterations. Converged values at iteration 1,000:

Quantity	Value
Fx (drag)	−393.5 N
Fz (lift)	+713.5 N
My about CFD origin	+1,017.9 N·m

Non-dimensionalized, with the pitching moment shifted to the CG (x_CG = 1.788 m from the CFD origin):

• C_D0 = 0.112
• C_L0 = 0.203
• C_m0 = −0.163 (about CG)

Positive lift at zero AoA is expected — the wing has camber and a slight built-in incidence. The slightly negative C_m0 (nose-down tendency at trim) is a healthy sign for a tailless design.

Run 2 — AoA = +5°

Figure 5a — Force convergence at α = 5°. Lift (orange) climbs to +1,291 N — an 81% increase over the baseline.

 

 

Figure 5b — Moment convergence at α = 5°. Total moment y (red) converges to +1,710 N·m — a 68% jump from Run 1, indicating growing pitch instability.

 

 

Quantity	Value
Wind-axis lift L	+1,318.3 N
My about CG	−597.8 N·m
CL5	0.376
Cm5	−0.487

Finite-difference slopes from Run 1 → Run 2:

C_Lα (0°→5°) = 1.98 /rad
C_mα (0°→5°) = +3.72 /rad (ENU/Gazebo sign convention)

A positive C_mα means the aircraft is pitch unstable: a nose-up perturbation generates additional nose-up moment, amplifying rather than damping the disturbance. This is the geometric consequence of the center of pressure sitting ahead of the center of gravity with no tail surface to counteract it.

Run 3 — AoA = +15°

Figure 6a — Force convergence at α = 15°. Lift (orange) reaches +2,584 N. Note the slightly longer settling time vs. Runs 1–2, consistent with a more complex separated-flow transient.

Figure 6b — Moment convergence at α = 15°. Total moment y (red) converges to +3,255 N·m. The slope from Run 2 to Run 3 is nearly identical to Run 1–2 — no nonlinear relief.

The 15° run addresses one specific question: does the instability soften at higher AoA? Some tailless platforms benefit from leading-edge vortex attachment at moderate angles that shifts the center of pressure rearward, reducing (or even reversing) the pitch instability. If MARID has this effect, the SAS requirement relaxes at the extremes of the flight envelope.

Quantity	Value
Wind-axis lift L	+2,577.7 N
My about CG	−1,365.3 N·m
CL15	0.735
Cm15	−1.113

Finite-difference slopes from Run 2 → Run 3:

C_Lα (5°→15°) = 2.06 /rad
C_mα (5°→15°) = +3.58 /rad (ENU/Gazebo sign convention)

The lift curve is nearly linear (2.06 vs. 1.98 /rad) — no stall signature in this AoA range. More critically, the instability barely changed: C_mα dropped only from 3.72 to 3.58 /rad. There is no nonlinear relief at high AoA.

Summary of Stability Derivatives

AoA range	CLα (/rad)	Cmα ENU (/rad)
0° → 5°	1.98	+3.72
5° → 15°	2.06	+3.58

The pitch instability is uniform and unrelieved across the tested AoA range. The lift curve slope is consistent with thin-airfoil theory corrected for finite span (theoretical 2π/rad ≈ 6.28 /rad for infinite span, reduced to ~2 /rad by the aspect ratio AR = 6.5 at span efficiency η = 0.97).

Flight Control Implication

An airframe with C_mα > 0 has an unstable real pitch pole. A rough estimate of the divergence time constant is:

τ ≈ √( I_y / (q·S·c·C_mα) )

With I_y = 132.1 kg·m² (full-assembly CAD value from Onshape), q·S·c = 1,226 N·m at 60 m/s, and C_mα = 3.6 /rad:

τ ≈ 173 ms

A 173 ms pitch divergence time constant is too fast for open-loop human control but sits comfortably within the bandwidth of a 50 Hz SAS loop, which has roughly 8–9 control cycles before the divergence reaches a meaningful amplitude. This is tighter than it sounds: with typical actuator lag and sensor latency, there is little phase margin to spare. A Stability Augmentation System (SAS) is still non-negotiable above approximately 40 m/s — below that speed the lower dynamic pressure slows the divergence enough to give the controller ample headroom.

MARID's pitch controller uses a cascaded attitude loop with feedforward from the ESKF (Error-State Kalman Filter) state estimator to provide the required stabilization. These CFD-derived coefficients are now embedded in the Gazebo aerodynamic model via the AdvancedLiftDrag plugin in marid_alt_gazebo.xacro:

• C_Lα = 2.06 /rad
• C_mα,xacro = 2.97 /rad (the plugin's moment arm already contributes 0.614 /rad, so the direct coefficient is reduced to bring the total to 3.58 /rad)

What's Next

Three runs at a single speed and a single plane are a starting point, not a final answer:

• Lateral runs at non-zero sideslip to extract C_Yβ, C_lβ (dihedral effect), C_nβ (weathercock stability)
• Negative AoA to verify symmetric behavior
• Higher-Re runs at cruise altitude conditions
• Wind tunnel cross-check planned for the full-scale prototype phase