气动面元法求解器中的下标越界错误

编程语言 2026-07-09

我得到一个用Python编写的空气动力学涡面板法求解器，我知道它是正确的。我已经尝试把它转换成C 以提高速度，但我的代码并不总是能正常工作。对于合理形状的翼型（如NACA 2412），C程序给出的结果与Python完全一致。然而，对于形状奇特的翼型，如NACA 9999，C程序给出的升力/阻力比是错误的。我发现问题出在翼型生成阶段，即最终的面板数值略有不正确。这会被带入后续计算，导致最终结果也不正确。代码量相当大（我已经无法再缩减）。正确的Python代码是：

import numpy as np
import math

class Panel:
    """
    A single panel (flat segment) on the airfoil surface.
    """
    def __init__(self, xa, ya, xb, yb):
        """
        Initialize a panel from point A to point B.

        Parameters:
        -----------
        xa, ya : float
            Coordinates of start point
        xb, yb : float
            Coordinates of end point
        """
        # print(f"{xa} {ya} {xb} {yb}")
        self.xa = xa
        self.ya = ya
        self.xb = xb
        self.yb = yb

        # Calculate center point (control point where we enforce boundary condition)
        self.xc = (xa + xb) / 2
        self.yc = (ya + yb) / 2

        # Calculate panel length
        self.length = np.sqrt((xb - xa)**2 + (yb - ya)**2)

        # Calculate tangent vector (along the panel)
        # Unit vector pointing from A to B
        self.tx = (xb - xa) / self.length
        self.ty = (yb - ya) / self.length

        # Calculate normal vector (perpendicular to panel, pointing outward)
        # For counterclockwise numbering, rotate tangent 90° to the left
        self.nx = -self.ty
        self.ny = self.tx

        # Calculate panel angle (for later use)
        self.beta = np.arctan2(yb - ya, xb - xa)


class Airfoil:
    def __init__(self, thickness, max_camber_dev, max_camber_location):
        self.name = "NACA ????"
        self.thickness = thickness
        self.max_camber_dev = max_camber_dev
        self.max_camber_location = max_camber_location
        self.generateAirfoil()
        self.create_airfoil_coordinates()
        self.create_panels()

    def generateAirfoil(self, num_points=100) -> None:
        # Create x coordinates from 0 to 1 (leading edge to trailing edge)
        # Use cosine spacing for better resolution at leading/trailing edges
        beta = np.linspace(0, np.pi, num_points)
        x_points = (1 - np.cos(beta)) / 2  # This gives more points near edges

        # Calculate thickness distribution (symmetric airfoil shape)
        # This formula defines the basic teardrop shape
        yt = 5 * self.thickness * (0.2969 * np.sqrt(x_points)
                                   - 0.1260 * x_points
                                   - 0.3516 * x_points**2
                                   + 0.2843 * x_points**3
                                   - 0.1015 * x_points**4)

        # Calculate camber line (this creates the curve/bend in the airfoil)
        yc = np.zeros_like(x_points)  # Start with straight line
        dyc_dx = np.zeros_like(x_points)  # Slope of camber line

        if self.max_camber_dev > 0:  # Only add camber if m > 0 (otherwise it's symmetric)
            # Front section (0 to p)
            mask_front = x_points < self.max_camber_location
            yc[mask_front] = self.max_camber_dev / self.max_camber_location**2 * (2 * self.max_camber_location * x_points[mask_front] - x_points[mask_front]**2)
            dyc_dx[mask_front] = 2 * self.max_camber_dev / self.max_camber_location**2 * (self.max_camber_location - x_points[mask_front])

            # Rear section (p to 1)
            mask_rear = x_points >= self.max_camber_location
            yc[mask_rear] = self.max_camber_dev / (1 - self.max_camber_location)**2 * ((1 - 2*self.max_camber_location) + 2*self.max_camber_location*x_points[mask_rear] - x_points[mask_rear]**2)
            dyc_dx[mask_rear] = 2 * self.max_camber_dev / (1 - self.max_camber_location)**2 * (self.max_camber_location - x_points[mask_rear])

        # Calculate angle of camber line
        theta = np.arctan(dyc_dx)

        # Upper surface coordinates
        # (camber line + thickness perpendicular to camber)
        self.x_upper = x_points - yt * np.sin(theta)
        self.y_upper = yc + yt * np.cos(theta)

        # Lower surface coordinates
        # (camber line - thickness perpendicular to camber)
        self.x_lower = x_points + yt * np.sin(theta)
        self.y_lower = yc - yt * np.cos(theta)

    def create_airfoil_coordinates(self) -> None:
        """
        Create airfoil as a closed loop, starting from trailing edge,
        going over the top, around the leading edge, and back underneath.
        """
        # Reverse upper surface (to go from trailing edge to leading edge)
        x_u = self.x_upper[::-1]
        y_u = self.y_upper[::-1]

        # Concatenate: upper surface (TE to LE) + lower surface (LE to TE)
        # Remove duplicate point at trailing edge
        x = np.concatenate([x_u, self.x_lower[1:]])
        y = np.concatenate([y_u, self.y_lower[1:]])

        self.x_points = x
        self.y_points = y

    def create_panels(self) -> None:
        """
        Create panels from airfoil coordinates.
        """
        panels = []

        # Create one panel between each pair of consecutive points
        num_panels = len(self.x_points) - 1
        print(f"Num panels: {num_panels}")

        for i in range(num_panels):
            panel = Panel(self.x_points[i], self.y_points[i], self.x_points[i+1], self.y_points[i+1])
            panels.append(panel)

        self.panels = panels

并且C 代码是

#include <math.h>
#include <stdio.h>
#include <string.h>

#define NUM_POINTS 100
#define NUM_PANELS (NUM_POINTS * 2 - 2)
#define HALF_NUM_PANELS (NUM_PANELS/2)

typedef struct {
    double xa, ya;
    double xb, yb;
    double xc, yc;
    double length;
    double tx, ty;
    double nx, ny;
    double beta;
    double sigma;
    double vt, cp;
} Panel;

typedef struct {
    char name[64];

    double thickness;
    double max_camber_dev;
    double max_camber_location;

    double x_points[NUM_POINTS*2-1];
    double y_points[NUM_POINTS*2-1];

    Panel panels[NUM_PANELS];
} Airfoil;

static void init_panel(Panel* p,
                       double xa, double ya,
                       double xb, double yb) {

    p->xa = xa;
    p->ya = ya;
    p->xb = xb;
    p->yb = yb;

    p->xc = (xa + xb) / 2.0;
    p->yc = (ya + yb) / 2.0;

    p->length = sqrt((xb - xa)*(xb - xa) +
                     (yb - ya)*(yb - ya));

    p->tx = (xb - xa) / p->length;
    p->ty = (yb - ya) / p->length;

    p->nx = -p->ty;
    p->ny = p->tx;

    p->beta = atan2(yb - ya, xb - xa);
}

void init_airfoil_empty(Airfoil *a) {
    memset(a, 0, sizeof(Airfoil));
}

void create_panels(Airfoil* a) {
    printf("Num panels: %d\n", NUM_PANELS);
    for (int i = 0; i < NUM_PANELS; i++) {
        init_panel(&a->panels[i],
                   a->x_points[i],
                   a->y_points[i],
                   a->x_points[i + 1],
                   a->y_points[i + 1]);
    }
}

void init_airfoil_values(Airfoil *a,
                         double thickness,
                         double max_camber_dev,
                         double max_camber_location) {

    init_airfoil_empty(a);

    strcpy(a->name, "Generated NACA");

    a->thickness = thickness;
    a->max_camber_dev = max_camber_dev;
    a->max_camber_location = max_camber_location;

    double x_points[NUM_POINTS];
    double yt[NUM_POINTS];
    double yc[NUM_POINTS] = {0};
    double dyc_dx[NUM_POINTS] = {0};
    for (int i = 0; i < NUM_POINTS; i ++) {
        x_points[i] = (1 - cos(M_PI * i / (NUM_POINTS - 1))) / 2;
        yt[i] = 5 * thickness * (0.2969 * sqrt(x_points[i])
                              - 0.1260 * x_points[i]
                              - 0.3516 * x_points[i] * x_points[i]
                              + 0.2843 * x_points[i] * x_points[i] * x_points[i]
                              - 0.1015 * x_points[i] * x_points[i] * x_points[i] * x_points[i]);
    }

    if (max_camber_dev > 0) {
        int i;
        for (i = 0; x_points[i] < max_camber_location; i++) {
            yc[i] = max_camber_dev / (max_camber_location * max_camber_location) * (2 * max_camber_location * x_points[i] - (x_points[i] * x_points[i]));
            dyc_dx[i] = 2 * max_camber_dev / (max_camber_location * max_camber_location) * (max_camber_location - x_points[i]);
        }
        for (; i < NUM_POINTS; i++) {
            yc[i] = max_camber_dev / ((1 - max_camber_location) * (1 - max_camber_location)) * ((1 - 2*max_camber_location) + 2*max_camber_location*x_points[i] - x_points[i] * x_points[i]);
            dyc_dx[i] = 2 * max_camber_dev / ((1 - max_camber_location) * (1 - max_camber_location)) * (max_camber_location - x_points[i]);
        }
    }

    double theta[NUM_POINTS];
    double x_upper[NUM_POINTS];
    double y_upper[NUM_POINTS];
    double x_lower[NUM_POINTS];
    double y_lower[NUM_POINTS];
    for (int i = 0; i < NUM_POINTS; i++) {
        theta[i] = atan(dyc_dx[i]);
        x_upper[i] = x_points[i] - yt[i] * sin(theta[i]);
        y_upper[i] = yc[i] + yt[i] * cos(theta[i]);
        x_lower[i] = x_points[i] + yt[i] * sin(theta[i]);
        y_lower[i] = yc[i] - yt[i] * cos(theta[i]);
    }

    for (int i = 0; i < NUM_POINTS; i++) {
        a->x_points[i] = x_upper[NUM_POINTS-i-1];
        a->y_points[i] = y_upper[NUM_POINTS-i-1];
    }
    for (int i = 1; i < NUM_POINTS; i++) {
        a->x_points[NUM_POINTS+i-1] = x_lower[i];
        a->y_points[NUM_POINTS+i-1] = y_lower[i];
    }

    create_panels(a);
}

对于非寻常翼型，C版本的最终面板略有差异。准确地说，在C 中，最终面板的不同参数是 0.990662711025 -0.004897731951 0.990412286936 -0.004747205633 0.8570821543420353 -0.5151797557245585，它们的顺序是 panel.xc, panel.yc, panel.xa, panel.ya, panel.tx, panel.ty。在Python中，它们是 0.990662711025 -0.004897731951 0.990412286936 -0.004747205633 0.8570821543418085 -0.5151797557249356。确实存在一些差异，但非常小。这些差异只有在导出到大约18位小数时才会显现。当我把C 的最终面板导入到许多双精度浮点数的Python环境中时，Python也会给出错误的答案，因此这一定与翼型构造有关。在此先行致谢，感谢任何我能得到的帮助。