滑翔（Gliding）是人类最古老的飞行方式之一。从纸飞机到翼装飞行，其核心都依赖于空气动力学中的升力与阻力平衡。本文将用 Java 构建一个滑翔模拟器，重点讲解升力系数计算、空气阻力建模以及基于梯度下降的最优滑翔路径规划算法。
一、空气动力学基础模型
滑翔体在空气中运动时，主要受三种力作用：
- 重力 $G = mg$，方向竖直向下
- 升力 $L = \frac{1}{2} \rho v^2 S C_L$，方向垂直于来流速度向上
- 阻力 $D = \frac{1}{2} \rho v^2 S C_D$，方向与速度相反
其中 $\rho$ 为空气密度，$v$ 为空速，$S$ 为机翼面积，$C_L$ 和 $C_D$ 分别为升力系数与阻力系数。二者通过极曲线（polar curve）关联：
$$C_D = C_{D0} + k \cdot C_L^2$$
$C_{D0}$ 为零升阻力系数，$k$ 为诱导阻力因子。这个二次关系是路径优化问题的关键约束。
二、核心数据结构
/**
* 三维向量类，描述滑翔体的位置、速度、受力
*/
public class Vector3D {
public double x, y, z;
public Vector3D(double x, double y, double z) {
this.x = x; this.y = y; this.z = z;
}
public Vector3D add(Vector3D v) {
return new Vector3D(x + v.x, y + v.y, z + v.z);
}
public Vector3D scale(double s) {
return new Vector3D(x * s, y * s, z * s);
}
public double length() {
return Math.sqrt(x*x + y*y + z*z);
}
public Vector3D normalize() {
double len = length();
if (len < 1e-9) return new Vector3D(0, 0, 0);
return new Vector3D(x/len, y/len, z/len);
}
@Override
public String toString() {
return String.format("(%.2f, %.2f, %.2f)", x, y, z);
}
}
/**
* 滑翔机物理参数
*/
public class Glider {
// 位置与速度（米）
public Vector3D position;
public Vector3D velocity;
// 物理参数
public static final double MASS = 5.0; // 质量 kg
public static final double WING_AREA = 0.8; // 机翼面积 m²
public static final double RHO = 1.225; // 海平面空气密度 kg/m³
public static final double CD0 = 0.015; // 零升阻力系数
public static final double K = 0.045; // 诱导阻力因子
public static final double G = 9.81; // 重力加速度
// 控制量：当前攻角（Angle of Attack，弧度）
public double angleOfAttack;
public Glider(double x, double y, double z) {
this.position = new Vector3D(x, y, z);
this.velocity = new Vector3D(15, 0, -2); // 初始向前下方滑翔
this.angleOfAttack = Math.toRadians(5);
}
}
三、升力与阻力计算
根据攻角 $\alpha$ 计算升力系数 $C_L$，再通过极曲线求 $C_D$，最终得出气动力矢量。
public class Aerodynamics {
/**
* 根据攻角计算升力系数（采用小攻角线性近似）
* @param alpha 攻角（弧度）
* @return 升力系数
*/
public static double liftCoefficient(double alpha) {
// CL = 2π * alpha（薄翼理论，小攻角范围有效）
return 2 * Math.PI * alpha;
}
/**
* 根据升力系数计算阻力系数（极曲线公式）
* @param cl 升力系数
* @return 阻力系数
*/
public static double dragCoefficient(double cl) {
return Glider.CD0 + Glider.K * cl * cl;
}
/**
* 计算当前状态下滑翔机受到的气动力
* @param glider 滑翔机状态
* @return [升力矢量, 阻力矢量]
*/
public static Vector3D[] computeForces(Glider glider) {
double v = glider.velocity.length();
if (v < 0.1) v = 0.1; // 防止除零
// 速度方向单位向量
Vector3D vDir = glider.velocity.normalize();
// 升力系数与阻力系数
double cl = liftCoefficient(glider.angleOfAttack);
double cd = dragCoefficient(cl);
// 动压 q = 0.5 * rho * v^2
double q = 0.5 * Glider.RHO * v * v;
// 升力大小
double liftMag = q * Glider.WING_AREA * cl;
// 阻力大小
double dragMag = q * Glider.WING_AREA * cd;
// 升力方向：垂直于速度方向且在竖直平面内
// 简化处理：假设滑翔机对称飞行，升力方向 = 归一化(速度水平分量反向, 0, 速度水平大小)
double vh = Math.sqrt(vDir.x * vDir.x + vDir.y * vDir.y);
Vector3D liftDir;
if (vh < 1e-6) {
liftDir = new Vector3D(0, 0, 1);
} else {
liftDir = new Vector3D(-vDir.x * vDir.z / vh, -vDir.y * vDir.z / vh, vh).normalize();
}
Vector3D lift = liftDir.scale(liftMag);
Vector3D drag = vDir.scale(-dragMag); // 阻力与速度反向
return new Vector3D[]{lift, drag};
}
}
四、物理积分与状态更新
采用欧拉法进行每帧数值积分，将加速度积分得到速度，速度积分得到位移。
public class GliderSimulator {
private Glider glider;
public static final double DT = 0.05; // 时间步长 50ms
public GliderSimulator(Glider g) {
this.glider = g;
}
/**
* 单步物理积分（欧拉法）
*/
public void step() {
Vector3D[] forces = Aerodynamics.computeForces(glider);
Vector3D lift = forces[0];
Vector3D drag = forces[1];
Vector3D gravity = new Vector3D(0, 0, -Glider.MASS * Glider.G);
// 合力
Vector3D netForce = lift.add(drag).add(gravity);
// 加速度 a = F/m
Vector3D acceleration = netForce.scale(1.0 / Glider.MASS);
// 欧拉积分：v = v + a*dt, p = p + v*dt
glider.velocity = glider.velocity.add(acceleration.scale(DT));
glider.position = glider.position.add(glider.velocity.scale(DT));
}
/**
* 模拟滑翔轨迹（固定攻角）
* @param steps 模拟步数
* @return 轨迹点列表
*/
public java.util.List<Vector3D> simulateTrajectory(int steps) {
java.util.List<Vector3D> traj = new java.util.ArrayList<>();
for (int i = 0; i < steps; i++) {
step();
traj.add(new Vector3D(glider.position.x, glider.position.y, glider.position.z));
if (glider.position.z <= 0) break; // 触地
}
return traj;
}
}
五、最优滑翔路径：梯度下降法
滑翔的核心优化目标是：给定起始高度，如何选择攻角序列使水平飞行距离最远？ 这是一个连续决策问题，我们将其离散化为 $N$ 个阶段，每阶段选择一个攻角 $\alpha_i$，然后用梯度下降迭代优化。
public class OptimalGlidePath {
/**
* 计算一条攻角序列对应的滑翔水平距离
* @param alphas 攻角序列（弧度）
* @param startPos 起始位置
* @param startVel 起始速度
* @return 最终水平距离（x方向）
*/
public static double evaluateDistance(double[] alphas, Vector3D startPos, Vector3D startVel) {
Glider g = new Glider(startPos.x, startPos.y, startPos.z);
g.velocity = new Vector3D(startVel.x, startVel.y, startVel.z);
GliderSimulator sim = new GliderSimulator(g);
int stepsPerSegment = 20; // 每段攻角持续20步
for (int seg = 0; seg < alphas.length; seg++) {
g.angleOfAttack = alphas[seg];
for (int s = 0; s < stepsPerSegment; s++) {
sim.step();
if (g.position.z <= 0) return g.position.x;
}
}
return g.position.x;
}
/**
* 使用梯度下降优化攻角序列
* @param segments 离散段数
* @param iterations 迭代次数
* @param lr 学习率
* @return 最优攻角序列
*/
public static double[] optimize(int segments, int iterations, double lr,
Vector3D startPos, Vector3D startVel) {
double[] alphas = new double[segments];
// 初始化为小攻角
for (int i = 0; i < segments; i++) alphas[i] = Math.toRadians(4);
double epsilon = Math.toRadians(0.5); // 数值微分步长
for (int iter = 0; iter < iterations; iter++) {
double baseScore = evaluateDistance(alphas, startPos, startVel);
double[] gradients = new double[segments];
// 对每个段进行数值微分求梯度
for (int i = 0; i < segments; i++) {
alphas[i] += epsilon;
double plusScore = evaluateDistance(alphas, startPos, startVel);
alphas[i] -= epsilon;
gradients[i] = (plusScore - baseScore) / epsilon;
}
// 梯度上升（距离越大越好）
for (int i = 0; i < segments; i++) {
alphas[i] += lr * gradients[i];
// 限制攻角在合理范围 [0°, 15°]
alphas[i] = Math.max(0, Math.min(Math.toRadians(15), alphas[i]));
}
if (iter % 50 == 0) {
System.out.printf("Iter %d: distance=%.2f m%n", iter, baseScore);
}
}
return alphas;
}
}
六、完整演示程序
public class GlideDemo {
public static void main(String[] args) {
// 起始条件：高度100米，水平初速15m/s
Vector3D startPos = new Vector3D(0, 0, 100);
Vector3D startVel = new Vector3D(15, 0, -1);
System.out.println("=== 固定攻角滑翔模拟 ===");
Glider g1 = new Glider(startPos.x, startPos.y, startPos.z);
g1.velocity = new Vector3D(startVel.x, startVel.y, startVel.z);
g1.angleOfAttack = Math.toRadians(5);
GliderSimulator sim1 = new GliderSimulator(g1);
java.util.List<Vector3D> traj1 = sim1.simulateTrajectory(2000);
System.out.println("固定5°攻角轨迹点数: " + traj1.size());
System.out.println("最终位置: " + g1.position);
System.out.println("水平飞行距离: " + String.format("%.2f", g1.position.x) + " m");
System.out.println("\n=== 梯度下降优化攻角序列 ===");
double[] optimalAlphas = OptimalGlidePath.optimize(
10, // 10段攻角
200, // 200次迭代
0.001, // 学习率
startPos, startVel
);
System.out.println("\n最优攻角序列（度）:");
for (int i = 0; i < optimalAlphas.length; i++) {
System.out.printf(" 段%d: %.2f°%n", i, Math.toDegrees(optimalAlphas[i]));
}
// 用最优攻角重新模拟
Glider g2 = new Glider(startPos.x, startPos.y, startPos.z);
g2.velocity = new Vector3D(startVel.x, startVel.y, startVel.z);
GliderSimulator sim2 = new GliderSimulator(g2);
int stepsPerSegment = 20;
for (int seg = 0; seg < optimalAlphas.length; seg++) {
g2.angleOfAttack = optimalAlphas[seg];
for (int s = 0; s < stepsPerSegment; s++) {
sim2.step();
if (g2.position.z <= 0) break;
}
}
System.out.println("\n优化后水平飞行距离: " + String.format("%.2f", g2.position.x) + " m");
}
}
七、算法复杂度分析
| 模块 | 复杂度 | 说明 |
|---|---|---|
| 气动力计算 | O(1) | 固定公式运算 |
| 单步物理积分 | O(1) | 向量加减乘除 |
| 轨迹模拟 | O(T) | T 为模拟步数 |
| 单轮距离评估 | O(T) | 完整模拟一次 |
| 梯度下降优化 | O(N × T × I) | N段攻角，I次迭代，每轮需N次数值微分 |
以 10 段攻角、200 次迭代为例，总共约 10 × 200 × 200 = 40 万步物理积分。Java 单线程可在 1 秒内完成，满足离线规划需求。若需实时规划，可改用遗传算法或模型预测控制（MPC）降低计算量。
八、扩展方向
- 热气流利用：引入上升气流场，目标变为最大化滞空时间或飞行距离，需使用动态规划或A*搜索在气流网格中寻找最优路径。
- 三维航迹规划：加入转向控制，利用Dubins曲线或B样条生成平滑轨迹。
- 风切变建模：在不同高度层设置不同风速，模拟真实大气环境。
滑翔模拟的精髓在于将连续物理过程离散化，并用梯度下降寻找最优控制序列。这套方法论不仅适用于飞行模拟，也可迁移到机器人运动规划与自动驾驶轨迹优化领域。
