Моделирование обтекания тела Ахмеда в OpenFOAM

1. Математическая постановка задачи

Рассмотрим стационарное обтекание тела Ахмеда потоком вязкой несжимаемой жидкости. Набегающий поток со скоростью U∞ направлен вдоль оси Ox. Тело Ахмеда представляет собой упрощенную геометрическую модель легкового автомобиля.

Параметр	Обозначение	Значение
Скорость набегающего потока	U∞	40 м/с
Характерная длина тела	L	1.044 м
Характерная площадь	Aref	0.112032 м²
Кинематическая вязкость воздуха	ν	1.5·10⁻⁵ м²/с
Плотность воздуха	ρ	1.225 кг/м³

Уравнения Навье–Стокса (RANS):

        (U · ∇)U = −(1/ρ)∇P + ∇·(ν∇U) − ∇·(u'u')

        ∇·U = 0

Модель турбулентности k-ω SST:

        D(ρk)/Dt = P̃ + ∇·(ρDk∇k) − ρβ*ωk

        D(ρω)/Dt = (γ/νt)P̃ + ∇·(ρDω∇ω) − ρβω² − ρ(1−F1)CDkω

2. Создание проекта OpenFOAM

Структура проекта ahmedBody:

./ahmedBody/
    |-- 0/
    |   |-- k
    |   |-- nut
    |   |-- omega
    |   |-- p
    |   `-- U
    |-- constant/
    |   |-- triSurface/
    |   |   `-- ahmedBody.stl
    |   |-- transportProperties
    |   `-- turbulenceProperties
    `-- system/
        |-- blockMeshDict
        |-- controlDict
        |-- decomposeParDict
        |-- fvSchemes
        |-- fvSolution
        |-- snappyHexMeshDict
        `-- surfaceFeatureExtractDict

Команды для запуска расчета:

openfoam2412
surfaceFeatureExtract
blockMesh
snappyHexMesh -overwrite
simpleFoam
paraFoam

Параллельный запуск:

decomposePar
mpirun -np 6 simpleFoam -parallel
reconstructPar

3. Настройка расчетной сетки

surfaceFeatureExtractDict:

ahmedBody.stl
{
    extractionMethod    extractFromSurface;
    extractFromSurfaceCoeffs
    {
        includedAngle   150;
    }
    writeObj            yes;
}

blockMeshDict (фрагмент):

scale 1;
vertices
(
    (-4 0 -1)   // 0
    (6 0 -1)    // 1
    (6 4 -1)    // 2
    (-4 4 -1)   // 3
    (-4 0 1)    // 4
    (6 0 1)     // 5
    (6 4 1)     // 6
    (-4 4 1)    // 7
);
blocks
(
    hex (0 1 2 3 4 5 6 7) (50 20 10) simpleGrading (1 1 1)
);

snappyHexMeshDict (фрагмент):

castellatedMesh true;
snap            true;
addLayers       true;

geometry
{
    ahmedBody.stl { type triSurfaceMesh; name ahmedBody; }
    refinementBoxBig { type searchableBox; min (-2.5 0.0 -0.7); max (3.5 2.0 0.7); }
    refinementBox    { type searchableBox; min (-1.5 0.0 -0.5); max (0.5 0.7 0.5); }
}

4. Настройка граничных условий

Поле давления (0/p):

dimensions      [0 2 -2 0 0 0 0];
internalField   uniform 0;
boundaryField
{
    inlet       { type zeroGradient; }
    outlet      { type fixedValue; value uniform 0; }
    ground      { type zeroGradient; }
    sides       { type slip; }
    top         { type slip; }
    ahmedBody   { type zeroGradient; }
}

Поле скорости (0/U):

dimensions      [0 1 -1 0 0 0 0];
internalField   uniform (40 0 0);
boundaryField
{
    inlet       { type fixedValue; value $internalField; }
    outlet      { type inletOutlet; inletValue uniform (0 0 0); value $internalField; }
    ground      { type fixedValue; value uniform (40 0 0); }
    sides       { type slip; }
    top         { type slip; }
    ahmedBody   { type fixedValue; value uniform (0 0 0); }
}

Поля турбулентности (k, omega, nut):

// k
dimensions [0 2 -2 0 0 0 0]; internalField uniform 0.24;
boundaryField { inlet fixedValue; outlet inletOutlet; ground kqRWallFunction; ahmedBody kqRWallFunction; sides slip; top slip; }

// omega
dimensions [0 0 -1 0 0 0 0]; internalField uniform 12.78;
boundaryField { inlet fixedValue; outlet inletOutlet; ground omegaWallFunction; ahmedBody omegaWallFunction; sides slip; top slip; }

// nut
dimensions [0 2 -1 0 0 0 0]; internalField uniform 0;
boundaryField { inlet calculated; outlet calculated; ground nutkWallFunction; ahmedBody nutkWallFunction; sides calculated; top calculated; }

5. Настройка методов дискретизации

controlDict:

application     simpleFoam;
startFrom       latestTime;
startTime       0;
stopAt          endTime;
endTime         2000;
deltaT          1;
writeControl    timeStep;
writeInterval   100;

fvSchemes:

ddtSchemes { default steadyState; }
gradSchemes { default Gauss linear; }
divSchemes
{
    div(phi,U)      bounded Gauss GammaV 0.3;
    div(phi,k)      bounded Gauss upwind;
    div(phi,omega)  bounded Gauss upwind;
}
laplacianSchemes { default Gauss linear corrected; }

fvSolution:

p { solver PBiCGStab; preconditioner GAMG; tolerance 1e-05; relTol 0.001; }
U { solver PBiCG; preconditioner DILU; tolerance 1e-06; relTol 0.001; }
k { solver PBiCGStab; preconditioner DILU; tolerance 1e-05; relTol 0.1; }
omega { solver PBiCGStab; preconditioner DILU; tolerance 1e-05; relTol 0.1; }
SIMPLE { nNonOrthogonalCorrectors 0; consistent yes; }
relaxationFactors { fields { p 0.5; } equations { U 0.7; k 0.7; omega 0.7; } }

6. Коэффициенты аэродинамических сил

Формулы:

        CD = Fx / (0.5 · ρ · U∞² · Aref)

        CL = Fy / (0.5 · ρ · U∞² · Aref)

Python-скрипт для построения графиков:

import pandas as pd
import matplotlib.pyplot as plt

df = pd.read_csv('postProcessing/forces/0/coefficient.dat', comment='#', delim_whitespace=True,
                 names=['Time','Cd','Cd(f)','Cd(r)','Cl','Cl(f)','Cl(r)',
                        'CmPitch','CmRoll','CmYaw','Cs','Cs(f)','Cs(r)'])

df_cut = df[df['Time'] >= 50]
cd_mean = df_cut['Cd'].iloc[-100:].mean()
cl_mean = df_cut['Cl'].iloc[-100:].mean()

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8))
ax1.plot(df_cut['Time'], df_cut['Cd'], 'b-', label='CD')
ax1.axhline(y=cd_mean, color='r', linestyle='--', label=f'CD mean: {cd_mean:.4f}')
ax1.legend()
ax2.plot(df_cut['Time'], df_cut['Cl'], 'b-', label='CL')
ax2.axhline(y=cl_mean, color='r', linestyle='--', label=f'CL mean: {cl_mean:.4f}')
ax2.legend()
plt.savefig('coefficients.png', dpi=300)
plt.show()

print(f"CD: {cd_mean:.6f}, CL: {cl_mean:.6f}")

7. Результаты

Источник	CD	CL
Эксперимент (Ahmed et al., 1984)	0.299	0.345
Настоящий расчет (OpenFOAM)	0.313	0.334

✓ Результаты показывают хорошее согласие с экспериментом. Отличия составляют менее 5%.

8. Переход к реальному автомобилю

Отработанная методология переносится на реалистичную геометрию автомобиля. Замена STL-геометрии, все настройки решателя и схемы остаются неизменными. Расчет на детальной сетке (>10 млн ячеек) требует значительных ресурсов.

💡 Примечание: Пользователь, освоивший методологию на теле Ахмеда, способен перенести полученные навыки на расчеты любых автомобильных форм.

startup