Цифровой двойник на примере гидроамортизатора стойки шасси самолета: от физико-математической модели к численному ядру

Это четвертая и замыкающая статья первого цикла о цифровых двойниках, CFD / вычислительной гидрогазодинамике и инженерном моделировании. В первых трех публикациях были последовательно рассмотрены вычислительная база CFD и HPC / высокопроизводительных вычислений, место цифровых двойников в современной инженерии и структурная декомпозиция цифрового двойника изделия на примере системы шасси коммерческого самолета. В этой статье цикл замыкается переходом к предметной физике: от архитектуры цифрового двойника — к физико-математической модели и численному ядру цифрового двойника гидроамортизатора. Одновременно эта статья открывает следующий контур публикаций — о развитии цифровых двойников различных изделий и технологий в эпоху супервычислений и AI / искусственного интеллекта.

Предыдущие статьи цикла:

• Статья 1. Развитие математического моделирования и цифровых моделей. CFD-2030
• Статья 2. Цифровые двойники в современной инженерии
• Статья 3. Структурное построение цифрового двойника изделия

Аннотация

Статья является четвёртой публикацией цикла и непосредственно продолжает статью 3, в которой гидроамортизатор был выделен как самостоятельный двойник одной стойки шасси. Настоящая работа переводит этот вывод в предметную плоскость: от архитектурного уровня — к физико-математической постановке, численной схеме и тестовому сценарию верификации и валидации.

В статье предложена переработанная физико-математическая основа цифрового двойника гидроамортизатора для высокоскоростного демпфирования в частично заполненном цилиндре со свободным газовым объёмом. Исходная трёхмерная постановка задачи высокоскоростного демпфирования была основана на уравнениях Навье–Стокса, VOF методе объёмных долей жидкости, RANS осреднённых по Рейнольдсу уравнениях Навье–Стокса и деформируемой подвижной сетке ^[1]. В настоящей работе физическая модель расширена до неизотермической системы Навье–Стокса–Фурье, учитывающей баланс энергии, тепломассоперенос, слабую сжимаемость газа и жидкости, температурно-зависимую динамическую вязкость, теплопроводность и обратную тепловую связь через гидродинамическую силу.

Математическое моделирование в CFD зависит от многих первоначально заданных параметров, используемых вычислительных мощностей и экономической целесообразности. В данной статье рассматривается подход к моделированию, основанный на деформации расчётной сетки без изменения её топологии. Движение поршня учитывается через перемещение узлов сетки, скорости граней и относительные массовые потоки через деформируемые контрольные объёмы. Для перемещения внутренних узлов сетки применяется метод обратного взвешивания по расстоянию IDW ^[8]. Для предотвращения чрезмерной деформации ячеек узлы на неподвижных поверхностях перемещаются тангенциально по поверхности, сохраняя исходную геометрию кожуха и направляющих элементов.

Течение рассматривается как слабосжимаемое; диапазон чисел Рейнольдса в характерных дросселирующих зазорах составляет $10^5$–$10^6$. Межфазная граница описывается односкоростной VOF-моделью. Испарение, кипение и конденсация исключаются; межфазное взаимодействие в базовом расчёте задаётся тепловым обменом, а растворение воздуха в жидкости может включаться отдельным скаляром. Для замыкания осреднённой системы уравнений используется SST $k$–$\omega$ модель переноса касательных напряжений как базовая RANS-модель ^{[6, 7]}.

Для численного интегрирования предложен полностью неявный конечно-объёмный алгоритм на деформируемой сетке с сохранением топологии. По времени применяется трёхслойная схема второго порядка BDF2. Конвективные члены в уравнениях движения, энергии и переноса турбулентных параметров аппроксимируются линейной противоточной схемой LUD. Уравнение переноса объёмной доли фаз дискретизируется высокоточной схемой захвата границы фаз HRIC, уменьшающей численную диффузию границы фаз. Учёт силы тяжести в давлении и потоках выполняется через поправку объёмных сил ^{[13, 14]}. В результате получена согласованная система уравнений, алгоритм её решения, процедура проверки сходимости и программа верификации и валидации цифрового двойника.

Ключевые слова: гидроамортизатор, высокоскоростное демпфирование, цифровой двойник, Навье–Стокса–Фурье, VOF / Volume of Fluid / метод объёмных долей жидкости, SST $k$–$\omega$ / Shear Stress Transport / перенос касательных напряжений, деформируемая сетка без изменения топологии, IDW / Inverse Distance Weighting / метод обратного взвешивания по расстоянию, fully implicit / полностью неявная схема, BDF2, LUD, HRIC, weakly compressible / слабосжимаемая среда, тепломассоперенос, температурная вязкость.

1. Введение и цель работы

В статье 3 было показано, что гидроамортизатор следует рассматривать не как изолированную деталь, а как самостоятельный двойник составной части одной стойки шасси. Именно поэтому следующей логической ступенью становится отдельная статья, посвящённая уже самому гидроамортизатору: его физике, математической модели, численной реализации и роли в цифровом двойнике стойки и всей системы шасси. Если читать цикл последовательно, то эта публикация опирается на три предыдущих материала: о развитии CFD / вычислительной гидрогазодинамики и HPC / высокопроизводительных вычислений, о цифровых двойниках в современной инженерии и о структурной декомпозиции цифрового двойника изделия на примере системы шасси.

Гидроамортизатор преобразует кинетическую энергию подвижных частей в работу течения жидкости через зазоры, канавки и отверстия, а затем — в тепло. Для высокоскоростного режима характерны большие градиенты давления, турбулентные течения, перемещение свободной поверхности, вовлечение газа, локальное распыление и выраженная чувствительность к начальному уровню заполнения. В исходной трёхмерной работе задача была сформулирована как численное моделирование частично заполненного гидравлического цилиндра с дросселирующим профилем; свободная поверхность отслеживалась VOF-методом, а подвижные части учитывались на деформируемой сетке ^[1].

Изотермическое приближение удобно для первого инженерного расчёта закона движения поршня, но оно исключает важный механизм обратной связи: скорость и сдвиг создают вязкую и турбулентную диссипацию; диссипация повышает температуру; температура меняет динамическую вязкость; вязкость меняет распределение давления, расход через дросселирующие зазоры и гидродинамическую силу на поршне. Для гидравлических масел эта цепочка критична, поскольку вязкость обычно быстро уменьшается с ростом температуры. Для воздуха, напротив, динамическая вязкость по закону Сазерленда растёт с температурой ^{[17, 18]}.

Цель статьи 4 — дать инженерно применимую постановку цифрового двойника гидроамортизатора гражданского самолёта на основе уравнений Навье–Стокса–Фурье, VOF и SST $k$–$\omega$, и разработать численную часть с учётом следующих принципиальных требований:

1. движение подвижных частей учитывать через деформацию конечно-объёмной сетки без изменения её топологии;
2. внутренние узлы сетки перемещать методом IDW;
3. узлы, лежащие на неподвижных поверхностях, смещать тангенциально вдоль поверхности так, чтобы геометрия кожуха и направляющих не нарушалась;
4. для подсистемы Навье–Стокса использовать полностью неявный метод решения;
5. по времени использовать трёхслойную схему второго порядка;
6. для конвекции в уравнении импульса и уравнениях SST применять LUD;
7. для переноса объёмной доли фаз применять HRIC;
8. силу тяжести учитывать через согласованную поправку объёмных сил в интерполяции давления и массовых потоков.

2. Постановка задачи и корректировка базовой модели

На уровне вывода физических уравнений конкретная геометрия внутренней профилированной поверхности действительно не определяет вид законов сохранения. Она задаёт область $\Omega(t)$, границы $\partial \Omega(t)$, закон движения поршня, распределение зазоров и сеточную реализацию, но не меняет структуру уравнений массы, импульса и энергии. Поэтому в статье геометрия остаётся параметрической, а конкретный профиль переносится в сеточную модель и расчётный эксперимент.

1. VOF при слабой сжимаемости. Уравнение

   $$\frac{\partial \alpha_q}{\partial t}+\nabla\cdot(\alpha_q\mathbf u)=0$$

   корректно только при постоянной плотности фазы и нулевой дивергенции скорости. Для газа с $\rho_g=p_{abs}/(R_gT)$ и для слабосжимаемой жидкости физически корректной является запись

   $$\frac{\partial(\alpha_q\rho_q)}{\partial t}+\nabla\cdot(\alpha_q\rho_q\mathbf u)=\dot m_q.$$

2. Уравнение состояния газа для абсолютного давления. В формуле $\rho_g=p_{abs}/(R_gT)$ нельзя подставлять манометрическое давление.

3. Уравнение энергии. Если $\mathbf q=-\lambda_{eff}\nabla T$, то приток тепла в объём записывается как $\nabla\cdot(\lambda_{eff}\nabla T)$.

4. Тепловые источники согласованы с RANS-замыканием. Молекулярная диссипация $\tau_\mu:\nabla \mathbf u$ и турбулентный переход энергии $\beta^*\rho k\omega$ должны входить в энергобаланс без двойного счёта ^{[6, 7]}.

5. Изотропная часть турбулентного напряжения согласована с давлением. Член $-\tfrac{2}{3}\rho k I$ может быть либо оставлен в напряжениях, либо поглощён в модифицированное давление; смешивать оба подхода нельзя.

6. Деформируемая конечно-объёмная сетка удовлетворяет геометрическому закону сохранения. На движущейся сетке потоки через грани должны вычисляться через относительную скорость $\mathbf u-\mathbf u_g$, где $\mathbf u_g$ — скорость грани, восстановленная по движению узлов ^{[4, 9, 10]}.

7. Угол положения цилиндра от $0^\circ$ до $90^\circ$ входит через локальный вектор тяжести. Меняется не вид уравнений, а вектор $\mathbf g_{loc}$ и начальная гидростатика.

2.1. Сокращения и их смысл

• NSF / Navier–Stokes–Fourier / уравнения Навье–Стокса–Фурье: связанная система законов сохранения массы, импульса и энергии с теплопроводностью по Фурье.
• VOF / Volume of Fluid / метод объёмных долей жидкости: метод фиксации границы фаз, в котором положение границы фаз описывается объёмной долей одной из фаз ^[2].
• RANS / Reynolds-Averaged Navier–Stokes / осреднённые по Рейнольдсу уравнения Навье–Стокса: инженерный подход к моделированию турбулентности через осреднение и дополнительные замыкающие уравнения.
• SST / Shear Stress Transport / перенос касательных напряжений: семейство моделей $k$–$\omega$, устойчивое в пристеночной области и при неблагоприятном градиенте давления ^{[6, 7]}.
• DMS / deforming mesh strategy / стратегия деформируемой сетки: способ учёта движения границ через перемещение узлов сетки и скоростей граней без изменения топологии сеточной структуры.
• IDW / Inverse Distance Weighting / метод обратного взвешивания по расстоянию: алгебраический метод распределения граничных смещений на внутренние узлы сетки ^[8].
• GCL / Geometric Conservation Law / геометрический закон сохранения: дискретное условие согласованности изменения объёма ячейки со скоростями её граней на движущейся сетке ^{[9, 10]}.
• EOS / Equation of State / уравнение состояния: связь между давлением, температурой, плотностью и составом фазы.
• BFW / body-force-weighted pressure interpolation / интерполяция давления с поправкой объёмных сил: способ вычисления граневого давления и массового потока, сохраняющий гидростатический баланс при наличии тяжести ^{[13, 14]}.
• CFD / Computational Fluid Dynamics / вычислительная гидрогазодинамика: численное решение уравнений течения и теплопереноса.
• ROM / Reduced-Order Model / редуцированная модель: быстродействующая упрощённая модель, аппроксимирующая результаты высокоточных расчётов.
• CSF / Continuum Surface Force / континуальное поверхностное усилие: объёмная запись сил поверхностного натяжения в многофазных моделях типа VOF ^[3].
• LTE / Local Thermal Equilibrium / локальное тепловое равновесие: допущение о равенстве температур фаз в пределах расчётной ячейки.
• ODE / Ordinary Differential Equation / обыкновенное дифференциальное уравнение: уравнение движения поршня как функции времени.
• PDE / Partial Differential Equation / уравнение в частных производных: распределённое уравнение для полей скорости, давления, температуры и фазовой доли.
• GCI / Grid Convergence Index / индекс сеточной сходимости: оценка сеточной неопределённости по результатам расчётов на нескольких сетках ^{[25, 26]}.
• HRIC / High Resolution Interface Capturing / высокоточная схема фиксации границы фаз: компрессивная ограниченная схема дискретизации переноса объёмной доли фаз ^[5].
• LUD / Linear Upwind Differencing / линейная противоточная схема: линейная реконструкция значения на грани по направлению потока, обеспечивающая второй порядок точности на гладких решениях ^{[15, 16]}.
• BDF2 / Backward Differentiation Formula 2 / назад-разностная формула второго порядка: трёхслойная полностью неявная схема по времени второго порядка точности.
• MMS / Method of Manufactured Solutions / метод аналитически изготовленных решений: способ верификации кода путём сравнения с искусственно построенным точным решением.

2.2. Основные обозначения и их физико-математический смысл

2.2.1. Геометрия, сетка и кинематика

$\Omega(t)$ — расчётная область, занятая жидкостью и газом в момент времени $t$.

$\partial\Omega(t)$ — граница расчётной области.

$\Gamma_p$ — часть границы, соответствующая поверхности поршня.

$\Gamma_s$ — неподвижные поверхности: кожух, направляющие стержни и другие неподвижные элементы.

$\mathcal T_h^n$ — топология конечно-объёмной сетки на слое $n$; в работе принимается $\mathcal T_h^{n+1}\equiv \mathcal T_h^n$.

$V_P^n$, $V_P^{n+1}$ — объём контрольного объёма на текущем и следующем временных слоях.

$\mathbf x_i^n$, $\mathbf x_i^{n+1}$ — координаты узла сетки на соседних временных слоях.

$\Delta \mathbf x_i$ — смещение узла сетки за шаг времени.

$\mathbf u_g$ — скорость движения сетки; на грани обозначается $\mathbf u_{g,f}$.

$\phi_{g,f}=\mathbf u_{g,f}\cdot \mathbf S_f$ — объёмный поток, создаваемый движением грани.

$\mathbf S_f$ — вектор площади грани конечного объёма; направлен по внешней нормали.

$\mathbf n$ — единичная внешняя нормаль к поверхности или грани.

$\mathbf e_z$ — единичный вектор вдоль оси движения поршня.

$\Pi_{\Gamma_s}$ — оператор проекции смещённого узла на точную неподвижную поверхность.

$d_{ib}$ — расстояние от внутреннего узла $i$ до граничного узла $b$.

$w_{ib}$ — вес IDW-интерполяции.

$p_{IDW}$ — показатель степени в формуле IDW.

$\chi_b$ — коэффициент тангенциального перераспределения узла на неподвижной стенке.

$V_0$ — полный внутренний объём гидроамортизатора в начальном состоянии.

$S_p$ — поверхность поршня, по которой интегрируется гидродинамическая сила.

$\beta$ — угол положения гидроамортизатора относительно вертикали; рассматривается диапазон $\beta\in[0^\circ,90^\circ]$.

$\mathbf g_0$, $\mathbf g_{loc}$ — глобальный и локальный векторы ускорения силы тяжести.

$R(\beta)$ — матрица поворота, переводящая глобальную систему координат в локальную, связанную с гидроамортизатором.

$\mathbf u$ — массово-осреднённая скорость смеси.

$\mathbf u_w$ — скорость стенки; на неподвижной стенке равна нулю, на поршне равна $V_p\mathbf e_z$.

$\nabla$, $\nabla\cdot$ — оператор градиента и дивергенции.

$D/Dt$ — материальная производная.

$\Delta t$, $\Delta x$ — шаг по времени и характерный шаг сетки.

2.2.2. Фазы, смесь и термодинамика

$\alpha_\ell$, $\alpha_g$ — объёмные доли жидкой и газовой фаз.

$q\in\{\ell,g\}$ — индекс фазы.

$\dot m_q$ — объёмный источник массы фазы; в базовом расчёте $\dot m_q=0$.

$\rho_q$, $\rho$ — плотность фазы и плотность смеси.

$p$, $p_{abs}$, $p_{ref}$ — манометрическое давление, абсолютное давление и опорное давление.

$T$, $T_q$, $T_0$, $T_w$, $T_{ext}$ — температура смеси, фазы, начальная температура, температура стенки и внешней среды.

$h$, $h_q$ — удельная энтальпия смеси и фаз.

$\mathbf q_{eff}$ — эффективный тепловой поток.

$c_p$, $c_{p,q}$ — удельная теплоёмкость при постоянном давлении для смеси и фазы.

$\lambda$, $\lambda_q$, $\lambda_t$, $\lambda_{eff}$ — молекулярная теплопроводность, фазовая теплопроводность, турбулентная теплопроводность и эффективная теплопроводность.

$\mu$, $\mu_q$, $\mu_t$ — динамическая вязкость смеси, фазы и турбулентная вязкость.

$R_g$ — удельная газовая постоянная воздуха.

$K_\ell$ — модуль объёмной упругости жидкости.

$\beta_\ell$ — коэффициент объёмного теплового расширения жидкости.

$\chi_Y$ — коэффициент влияния растворённого газа на плотность жидкости.

$\sigma(T)$ — коэффициент поверхностного натяжения.

$\nu$, $\nu_t$ — кинематическая и турбулентная кинематическая вязкость.

$\nu_{40}$, $\nu_{100}$ — кинематическая вязкость жидкости при 40 °C и 100 °C.

$A$, $B$ — коэффициенты в уравнении Вальтера / ASTM D341 / уравнении вязкость–температура.

$S_\mu$ — константа Сазерленда.

$h_w$ — коэффициент теплоотдачи на внешней границе.

2.2.3. Напряжения, турбулентность и межфазные силы

$I$ — единичный тензор второго ранга.

$\mathbf S$ — тензор скоростей деформации.

$\tau_\mu$ — тензор молекулярных вязких напряжений.

$\tau_{eff}$ — эффективный тензор напряжений.

$k$, $\omega$ — турбулентная кинетическая энергия и удельная скорость диссипации.

$P_k$, $\tilde P_k$ — скорость производства турбулентной кинетической энергии и её ограниченное значение.

$Pr_t$, $Sc_t$ — турбулентные числа Прандтля и Шмидта.

$\mathbf f_\sigma$ — объёмная сила поверхностного натяжения.

$\mathbf n_i$ — единичная нормаль к границе фаз.

$\kappa$ — кривизна границы фаз.

$\Phi_\mu$, $\Phi_t$ — молекулярная и турбулентная диссипация в энергобалансе.

$Q_{\ell g}$ — межфазный тепловой поток между жидкостью и газом.

$a_i$ — удельная межфазная площадь.

$a_1$, $\beta^*$, $\beta_\omega$, $\gamma$, $\sigma_k$, $\sigma_\omega$, $\sigma_{\omega 2}$, $F_1$, $F_2$ — коэффициенты и функции смешения модели SST.

2.2.4. Растворённый газ, динамика поршня и численные показатели

$Y_d$ — массовая доля растворённого воздуха в жидкости.

$D_\ell$ — молекулярный коэффициент диффузии растворённого газа.

$S_d$ — источник массы растворённого компонента.

$k_L$ — коэффициент массоотдачи в жидкости.

$C_d$, $C^*$ — текущая и равновесная концентрации растворённого газа.

$H_s^{cp}(T)$ — константа Генри в форме «концентрация/парциальное давление» ^[27].

$z_p$, $V_p$ — координата и скорость поршня вдоль оси.

$M$ — приведённая масса движущихся частей.

$F_s(z_p)$ — сила упругого элемента.

$F_h(t)$ — гидродинамическая сила на поршне.

$F_{fric,0}(z_p)$ — сухое или конструктивное трение.

$\eta$ — коэффициент механического вязкого сопротивления.

$W_h$ — работа гидродинамического демпфирования.

$Re$, $We$, $Fr$, $Ma$, $Pr$ — числа Рейнольдса, Вебера, Фруда, Маха и Прандтля.

$Co$, $Co_\alpha$ — число Куранта для потока и границы фаз.

$y^+$ — безразмерное пристенное расстояние.

$R_\phi$ — невязка дискретного уравнения для переменной $\phi$.

$E_m$, $E_E$ — относительный дисбаланс массы и энергии.

$p_{obs}$, $Q_{ext}$, $GCI$, $F_{GCI}$ — наблюдаемый порядок сходимости, экстраполированное значение, индекс сеточной сходимости и коэффициент запаса.

$U_{exp}$, $U_{num}$ — экспериментальная и численная неопределённости.

$p'$ — поправка к давлению на текущей итерации.

$\psi=(\partial \rho/\partial p)_{T,\alpha,Y}$ — коэффициент сжимаемости смеси.

3. Цифровой двойник гидроамортизатора

Цифровой двойник в данной задаче не сводится к одной CFD-модели. Это связанная система: физический объект, измерения, высокоточная расчётная модель, идентификация параметров, быстрый вычислительный заменитель и контур принятия инженерных решений. NIST / National Institute of Standards and Technology / Национальный институт стандартов и технологий США трактует цифровой двойник как вычислительное представление физической системы, используемое для прогнозирования, анализа и поддержки решений ^[21]. ISO 23247 задаёт базовый каркас цифрового двойника для производственных систем ^[22], а работа Glaessgen и Stargel связывает цифровой двойник с высокоточным моделированием, обновлением по данным и прогнозированием состояния физического объекта ^[23].

Для гидроамортизатора цифровой двойник целесообразно строить в три уровня.

Уровень A — высокоточная физическая модель.
Это трёхмерный нестационарный расчёт NSF–VOF–SST на деформируемой конечно-объёмной сетке без изменения топологии. Он используется для проектных расчётов, анализа критических режимов, оценки влияния уровня заполнения, угла положения, температуры, газового объёма, профиля канавок и параметров трения.

Уровень B — идентификационная модель.
Она сопоставляет расчёт с измерениями $z_p(t)$, $V_p(t)$, $p(t)$, $T(t)$, положением свободной поверхности и объёмом газа, после чего уточняет параметры $\mu(T,p)$, $h_w$, $\eta$, начальное газосодержание и коэффициенты растворения.

Уровень C — быстрый двойник.
На основе матрицы высокоточных расчётов строится ROM / Reduced-Order Model / редуцированная модель, быстро оценивающая $V_{extr}$, $V_{end}$, $p_{max}$, $T_{max}$, работу демпфирования и эксплуатационные допуски по уровню заполнения и углу установки.

4. Методологический обзор применяемых моделей

4.1. Исходный класс трёхмерных моделей высокоскоростного демпфирования

Исходная постановка задачи высокоскоростного демпфирования в гидроамортизаторе включала трёхмерную конечно-объёмную схему, VOF-описание свободной поверхности, RANS-модель турбулентности и деформируемую сетку, движущуюся вслед за поршнем ^[1]. Сильная сторона этого класса моделей — прямой расчёт гидродинамической силы на поршне по полю давления и вязких напряжений. Именно этот результат затем связывает CFD-часть с ODE-уравнением движения поршня и с цифровым двойником устройства.

4.2. VOF для свободной поверхности, расслоения и вовлечения газа

Классический VOF-метод Хирта и Николса ^[2] остаётся базовым инструментом для описания свободной поверхности, расслоения, перетекания газа и образования крупных газовых карманов. В задачах VOF хорошо описывает: начальное разделение жидкость–газ, деформацию свободной поверхности при разгоне потока, перенос газовых карманов через верхнюю часть зазора и вытеснение жидкости через нижнюю часть.

Для инженерной постановки важен не только сам факт использования VOF, но и выбор схемы переноса границы фаз. Геометрические схемы типа PLIC кусочно-линейной реконструкции границы обладают высокой точностью, однако для сложной движущейся сетки и полностью неявной схемы практичнее использовать алгебраическое отслеживание и фиксацию границ фаз — в частности HRIC ^[5]. В этой работе HRIC принимается как базовый вариант для дискретизации уравнения переноса объёмной доли.

4.3. RANS и выбор базовой модели SST $k$–$\omega$

При $Re=10^5$–$10^6$ прямое разрешение всех турбулентных масштабов в полной трёхмерной геометрии в данной статье для гидроамортизатора самолета практически нецелесообразно. Поэтому инженерный базовый уровень цифрового двойника должен оставаться в классе RANS-моделей. Для пристеночных течений в узких зазорах, канавках и зонах возможного отрыва оптимальным базовым выбором является SST $k$–$\omega$ Ментера ^{[6, 7]}. Она объединяет точность $k$–$\omega$ в пристеночной зоне и устойчивость $k$–$\varepsilon$ в удалённой области потока.

В настоящей статье SST используется как единственное базовое замыкание результирующей системы уравнений. Модель $k$–$\varepsilon$ рассматривается лишь как вспомогательный вариант для анализа чувствительности, но не входит в итоговую систему.

4.4. Деформируемая сетка без изменения топологии и метод IDW

В статье ^[1] использовалась произвольная лагранжево–эйлерова формулировка ^{[4, 5]}. В данном примере для практической реализации цифрового двойника гидроамортизатора самолета удобнее перейти к явной стратегии движения узлов сетки без изменения топологии. При таком подходе:

• топология сетки не перестраивается;
• меняются только координаты узлов и геометрия ячеек;
• консервативность обеспечивается через относительные потоки через движущиеся грани и выполнение GCL;
• качество ячеек поддерживается алгебраической интерполяцией смещений по IDW ^[8].

Для внутренних узлов сетки применяются веса, обратно пропорциональные степени расстояния до граничных узлов. Для узлов на неподвижных поверхностях смещение ограничивается касательной плоскостью, а затем узел проецируется обратно на исходную поверхность. Это позволяет двигать сетку вслед за поршнем, не искажая геометрию кожуха и направляющих.

4.5. Полностью неявная конечно-объёмная схема

Для цифрового двойника важны не только физическая корректность модели, но и численная устойчивость на малых шагах времени и при больших отношениях плотностей и вязкостей. Поэтому результирующая схема формулируется как полностью неявная конечно-объёмная схема на коллокированной сетке ^{[15, 16]}. По времени используется BDF2, конвекция импульса и турбулентных величин аппроксимируется LUD, а граница фаз — HRIC.

Для связи давления и скорости предлагается блочно-сопряжённое решение уравнений импульса и непрерывности. Чтобы избежать давления «в клетку» и сохранить гидростатический баланс, массовые потоки на гранях вычисляются с использованием согласованной импульсно-взвешенной интерполяции типа Rhie–Chow ^{[11, 12]} с поправкой объёмных сил ^{[13, 14]}.

4.6. Различия базового ^[1] и предлагаемого численного ядра

4.7. Альтернативы базовой RANS-постановке

Для большинства инженерных сценариев цифрового двойника RANS остаётся достаточным уровнем детализации. Однако при необходимости исследования тонкой нестационарной структуры течений, распыла и крупномасштабных вихрей в газовом объёме уровень A цифрового двойника может быть расширен до LES-моделирования больших вихрей или гибридных методов DDES задержанного моделирования оторвавшихся вихрей ^{[28, 29]}. Эти модели не входят в базовую результирующую систему уравнений, но полезны как эталонный инструмент для упрощённых геометрий и для анализа чувствительности.

5. Физические явления в гидроамортизаторе со свободным объёмом

В рабочем цикле одновременно действуют несколько групп процессов.

Гидродинамическое дросселирование. Поршень вытесняет жидкость через профилированные зазоры и отверстия. Локальная скорость в зазоре существенно выше средней скорости поршня. Инерционные силы доминируют над вязкими, но вязкость определяет пристенные потери, турбулентную генерацию и тепловой источник.

Газовая подушка и слабая сжимаемость. Свободный газ занимает 1–10 % объёма и сжимается при изменении камер. Даже при малом числе Маха давление газа может сильно меняться из-за изменения объёма. Жидкость рассматривается как слабосжимаемая через модуль объёмной упругости.

Свободная поверхность и расслоение. При ненулевом угле наклона начальное распределение жидкости и газа задаётся гидростатикой. При движении поршня возникают волны, газовые карманы, перетекание газа через верхнюю часть зазора и жидкостные струи через нижнюю часть. Именно эта картина объясняет чувствительность к уровню заполнения и углу.

Турбулентность и распыление. Для оценки режимов полезны безразмерные числа

$$Re=\frac{\rho U L}{\mu},\qquad We=\frac{\rho U^2 L}{\sigma},\qquad Fr=\frac{U}{\sqrt{gL}},\qquad Ma=\frac{U}{c},\qquad Pr=\frac{\mu c_p}{\lambda}.$$

Высокий $Re$ требует турбулентного замыкания. Большой $We$ означает, что инерция способна разрушать поверхность струи быстрее, чем поверхностное натяжение её стабилизирует. Малый $Ma$ позволяет использовать алгоритм на основе давления, но не отменяет изменения давления в газовой подушке.

Тепловая обратная связь. Механическая работа торможения

$$W_h(t)=\int_{z_p(0)}^{z_p(t)} F_h(z)\,dz$$

частично запасается в газе, частично уходит в стенки, а частично превращается в тепло жидкости через диссипацию. Если $\mu_\ell(T)$ падает, уменьшается гидравлическое сопротивление зазора и меняется форма кривой $V_p(z_p)$.

6. Допущения, расчётные параметры и область применимости

Базовая модель принимает следующие допущения.

• Две фазы: жидкость $\ell$ и газ $g$.
• Начальная объёмная доля жидкости составляет $0.90 \le V_\ell/V_0 \le 0.99$, свободного газа — $0.01 \le V_g/V_0 \le 0.10$.
• Течение слабосжимаемое: газ описывается уравнением состояния идеального газа, жидкость — слабосжимаемой моделью с модулем объёмной упругости.
• Межфазная модель — односкоростная VOF: $\mathbf u_\ell=\mathbf u_g=\mathbf u$ на масштабе расчётной ячейки.
• Турбулентность — SST $k$–$\omega$ как базовое RANS-замыкание.
• Испарение, конденсация, кипение и кавитация не включаются. Межфазное взаимодействие в базовом расчёте — только теплообмен. Растворение воздуха в жидкости может быть включено отдельным скаляром.
• Положение цилиндра в пространстве изменяется в диапазоне $\beta\in[0^\circ,90^\circ]$ и входит в задачу через $\mathbf g_{loc}$.
• Геометрия произвольная трёхмерная; на этапе вывода уравнений она задаётся параметрически через область и законы движения границ.

7. Физико-математическая модель

Все обозначения, индексы и англоязычные сокращения, используемые ниже, расшифрованы в разделе 2.

7.1. Деформация расчётной сетки без изменения топологии

Топология сетки не меняется:

$$\mathcal T_h^{n+1}\equiv \mathcal T_h^n.$$

Координаты узлов обновляются по правилу

$$\mathbf x_i^{n+1}=\mathbf x_i^n+\Delta \mathbf x_i.$$

Для узлов на поверхности поршня $\Gamma_p$ задаётся смещение

$$\Delta \mathbf x_b=\Delta z_p\,\mathbf e_z,\qquad b\in\Gamma_p.$$

Для узлов на неподвижных поверхностях $\Gamma_s$ вводится тангенциальное смещение с последующей проекцией на исходную геометрию:

$$\Delta \mathbf x_b= \chi_b\left(I-\mathbf n_b\otimes\mathbf n_b\right)\left(\Delta z_p\,\mathbf e_z\right), \qquad b\in\Gamma_s,$$

$$\mathbf x_b^{n+1}=\Pi_{\Gamma_s}\!\left(\mathbf x_b^n+\Delta \mathbf x_b\right).$$

Тем самым выполняется условие сохранения геометрии неподвижной стенки по нормали:

$$\mathbf n_b\cdot(\mathbf x_b^{n+1}-\mathbf x_b^n)=0.$$

Для внутренних узлов применяется IDW-интерполяция смещений:

$$\Delta \mathbf x_i= \frac{\sum\limits_{b\in\Gamma_p\cup\Gamma_s} w_{ib}\,\Delta \mathbf x_b} {\sum\limits_{b\in\Gamma_p\cup\Gamma_s} w_{ib}}, \qquad w_{ib}=\frac{1}{(d_{ib}+\epsilon_d)^{p_{IDW}}}.$$

Скорость узла и скорость грани вычисляются как

$$\mathbf u_{g,i}^{n+1}=\frac{\mathbf x_i^{n+1}-\mathbf x_i^n}{\Delta t}, \qquad \mathbf u_{g,f}^{n+1}=\mathcal I_f(\mathbf u_{g,i}^{n+1}),$$

где $\mathcal I_f$ — интерполяция на грань.

Эта часть модели не меняет непрерывные уравнения в частных производных в физической области, а определяет закон движения сетки и геометрию деформируемых контрольных объёмов.

7.2. Уравнения состояния и свойства фаз

Объёмные доли удовлетворяют условию

$$\alpha_\ell+\alpha_g=1,\qquad 0\le \alpha_q\le 1.$$

Смесевые величины определяются как

$$\rho=\alpha_\ell\rho_\ell+\alpha_g\rho_g, \qquad \rho h=\alpha_\ell \rho_\ell h_\ell+\alpha_g\rho_g h_g.$$

Для газа используется EOS / Equation of State / уравнение состояния идеального газа с абсолютным давлением:

$$\rho_g=\frac{p_{abs}}{R_g T}, \qquad p_{abs}=p_{ref}+p.$$

Для жидкости принимается слабосжимаемая параметрическая модель:

$$\rho_\ell= \rho_{\ell0} \left[ 1+\frac{p_{abs}-p_0}{K_\ell} -\beta_\ell (T-T_0) +\chi_Y (Y_d-Y_{d0}) \right].$$

Если перенос растворённого газа не рассматривается, член с $Y_d$ опускается.

Динамическая вязкость воздуха описывается законом Сазерленда ^{[17, 18]}:

$$\mu_g(T)=\mu_{g0}\left(\frac{T}{T_0}\right)^{3/2}\frac{T_0+S_\mu}{T+S_\mu}.$$

Для гидравлической жидкости при наличии двух паспортных точек $\nu_{40}$ и $\nu_{100}$ используется корреляция ASTM D341 ^[20]:

$$\log_{10}\!\left[\log_{10}(\nu+0.7)\right]=A-B\log_{10}T, \qquad \mu_\ell=\rho_\ell \nu \cdot 10^{-6},$$

где $\nu$ задаётся в cSt, а $T$ — в K. Для водных сред может использоваться формуляция IAPWS ^[19]; для рабочей жидкости цифрового двойника предпочтительна табличная функция $\mu_\ell(T,p)$, построенная по паспортным и стендовым данным.

Теплопроводность и теплоёмкость смеси задаются как

$$\lambda=\alpha_\ell\lambda_\ell+\alpha_g\lambda_g, \qquad \rho c_p=\alpha_\ell \rho_\ell c_{p,\ell}+\alpha_g\rho_g c_{p,g}.$$

Эффективная теплопроводность

$$\lambda_{eff}=\lambda+\lambda_t, \qquad \lambda_t=\frac{\mu_t c_p}{Pr_t}.$$

7.3. Баланс фазовой массы и VOF

Математически корректная форма для фазовой массы имеет вид

$$\frac{\partial(\alpha_q\rho_q)}{\partial t} +\nabla\cdot(\alpha_q\rho_q\mathbf u) =\dot m_q, \qquad q\in\{\ell,g\}.$$

В базовой постановке без испарения, конденсации и кавитации

$$\dot m_\ell=\dot m_g=0.$$

Отсюда видно, что уравнение для объёмной доли в виде

$$\frac{\partial \alpha_q}{\partial t}+\nabla\cdot(\alpha_q\mathbf u)=0$$

является частным случаем только при постоянной $\rho_q$ и отсутствии объёмного расширения. Именно поэтому результирующая модель строится на фазовой массе, а не на упрощённом транспорте $\alpha_q$.

Для сохранения резкой границы фаз в численной схеме допускается добавление чисто численного сжимающего потока, связанного с HRIC. Этот поток не является физическим массообменом между фазами и не должен нарушать ограниченность объёмной доли:

$$0\le \alpha_\ell \le 1.$$

7.4. Баланс массы смеси

Суммирование фазовых уравнений даёт уравнение неразрывности смеси:

$$\frac{\partial \rho}{\partial t}+\nabla\cdot(\rho \mathbf u)=0.$$

В слабосжимаемой постановке дивергенция скорости не обращается тождественно в нуль, а определяется изменением плотности:

$$\nabla\cdot \mathbf u = -\frac{1}{\rho}\frac{D\rho}{Dt}.$$

Линеаризация плотности для построения уравнения давления записывается как

$$d\rho= \left(\frac{\partial \rho}{\partial p}\right)_{T,\alpha,Y}dp + \left(\frac{\partial \rho}{\partial T}\right)_{p,\alpha,Y}dT + \left(\frac{\partial \rho}{\partial \alpha}\right)_{p,T,Y}d\alpha + \left(\frac{\partial \rho}{\partial Y_d}\right)_{p,T,\alpha}dY_d.$$

7.5. Баланс импульса

Уравнение импульса смеси в физической области:

$$\frac{\partial (\rho \mathbf u)}{\partial t} +\nabla\cdot(\rho\mathbf u\otimes \mathbf u) = -\nabla p +\nabla\cdot \tau_{eff} +\rho \mathbf g_{loc} +\mathbf f_\sigma.$$

Здесь

$$\tau_{eff} = 2(\mu+\mu_t)\mathbf S -\frac{2}{3}(\mu+\mu_t)(\nabla\cdot\mathbf u)I.$$

Если используется полная RANS-запись с явной изотропной турбулентной частью, то

$$\tau_{eff}^{RANS} = 2(\mu+\mu_t)\mathbf S -\frac{2}{3}(\mu+\mu_t)(\nabla\cdot\mathbf u)I -\frac{2}{3}\rho k I,$$

а давление должно быть переопределено последовательно во всех формулах.

Поверхностное натяжение в CSF-форме записывается как

$$\mathbf f_\sigma=\sigma(T)\,\kappa\,\nabla \alpha_\ell, \qquad \kappa=-\nabla\cdot\mathbf n_i, \qquad \mathbf n_i=\frac{\nabla \alpha_\ell}{|\nabla\alpha_\ell|+\epsilon}.$$

При необходимости термокапиллярный вклад может быть добавлен как дополнительный член, пропорциональный $\nabla \sigma(T)$, но в базовом расчёте он обычно вторичен по сравнению с инерцией и перепадом давления.

7.6. Уравнение энергии Навье–Стокса–Фурье

В однотемпературной LTE-форме уравнение энергии удобно записывать через энтальпию:

$$\frac{\partial(\rho h)}{\partial t} +\nabla\cdot(\rho \mathbf u h) = \nabla\cdot(\lambda_{eff}\nabla T) +\frac{D p_{abs}}{Dt} +\Phi_\mu+\Phi_t+S_h.$$

Молекулярная диссипация

$$\Phi_\mu=\tau_\mu:\nabla \mathbf u, \qquad \tau_\mu= 2\mu \mathbf S -\frac{2}{3}\mu (\nabla\cdot\mathbf u)I.$$

Для SST-модели турбулентный тепловой источник задаётся как

$$\Phi_t=\beta^* \rho k \omega.$$

Если требуется раздельный теплообмен между газом и жидкостью, система может быть расширена до двухтемпературной модели:

$$\frac{\partial(\alpha_q\rho_q h_q)}{\partial t} +\nabla\cdot(\alpha_q\rho_q \mathbf u h_q) = \nabla\cdot\!\left[\alpha_q(\lambda_q+\lambda_{t,q})\nabla T_q\right] +\alpha_q \Phi_q +Q_{qr},$$

$$Q_{\ell g}=h_{\ell g}\,a_i\,(T_g-T_\ell), \qquad Q_{g\ell}=-Q_{\ell g}.$$

В базовой результирующей модели, чтобы не усложнять цифровой двойник на первом этапе, используется однотемпературная LTE-постановка.

7.7. Замыкание турбулентности: SST $k$–$\omega$

Результирующая модель замыкается уравнениями SST $k$–$\omega$ ^{[6, 7]}:

$$\frac{\partial(\rho k)}{\partial t} +\nabla\cdot(\rho \mathbf u k) = \tilde P_k-\beta^*\rho \omega k +\nabla\cdot\left[(\mu+\sigma_k \mu_t)\nabla k\right],$$

$$\frac{\partial(\rho \omega)}{\partial t} +\nabla\cdot(\rho \mathbf u \omega) = \gamma\frac{\rho}{\mu_t}\tilde P_k -\beta_\omega \rho \omega^2 +\nabla\cdot\left[(\mu+\sigma_\omega \mu_t)\nabla \omega\right] + 2(1-F_1)\frac{\rho\sigma_{\omega 2}}{\omega}\nabla k\cdot\nabla \omega.$$

Турбулентная вязкость вычисляется как

$$\mu_t=\frac{\rho a_1 k}{\max(a_1\omega,\;|S|F_2)}, \qquad |S|=\sqrt{2S_{ij}S_{ij}}.$$

Ограниченное производство турбулентной энергии:

$$\tilde P_k = \min\!\left(P_k,\;10\beta^*\rho\omega k\right).$$

Эта модель особенно полезна для гидроамортизатора, потому что корректно описывает пристеночные напряжения и градиенты давления в профилированных дросселирующих зазорах.

7.8. Опциональный перенос растворённого воздуха

Если требуется описывать растворение воздуха в жидкости, вводится дополнительный скаляр:

$$\frac{\partial(\alpha_\ell \rho_\ell Y_d)}{\partial t} +\nabla\cdot(\alpha_\ell \rho_\ell \mathbf u Y_d) = \nabla\cdot(\alpha_\ell \rho_\ell D_\ell \nabla Y_d)+S_d.$$

Один из допустимых вариантов источника:

$$S_d = k_L a_i (C^*-C_d), \qquad C^*=H_s^{cp}(T)\,p_a.$$

Эта подсистема не является обязательной для базового расчёта, но полезна при моделировании длительных циклических режимов, когда газосодержание жидкости меняется от цикла к циклу.

7.9. Уравнение движения поршня и связь с гидродинамикой

Связь между CFD-моделью и механикой устройства задаётся системой

$$\frac{dz_p}{dt}=V_p,$$

$$M\frac{dV_p}{dt} = F_{ext}(t)-F_s(z_p)-F_{fric,0}(z_p)-\eta V_p-F_h(t).$$

Гидродинамическая сила определяется интегрированием напряжений по поверхности поршня:

$$F_h(t)= \int_{S_p} \left[ (-pI+\tau_{eff})\cdot\mathbf n \right]\cdot\mathbf e_z\,dS.$$

Именно эта сила формирует тормозное действие устройства, а через работу

$$W_h=\int F_h\,dz_p$$

связывает механическую и тепловую подсистемы.

8. Начальные и граничные условия

Начальные условия должны быть гидростатически согласованы с заданным углом положения цилиндра. В начальный момент задаются:

• распределение объёмной доли $\alpha_\ell(\mathbf x,0)$, соответствующее уровню заполнения;
• гидростатическое поле давления:

  $$\nabla p(\mathbf x,0)=\rho(\mathbf x,0)\,\mathbf g_{loc};$$

• начальная температура $T_0$;
• начальное положение и скорость поршня $z_p(0)$, $V_p(0)$;
• начальные поля $k$ и $\omega$.

Для стенок и поршня:

$$\mathbf u=\mathbf u_w, \qquad \mathbf u_w= \begin{cases} 0, & \text{неподвижная стенка},\\ V_p \mathbf e_z, & \text{поршень}. \end{cases}$$

Для движения сетки:

$$\Delta \mathbf x_b=\Delta z_p \mathbf e_z,\qquad b\in\Gamma_p,$$

$$\Delta \mathbf x_b= \chi_b(I-\mathbf n_b\otimes\mathbf n_b)(\Delta z_p \mathbf e_z),\qquad b\in\Gamma_s.$$

Тепловые условия выбираются по назначению расчёта:

$$\nabla T\cdot \mathbf n = 0 \quad \text{адиабатическая стенка},$$

$$T=T_w \quad \text{изотермическая стенка},$$

$$-\lambda_{eff}\nabla T\cdot\mathbf n = h_w (T-T_{ext}) \quad \text{конвективная стенка}.$$

Для VOF на стенке задаётся условие непроникновения; если поверхностное натяжение существенно, дополнительно задаётся контактный угол. Для SST возможны два режима:

• низкорейнольдсовая интеграция, когда первая ячейка лежит в вязком подслое и

  $$k_w=0, \qquad \omega_w \approx \frac{60\nu}{\beta_1 (\Delta y_1)^2};$$

• автоматическая комбинированная пристенная обработка, когда сетка неоднородна и часть поверхности может иметь $y^+>1$.

Для гидроамортизатора предпочтителен гибридный automatic wall treatment / автоматический комбинированный пристенный режим: в зазоре целевой $y^+\sim 1$, на удалённых стенках допустим $y^+\sim 30$–$100$.

9. Численная схема решения

9.1. Интегральная конечно-объёмная форма на деформируемом контрольном объёме

Для произвольной переносимой величины $\phi$ консервативная конечно-объёмная запись на деформируемом контрольном объёме $V_P(t)$ имеет вид

$$\frac{d}{dt}\int_{V_P(t)} \rho \phi\,dV + \sum_{f\in \partial V_P} \dot m_f^{rel}\,\phi_f = \sum_{f\in \partial V_P} \Gamma_{\phi,f}\,(\nabla \phi)_f\cdot \mathbf S_f + \int_{V_P(t)} S_\phi\,dV,$$

где относительный массовый поток через грань

$$\dot m_f^{rel} = \rho_f\left[(\mathbf u_f-\mathbf u_{g,f})\cdot\mathbf S_f\right].$$

Таким образом, движение сетки учитывается только через скорости граней $\mathbf u_{g,f}$ и изменение геометрии ячеек. Этот подход эквивалентен по физическому смыслу использованию движущихся контрольных объёмов, но удобнее для алгоритма цифрового двойника.

9.2. Трёхслойная схема BDF2 и геометрический закон сохранения

По времени для величины $\rho\phi V$ используется трёхслойная схема второго порядка:

$$\left. \frac{d}{dt}(\rho \phi V)_P \right|^{n+1} \approx \frac{ 3(\rho\phi V)_P^{n+1} - 4(\rho\phi V)_P^n + (\rho\phi V)_P^{n-1} }{ 2\Delta t }.$$

Для деформируемой сетки дополнительно должно выполняться дискретное GCL:

$$\frac{ 3V_P^{n+1} - 4V_P^n + V_P^{n-1} }{ 2\Delta t } = \sum_{f\in \partial V_P} \phi_{g,f}^{n+1}, \qquad \phi_{g,f}^{n+1}=\mathbf u_{g,f}^{n+1}\cdot\mathbf S_f^{n+1}.$$

Выполнение этого условия является обязательным тестом корректности движения сетки: неподвижное однородное поле не должно создавать фиктивный массовый поток только из-за деформации ячеек ^{[9, 10]}.

9.3. Дискретизация VOF-уравнения схемой HRIC

Для жидкой фазы на ячейке $P$:

$$\frac{ 3(\alpha_\ell \rho_\ell V)_P^{n+1} - 4(\alpha_\ell \rho_\ell V)_P^{n} + (\alpha_\ell \rho_\ell V)_P^{n-1} }{ 2\Delta t } + \sum_f \alpha_{\ell,f}^{HRIC}\rho_{\ell,f}^{n+1}\phi_{rel,f}^{n+1} =0,$$

где

$$\phi_{rel,f}^{n+1}=(\mathbf u_f^{n+1}-\mathbf u_{g,f}^{n+1})\cdot\mathbf S_f^{n+1}.$$

Значение объёмной доли на грани строится как

$$\alpha_{\ell,f}^{HRIC} = \mathcal H_{HRIC} \left( \alpha_{\ell,U}, \alpha_{\ell,D}, \nabla\alpha_\ell, \theta_f, Co_f \right),$$

где $U$ и $D$ — ячейки по относительному потоку, $\theta_f$ — угол между нормалью границы фаз и нормалью грани, а $Co_f$ — локальное число Куранта. В HRIC компрессивная ветвь автоматически ослабляется при неблагоприятном $Co_f$, что повышает устойчивость и сохраняет ограниченность:

$$0\le \alpha_\ell^{n+1}\le 1.$$

9.4. Дискретизация уравнения импульса схемой LUD

Для каждой компоненты скорости $u_r$, $r\in\{x,y,z\}$:

$$\frac{ 3(\rho u_r V)_P^{n+1} - 4(\rho u_r V)_P^n + (\rho u_r V)_P^{n-1} }{ 2\Delta t } + \sum_f \dot m_f^{rel,n+1}\,(u_r)_f^{LUD} = -\sum_f p_f^{n+1}S_{f,r} + \sum_f (\tau_{eff}\cdot\mathbf S_f)_r^{n+1} + (\rho g_r V)_P^{n+1} + (f_{\sigma,r}V)_P^{n+1}.$$

LUD-реконструкция значения на грани записывается как

$$\phi_f^{LUD} = \phi_U + (\nabla \phi)_U\cdot(\mathbf x_f-\mathbf x_U).$$

На гладком решении эта схема даёт второй порядок по пространству, а при необходимости ограничитель может применяться локально для предотвращения неограниченных осцилляций.

9.5. Дискретизация энергии

Энтальпийное уравнение дискретизируется аналогично:

$$\frac{ 3(\rho h V)_P^{n+1} - 4(\rho h V)_P^n + (\rho h V)_P^{n-1} }{ 2\Delta t } + \sum_f \dot m_f^{rel,n+1} h_f^{LUD} = \sum_f \lambda_{eff,f}^{n+1}(\nabla T)_f^{n+1}\cdot \mathbf S_f^{n+1} + \left( \frac{D p_{abs}}{Dt} +\Phi_\mu+\Phi_t+S_h \right)_P^{n+1}V_P^{n+1}.$$

9.6. Дискретизация SST-уравнений

Для $k$:

$$\frac{ 3(\rho k V)_P^{n+1} - 4(\rho k V)_P^n + (\rho k V)_P^{n-1} }{ 2\Delta t } + \sum_f \dot m_f^{rel,n+1} k_f^{LUD} = \sum_f \Gamma_{k,f}^{n+1}(\nabla k)_f^{n+1}\cdot\mathbf S_f^{n+1} + (\tilde P_k-\beta^*\rho\omega k)_P^{n+1}V_P^{n+1}.$$

Для $\omega$:

$$\frac{ 3(\rho \omega V)_P^{n+1} - 4(\rho \omega V)_P^n + (\rho \omega V)_P^{n-1} }{ 2\Delta t } + \sum_f \dot m_f^{rel,n+1} \omega_f^{LUD} = \sum_f \Gamma_{\omega,f}^{n+1}(\nabla \omega)_f^{n+1}\cdot\mathbf S_f^{n+1} + \left( \gamma\frac{\rho}{\mu_t}\tilde P_k -\beta_\omega\rho\omega^2 +D_\omega \right)_P^{n+1}V_P^{n+1}.$$

9.7. Полностью неявное решение подсистемы Навье–Стокса

На каждом временном слое система для $(u,v,w,p)$ собирается в блочно-сопряжённом виде:

$$\begin{bmatrix} A_u & 0 & 0 & G_x \\ 0 & A_v & 0 & G_y \\ 0 & 0 & A_w & G_z \\ D_x & D_y & D_z & C_p \end{bmatrix} \begin{bmatrix} u\\v\\w\\p \end{bmatrix}^{n+1} = \begin{bmatrix} b_u\\b_v\\b_w\\b_c \end{bmatrix}.$$

Здесь блоки $A_u,A_v,A_w$ содержат вклад BDF2, LUD-конвекции, диффузии и неявных источников, а блоки $G_r$ и $D_r$ отвечают за связь давления и скорости. Член $C_p$ содержит вклад сжимаемости смеси.

Для вычисления массового потока через грань применяется согласованная интерполяция с поправкой объёмных сил:

$$\dot m_f^{n+1} = \rho_f^{n+1} \left[ \mathbf u_{H,f}^{n+1} - d_f \left( (\nabla p)_f^{n+1} - \rho_f^{n+1}\mathbf g_{loc} \right) \right]\cdot\mathbf S_f^{n+1} - \rho_f^{n+1}\phi_{g,f}^{n+1}.$$

Здесь $\mathbf u_{H,f}$ — скорость на грани, восстановленная из неявного оператора импульса без давления, а $d_f$ — коэффициент, связанный с диагональю дискретизированного уравнения импульса. Такая формула важна по двум причинам:

1. устраняется расщепление давления и скорости на коллокированной сетке ^{[11, 12]};
2. сохраняется правильный гидростатический баланс: если $(\nabla p)_f=\rho_f\mathbf g_{loc}$, фиктивный массовый поток не возникает ^{[13, 14]}.

Непрерывность смеси записывается как

$$\frac{ 3(\rho V)_P^{n+1} - 4(\rho V)_P^n + (\rho V)_P^{n-1} }{ 2\Delta t } + \sum_f \dot m_f^{n+1}=0.$$

При линеаризации плотности

$$\rho^{n+1}=\rho^*+\psi p'$$

получается уравнение поправки давления:

$$\frac{3}{2\Delta t}\psi_P V_P\,p'_P - \sum_f \rho_f D_{u,f}(\nabla p')_f\cdot\mathbf S_f = -R_{m,P},$$

где $R_{m,P}$ — невязка дискретной непрерывности по предварительным потокам. В зависимости от линейного решателя цифровой двойник может использовать либо явное уравнение поправки давления, либо прямое блочно-сопряжённое решение без выделения $p'$ как отдельной стадии.

9.8. Внешние нелинейные итерации и адаптация шага времени

Полная схема нелинейна из-за зависимости свойств от температуры и давления, HRIC-реконструкции, турбулентной вязкости, движения поршня и сетки. Поэтому на каждом временном слое выполняются внешние итерации Пикара–Ньютона.

Рекомендуемые ограничения шага:

$$Co=\max_f\left|\frac{(\mathbf u_f-\mathbf u_{g,f})\cdot\mathbf S_f}{V_P}\Delta t\right|\le 0.5, \qquad Co_\alpha \le 0.25\text{–}0.5.$$

Дополнительно контролируются:

• положительность объёмов ячеек $V_P>0$;
• рост невязок на внешних итерациях;
• изменение температуры и вязкости за один шаг;
• качество деформированной сетки: / перекос / неортогональность / отношение сторон.

10. Алгоритм численного решения системы уравнений

1. Переход на новый временной слой.
   Известны $\alpha^n$, $p^n$, $\mathbf u^n$, $T^n$, $k^n$, $\omega^n$, $z_p^n$, $V_p^n$ и состояние слоя $n-1$ для BDF2.

2. Прогноз хода поршня.
   По ODE уравнению движения поршня вычисляется предварительное значение $z_p^{n+1,0}$, $V_p^{n+1,0}$.

3. Деформация сетки.
   Узлы на поршне смещаются на $\Delta z_p$. Узлы на неподвижных стенках перемещаются тангенциально и проектируются на исходную геометрию. Внутренние узлы перемещаются методом IDW. Пересчитываются объёмы ячеек $V_P^{n+1}$, площади граней $\mathbf S_f^{n+1}$ и скорости граней $\mathbf u_{g,f}^{n+1}$.

4. Проверка сетки и GCL.
   Если объёмы ячеек отрицательны или GCL не выполняется в требуемой точности, шаг времени уменьшается и деформация повторяется.

5. Внешняя нелинейная итерация $m=0,1,\dots$.

6. Обновление свойств среды.
   По полям $T^{m}$, $p^{m}$, $\alpha^{m}$ вычисляются $\rho_q$, $\mu_q$, $\lambda_q$, $c_{p,q}$, затем смешанные свойства $\rho$, $\mu$, $\lambda$, $\lambda_{eff}$.

7. Неявное решение фазовой массы.
   Решается VOF-уравнение в фазово-массовой форме с HRIC для $\alpha_f$. После решения обеспечивается ограниченность $0\le\alpha_\ell\le1$.

8. Сборка и решение блочно-сопряжённой подсистемы $(u,v,w,p)$.
   Конвекция — LUD, потоки на гранях — с поправкой объёмных сил и скоростью сетки.

9. Неявное решение уравнения энергии.
   Решается энтальпийное уравнение с обновлёнными свойствами и диссипацией.

10. Неявное решение SST $k$–$\omega$.
    Решаются два уравнения турбулентности, обновляется $\mu_t$.

11. Опциональный расчёт растворённого воздуха.
    Если используется модель $Y_d$, решается транспортное уравнение растворённого газа.

12. Интегрирование силы на поршне и коррекция ODE.
    По обновлённым полям вычисляется $F_h^{m+1}$ и корректируются $z_p^{m+1}$, $V_p^{m+1}$.

13. Проверка сходимости внешней итерации.
    Контролируются невязки уравнений, дисбаланс массы и энергии, изменение $F_h$, $z_p$ и максимальной температуры.

14. Принятие временного слоя.
    Если критерии выполнены, значения переносятся на следующий слой; в противном случае продолжаются внешние итерации либо уменьшается шаг времени.

11. Результирующая математическая модель

Искомый вектор базовой модели имеет вид

$$\mathcal U= \{\alpha_\ell,\ p,\ \mathbf u,\ T,\ k,\ \omega,\ z_p,\ V_p\},$$

а при необходимости расширяется скаляром $Y_d$.

Результирующая система состоит из:

$$\boxed{ \begin{aligned} &\alpha_\ell+\alpha_g=1,\\[2pt] &\frac{\partial(\alpha_q\rho_q)}{\partial t} +\nabla\cdot(\alpha_q\rho_q\mathbf u)=0,\qquad q\in\{\ell,g\},\\[2pt] &\frac{\partial \rho}{\partial t}+\nabla\cdot(\rho \mathbf u)=0,\\[2pt] &\frac{\partial (\rho \mathbf u)}{\partial t} +\nabla\cdot(\rho \mathbf u\otimes \mathbf u) = -\nabla p+\nabla\cdot\tau_{eff} +\rho \mathbf g_{loc}+\mathbf f_\sigma,\\[2pt] &\frac{\partial(\rho h)}{\partial t} +\nabla\cdot(\rho \mathbf u h) = \nabla\cdot(\lambda_{eff}\nabla T) +\frac{D p_{abs}}{Dt} +\Phi_\mu+\Phi_t+S_h,\\[2pt] &\frac{\partial(\rho k)}{\partial t} +\nabla\cdot(\rho \mathbf u k) = \tilde P_k-\beta^*\rho \omega k +\nabla\cdot[(\mu+\sigma_k\mu_t)\nabla k],\\[2pt] &\frac{\partial(\rho \omega)}{\partial t} +\nabla\cdot(\rho \mathbf u \omega) = \gamma\frac{\rho}{\mu_t}\tilde P_k -\beta_\omega\rho \omega^2 +\nabla\cdot[(\mu+\sigma_\omega\mu_t)\nabla\omega] +D_\omega,\\[2pt] &\rho_g=\frac{p_{abs}}{R_g T}, \qquad \rho_\ell=\rho_\ell(p_{abs},T,Y_d),\\[2pt] &\mu_g=\mu_g(T), \qquad \mu_\ell=\mu_\ell(T,p),\\[2pt] &\frac{dz_p}{dt}=V_p, \qquad M\frac{dV_p}{dt}=F_{ext}-F_s-F_{fric,0}-\eta V_p-F_h,\\[2pt] &\mathbf x_i^{n+1}=\mathbf x_i^n+\Delta \mathbf x_i, \qquad \Delta \mathbf x_i=\mathrm{IDW}(\Delta \mathbf x_b),\\[2pt] &\mathcal T_h^{n+1}\equiv\mathcal T_h^n. \end{aligned} }$$

Эта система дополняется начальными и граничными условиями, законом движения сетки, требованиями GCL и полностью неявной дискретизацией конечными объёмами на деформируемой сетке.

12. Метод проверки сходимости решения системы уравнений

12.1. Алгебраическая сходимость на временном слое

Для каждой дискретной переменной $\phi$ контролируется относительная невязка:

$$R_\phi^{(m)}= \frac{\|A_\phi \phi^{(m)}-b_\phi\|_1}{\|b_\phi\|_1+\epsilon_R}.$$

Рекомендуемые критерии принятия временного слоя:

• $R_p < 10^{-6}$;
• $R_{\alpha} < 10^{-6}$;
• $R_u,\ R_T,\ R_k,\ R_\omega < 10^{-5}$.

Для сильносвязанных режимов более информативны также интегральные дисбалансы массы и энергии.

12.2. Контроль баланса массы и энергии

Для полной массы смеси:

$$E_m= \frac{ \left| M^{n+1}-M^n+\int_{t_n}^{t_{n+1}}\sum_{\partial\Omega}\dot m\,dt \right| }{M^n+\epsilon_E}.$$

Для замкнутого объёма без фазового массопереноса и без утечек через границы должно выполняться

$$E_m \ll 1.$$

Энергетический дисбаланс оценивается как

$$E_E= \frac{ \left| \Delta E_{fluid} + Q_w - W_h \right| }{ |W_h|+\epsilon_E }.$$

12.3. Сеточная и временная сходимость

Для трёх сеток с отношением сгущения $r=h_2/h_1=h_3/h_2$ наблюдаемый порядок сходимости вычисляется как

$$p_{obs} = \frac{\ln\left|\dfrac{Q_3-Q_2}{Q_2-Q_1}\right|}{\ln r}.$$

Экстраполированное значение:

$$Q_{ext}=Q_1+\frac{Q_1-Q_2}{r^{p_{obs}}-1}.$$

Индекс сеточной сходимости:

$$GCI_{12} = F_{GCI} \frac{|Q_1-Q_2|}{|Q_1|(r^{p_{obs}}-1)}, \qquad F_{GCI}=1.25$$

для асимптотической области по Roache ^[25]. При неасимптотическом поведении используется более консервативный запас ^[26].

Для времени аналогичная процедура выполняется на последовательности шагов $\Delta t_1>\Delta t_2>\Delta t_3$.

12.4. Специальные тесты для деформируемой сетки

Для новой численной схемы необходимо выполнить три тестовые задачи.

1. GCL-тест.
Однородное неподвижное поле на движущейся сетке не должно создавать фиктивной массы или импульса.

2. Гидростатический тест на наклонённом цилиндре.
При $\mathbf u=0$ и $\nabla p=\rho \mathbf g_{loc}$ схема должна сохранять гидростатику без паразитных течений.

3. Тест замкнутого объёма с заданным ходом поршня.
При отсутствии межфазного массопереноса и утечек изменение полной массы в области допускается только в пределах машинной/дискретизационной погрешности.

12.5. Контроль границы фаз

Для VOF-решения необходимо контролировать:

• потерю фазовой массы;
• изменение объёма газовой камеры;
• толщину размытия границы фаз;
• величину паразитных скоростей в квазистатическом режиме.

13. Верификация и валидация численного эксперимента

13.1. Верификация кода

Для кода цифрового двойника выполним программу верификации кода:

1. MMS для однофазной weakly compressible NSF-подсистемы с переменной вязкостью.
2. Куэтт–Пуазейль с вязким нагревом как тест энергии и температурно-зависимой вязкости.
3. Скачок Лапласа для проверки CSF-реализации ^[3].
4. GCL-тест на деформируемой IDW-сетке ^{[9, 10]}.
5. Гидростатический тест для BFW-интерполяции потоков ^{[13, 14]}.
6. Тесты SST на эталонных пристеночных течениях по формулировке Ментера ^{[6, 7]}.

13.2. Верификация решения

Для каждого ключевого режима необходимо выполнить:

• расчёты минимум на трёх сетках;
• расчёты минимум на трёх временных шагах;
• контроль $E_m$, $E_E$ и GCI;
• сравнение интегральных выходных величин $F_h(t)$, $p_{max}$, $T_{max}$, $V_{end}$;
• контроль устойчивости решения к изменению параметров движения сетки ($p_{IDW}$, радиус влияния, коэффициент $\chi_b$).

13.3. Валидационный эксперимент

Минимальный набор измеряемых величин на стенде:

• перемещение и скорость поршня $z_p(t)$, $V_p(t)$;
• давление в камерах;
• температура жидкости вблизи зазора и температура стенки;
• начальный уровень жидкости и объём свободного газа;
• угол установки $\beta$;
• свойства жидкости до и после испытания: $\mu(T)$, плотность, газосодержание.

Валидационная метрика для целевой величины $Q$:

$$\delta_Q= \frac{|Q_{CFD}-Q_{exp}|} {U_{exp}+U_{num}+|Q_{exp}|\epsilon_E}.$$

Согласно ASME V&V 20 ^[24], численная модель считается согласованной с экспериментом, если расхождение не превышает суммарной неопределённости, а не просто «визуально похоже».

13.4. Роль валидации в цифровом двойнике

Верификация и валидация это центральный элемент цифрового двойника: без оценённой численной и экспериментальной неопределённости невозможно надёжно перенести модель с уровня A на уровень C. Именно поэтому в цифровом двойнике должны храниться не только расчётные кривые, но и метаданные об использованной сетке, временном шаге, модели турбулентности и статусе V&V.

14. Выводы

В статье получена физико-математическая модель цифрового двойника гидроамортизатора высокоскоростного демпфирования, в которой:

• исходная изотермическая постановка расширена до системы Навье–Стокса–Фурье с теплопереносом и температурно-зависимой вязкостью;
• межфазная подсистема записана в фазово-массовой форме VOF, корректной для слабосжимаемого газа и жидкости;
• турбулентность замкнута моделью SST $k$–$\omega$;
• движение поршня учтено через деформацию конечно-объёмной сетки без изменения её топологии;
• внутренние узлы перемещаются методом IDW, а узлы на неподвижных поверхностях — тангенциально по поверхности;
• предложен полностью неявный алгоритм интегрирования на основе BDF2, LUD, HRIC и поправки объёмных сил;
• сформулированы критерии сходимости, сеточной проверки и программа верификации и валидации.

Таким образом, предложенная модель пригодна как для прямого численного эксперимента, так и для построения инженерного цифрового двойника, связанного с измерениями, идентификацией параметров и расчётной моделью.

Использование цифрового двойника в проектировании позволяет:

• смоделировать и определить физические параметры работы устройства (скорость, давление, температуру и др.) в различных условиях;
• рассчитать и выбрать оптимальные конструктивные параметры устройства и свойства материалов;
• рассчитать условия эксплуатации и технического обслуживания разработанного устройства;
• визуализировать в динамике работу устройства при различных условиях и нагрузках;
• сократить сроки и стоимость разработки, исключив изготовление и испытание полномасштабных образцов и многочисленные натурные эксперименты;
• определить и выбрать оптимальные экономические параметры: стоимость разработки и производства, стоимость владения и стоимость жизненного цикла устройства.

Этим материалом замыкается первый цикл из четырёх статей: от вычислительной базы CFD и HPC — к инженерной рамке цифровых двойников, затем к иерархической декомпозиции изделия и, наконец, к детальной предметной постановке цифрового двойника гидроамортизатора. Но содержательно это не конец темы, а точка перехода: дальше логично рассматривать, как подобные подходы масштабируются на другие физические процессы, изделия, технологические процессы и отрасли, как меняются вычислительные практики под влиянием супервычислений и как AI / искусственный интеллект начинает работать не только как инструмент автоматизации, но и как слой анализа, идентификации параметров, ROM / моделей пониженного порядка и поддержки инженерных решений.

Список источников

1. Efremov V., Kozelkov A., Dmitriev S., Kurkin A., Kurulin V., Utkin D. Technology of 3D Simulation of High-Speed Damping Processes in the Hydraulic Brake Device. IntechOpen, 2018. doi.org/…
2. Hirt C.W., Nichols B.D. Volume of Fluid (VOF) Method for the Dynamics of Free Boundaries. Journal of Computational Physics, 39(1), 201–225, 1981. doi.org/…
3. Brackbill J.U., Kothe D.B., Zemach C. A Continuum Method for Modeling Surface Tension. Journal of Computational Physics, 100(2), 335–354, 1992. doi.org/…
4. Muzaferija S., Perić M. Computation of Free-Surface Flows Using the Finite-Volume Method and Moving Grids. Numerical Heat Transfer, Part B, 32(4), 369–384, 1997. doi.org/…
5. Muzaferija S., Perić M. Computation of Free-Surface Flows Using Interface-Tracking and Interface-Capturing Methods. In: Mahrenholtz O., Markiewicz M. (eds.), Nonlinear Water Wave Interaction. WIT Press / Computational Mechanics Publications, 1999, pp. 59–100.
6. Menter F.R. Two-Equation Eddy-Viscosity Turbulence Models for Engineering Applications. AIAA Journal, 32(8), 1598–1605, 1994. doi.org/…
7. Menter F.R., Kuntz M., Langtry R. Ten Years of Industrial Experience with the SST Turbulence Model. In: Hanjalić K., Nagano Y., Tummers M. (eds.), Turbulence, Heat and Mass Transfer 4. Begell House, 2003, pp. 625–632.
8. Witteveen J.A.S., Bijl H. Explicit Mesh Deformation Using Inverse Distance Weighting Interpolation. Proceedings of the 19th AIAA Computational Fluid Dynamics Conference, AIAA-2009-3996, 2009. doi.org/…
9. Thomas P.D., Lombard C.K. Geometric Conservation Law and Its Application to Flow Computations on Moving Grids. AIAA Journal, 17(10), 1030–1037, 1979. doi.org/…
10. Demirdžić I., Perić M. Space Conservation Law in Finite Volume Calculations of Fluid Flow. International Journal for Numerical Methods in Fluids, 8(9), 1037–1050, 1988. doi.org/…
11. Rhie C.M., Chow W.L. Numerical Study of the Turbulent Flow Past an Airfoil with Trailing Edge Separation. AIAA Journal, 21(11), 1525–1532, 1983. doi.org/…
12. Majumdar S. Role of Underrelaxation in Momentum Interpolation for Calculation of Flow with Nonstaggered Grids. Numerical Heat Transfer, 13(1), 125–132, 1988. doi.org/…
13. Bartholomew P., Denner F., Abdol-Azis M.H., Marquis A., van Wachem B.G.M. Unified Formulation of the Momentum-Weighted Interpolation for Collocated Variable Arrangements. Journal of Computational Physics, 375, 177–208, 2018. doi.org/…
14. Hanimann L., Mangani L., Casartelli E., Mauri R. A Consistent and Implicit Rhie–Chow Interpolation for Drag Forces in Coupled Multiphase Solvers. International Journal of Turbomachinery, Propulsion and Power, 6(2), 7, 2021. doi.org/…
15. Ferziger J.H., Perić M., Street R.L. Computational Methods for Fluid Dynamics. 4th ed. Springer, 2020.
16. Moukalled F., Mangani L., Darwish M. The Finite Volume Method in Computational Fluid Dynamics. Springer, 2016.
17. NASA Glenn Research Center. Viscosity. Official educational/technical resource on gas viscosity and Sutherland’s law.
18. Sutherland W. The Viscosity of Gases and Molecular Force. Philosophical Magazine, Series 5, 36(223), 507–531, 1893. doi.org/…
19. IAPWS. Release on the IAPWS Formulation 2008 for the Viscosity of Ordinary Water Substance (R12-08, revised 2018).
20. ASTM International. ASTM D341. Standard Practice for Viscosity-Temperature Equations and Charts for Liquid Petroleum or Hydrocarbon Products.
21. NIST. Digital Twins. National Institute of Standards and Technology, thematic program materials.
22. ISO. ISO 23247-1:2021. Automation Systems and Integration — Digital Twin Framework for Manufacturing — Part 1: Overview and General Principles.
23. Glaessgen E., Stargel D. The Digital Twin Paradigm for Future NASA and U.S. Air Force Vehicles. AIAA 2012-1818, 2012.
24. ASME. Standard for Verification and Validation in Computational Fluid Dynamics and Heat Transfer (ASME V&V 20).
25. Roache P.J. Perspective: A Method for Uniform Reporting of Grid Refinement Studies. Journal of Fluids Engineering, 116(3), 405–413, 1994. doi.org/…
26. Celik I.B., Ghia U., Roache P.J., Freitas C.J., Coleman H., Raad P.E. Procedure for Estimation and Reporting of Uncertainty Due to Discretization in CFD Applications. Journal of Fluids Engineering, 130(7), 078001, 2008. doi.org/…
27. Sander R. Compilation of Henry's Law Constants (Version 5.0.0) for Water as Solvent. Atmospheric Chemistry and Physics, 23, 10901–12440, 2023. doi.org/…
28. Sagaut P. Large Eddy Simulation for Incompressible Flows: An Introduction. 3rd ed. Springer, 2006.
29. Shur M.L., Spalart P.R., Strelets M.K., Travin A.K. A Hybrid RANS-LES Approach with Delayed-DES and Wall-Modelled LES Capabilities. International Journal of Heat and Fluid Flow, 29(6), 1638–1649, 2008. doi.org/…