Журнал вычислительной математики и математической физики, 2021, T. 61, № 4, стр. 658-665

Об одном подходе к решению пространственных задач гидродинамики с учетом упругих процессов

А. Ю. Круковский¹, Ю. А. Повещенко¹, В. О. Подрыга^1, *, П. И. Рагимли¹

¹ ИПМ РАН
125047 Москва, Миусская пл., 4, Россия

Поступила в редакцию 21.03.2020
После доработки 08.06.2020
Принята к публикации 18.10.2020

DOI: 10.31857/S0044466921040062

Аннотация

Построена конечно-разностная аппроксимация упругих сил на разнесенных лагранжевых сетках, основанная на методе опорных операторов. Для векторов смещений на нерегулярных сетках, на топологическую и геометрическую структуру которых наложены минимальные разумные ограничения, применительно к разностным схемам для задач теории упругости построены аппроксимации операций векторного анализа. С учетом энергетического баланса среды построенные семейства интегрально-согласованных аппроксимаций операций векторного анализа достаточны для дискретного моделирования этих процессов. Рассматриваются схемы, как использующие тензор напряжений в явном виде, так и разделяющие его на шаровую и сдвиговую компоненты (давление и девиатор). Последнее используется для построения однородных алгоритмов, применимых как для твердого тела, так и для испаренной фазы. При построении аппроксимаций используется линейная теория упругости. В явном виде получены результирующие силы в пространственной геометрии. Приведены расчеты распространения звуковых волн в алюминиевой пространственно-трехмерной ортогональной пластине вследствие торцевого удара. Библ. 13. Фиг. 2.

Ключевые слова: метод опорных операторов, трехмерные конечно-разностные схемы, свойство консервативности, лагранжева сетка разнесенного типа.

1. ВВЕДЕНИЕ

В работе использована идея построения численной методики для сквозного расчета совокупности процессов, возникающих при воздействии на вещество, находящееся первоначально в конденсированном состоянии, с высокоинтенсивным потоком энергии, переносимой излучением (см. [1]). Такого рода методики широко применяются в области физики высоких плотностей энергий. Например, взаимодействие лазерного излучения с материалом твердой мишени может протекать в различных режимах в зависимости от интенсивности воздействия и общей поглощенной энергии. Если интенсивность лазерного импульса достаточно велика, чтобы ионизировать испаренное вещество, то над поверхностью мишени возникает слой плазмы – “корона” или “подушка”, в которой поглощается основная доля лазерной мощности. Перенос энергии в плотное вещество мишени происходит за счет процессов электронной и лучистой теплопроводностей, а также из-за импульса отдачи аблированной плазмы. По веществу распространяется ударная волна, амплитуда которой уменьшается по мере продвижения в глубь мишени. На некотором расстоянии от поверхности вещества температура и давление за фронтом оказываются недостаточными, чтобы осуществить фазовый переход, и вещество остается в твердом состоянии. В этом случае для корректного описания динамики мишени необходим учет упругих сил и сдвиговых напряжений.

В настоящей работе мы рассматриваем реализации дискретной модели упругих сил, основанной на методе опорных операторов (см. [2]–[9]). Тестовые расчеты проведены на трехмерных разностных сетках, состоящих из шестигранников.

2. ПОСТАНОВКА ЗАДАЧИ И ФИЗИЧЕСКИЕ ПРИБЛИЖЕНИЯ

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

(2.1)

$\frac{{\partial \eta }}{{\partial t}} = \eta \operatorname{div} {\mathbf{v}},$

(2.2)

$\rho \frac{{\partial {\mathbf{v}}}}{{\partial t}} = \operatorname{div} \sigma = - \nabla P + \operatorname{div} \zeta ,$

(2.3)

$\rho \frac{{\partial \xi }}{{\partial t}} = \operatorname{tr} (\sigma {{t}_{{v}}}) = - P\operatorname{div} {\mathbf{v}} + \operatorname{tr} (\zeta {{t}_{{v}}}) - \operatorname{div} {\mathbf{W}} + q,$

(2.4)

$\rho \frac{\partial }{{\partial t}}\left( {\xi + \frac{{{{{\mathbf{v}}}^{2}}}}{2}} \right) = - \operatorname{div} (P{\mathbf{v}}) + \operatorname{div} (\zeta {\mathbf{v}}) - \operatorname{div} {\mathbf{W}} + q.$

Здесь t – время, $d{\text{/}}dt = \partial {\text{/}}\partial t + {\mathbf{v}}\nabla $ – субстанциональная производная, tr() – след тензора, ${\mathbf{v}}$ – вектор скорости, ${{t}_{{\mathbf{v}}}} = \frac{1}{2}\left( {\frac{{d{\mathbf{v}}}}{{d{\text{r}}}} + \nabla {\mathbf{v}}} \right)$ – симметризованный тензор скоростей деформаций, $\rho $ – плотность, $\eta = 1{\text{/}}\rho $ – удельный объем, $P$ – давление, $\xi $ – удельная внутренняя энергия, учитывающая деформационные процессы, ${\mathbf{W}}$ – поток тепла, $q$ – объемная плотность мощности сторонних источников.

Приведение тензора напряжений к виду с разделенными шаровой ($P$) и сдвиговой (девиатор – $\zeta $) компонентами:

$\sigma = - P\delta + \zeta ,\quad P = - \frac{1}{3}\operatorname{tr} (\sigma )$

($\delta $ – метрический тензор), является необходимым шагом для построения однородной разностной схемы, применимой как для твердого вещества, где вклады от диагонального ($ - P\delta $) и недиагонального ($\zeta $) членов (отвечающих за объемные и сдвиговые деформации соответственно) сравнимы, так и для плазмы, где сдвиговые напряжения отсутствуют, но для корректного описания зависимости давления от плотности необходимо использование сложных моделей (см. [10]). Таким образом, вклады сдвиговых деформаций ($\operatorname{div} \zeta $) аппроксимируются либо на недеформированной сетке (см. [11]), либо в условиях отсутствия сдвиговой компоненты (модуль сдвига $\mu = 0$), что покрывает всю область параметров для вещества.

Тензор напряжений при твердотельных деформациях имеет вид (см. [11])

$\sigma = 2\mu {{t}_{{\mathbf{U}}}} + \nu \operatorname{tr} ({{t}_{{\mathbf{U}}}})\delta ,\quad {{t}_{{\mathbf{U}}}} = \frac{1}{2}\left( {\frac{{d{\mathbf{U}}}}{{dr}} + \nabla {\mathbf{U}}} \right).$

Здесь $\mu $ (модуль сдвига) и $\nu $ – неотрицательные коэффициенты Ламе, ${\mathbf{U}}$ – вектор смещений относительно недеформированной сетки, ${{t}_{{\mathbf{U}}}}$ – симметризованный тензор смещений. В то же время в твердом теле, а также жидкости или плазме (в которых $\mu = 0$) всегда справедливо

$\sigma = - P\delta + \zeta ,\quad \zeta = 2\mu \left( {{{t}_{{\mathbf{U}}}} - \frac{1}{3}\operatorname{tr} ({{t}_{{\mathbf{U}}}})\delta } \right).$

При адиабатических деформациях внутренняя энергия единицы массы тела представима в виде (см. [11])

$\xi = \frac{{{{k}_{{ad}}}}}{{2\rho }}{{(\operatorname{tr} ({{t}_{{\mathbf{U}}}}))}^{2}} + \frac{1}{{4\mu \rho }}\operatorname{tr} {{(\zeta )}^{2}} = \varepsilon (\rho ,T) + \frac{1}{{4\mu \rho }}\operatorname{tr} {{(\zeta )}^{2}},$

где ${{k}_{{ad}}}$ – адиабатический модуль сжатия. Зависимость внутренней энергии от плотности и температуры $\varepsilon (\rho ,T)$ считается заданной.

3. МЕТОД ОПОРНЫХ ОПЕРАТОРОВ ДЛЯ АППРОКСИМАЦИИ УПРУГИХ ПРОЦЕССОВ

3.1. Метрические сетки метода опорных операторов

Для типов сеток, состоящих из ячеек ${(\Omega )}$, образованных узлами ${(\omega )}$, гранями ${(\sigma )}$ и ребрами ${(\lambda )}$, характерно наличие замкнутой сопряженной (“сдвинутой”) сетки, состоящей, например, из доменов ${d(\omega )}$ вокруг узлов $\omega $ ( фиг. 1).

Фиг. 1.

Построение базисов.

Грани узлового домена определяются метрическим оператором сетки $\sigma (\lambda ) = \sum\nolimits_{\varphi (\lambda )} {{{V}_{\varphi }}{\mathbf{e}}_{\varphi }^{'}} (\lambda )$ (см. также ниже). Базисы $\varphi (\lambda )$ здесь попарно входят в ячейки ${\Omega (\lambda )}$, примыкающие к ребру λ. Метрическая калибровка разностной сетки состоит в выборе объемов базисов с естественным условием нормировки $\sum\nolimits_{\varphi (\Omega )} {{{V}_{\varphi }}} = {{V}_{\Omega }}$. Она определяет конструкцию замкнутой сопряженной сетки для различных классов сеток. Это треугольно-четырехугольные 2D сетки, тетраэдральные, параллелепипедные, призматические (и т.д.), 3D сетки, а также их мортарные сшивки, адаптация (с введением новых узлов в ячейках $\Omega $, называемых дуальными) с сохранением самосопряженности и знакоопределенности соответствующих “дивергентно-градиентных” операций векторного анализа континуальных краевых задач. Дальнейшее изложение носит общий характер, конкретный выбор локальных базисных объемов ${{{V}_{\varphi }}}$ иллюстрируется на примере треугольно-четырехугольной 2D сетки.

В области О введем семейство нерегулярных разностных сеток. Ограничимся случаем, когда сетка состоит из треугольных и четырехугольных ячеек ${(\Omega )}$, базисов ${(\varphi )}$, узлов ${(\omega )}$, ребер ${(\lambda )}$ и связанных с ними ${(\sigma (\lambda ))}$ границами балансовых узловых доменов $d(\omega )$ (фиг. 1).

Базисы $\varphi $ создаются системой исходных (ковариантных) ортов ${\mathbf{e}}(\lambda )$, образованных ребрами. Под центрами ячеек $\Omega $ и ребер $\lambda $ будем понимать средние арифметические радиус-векторов узлов $\omega $, их образующих. Кривая, соединяющая эти центры (двух смежных через ребро ячеек или ячейку с граничным ребром $\partial \lambda $), представляет собой поверхность

$\sigma (\lambda ) = \sum\limits_{\varphi (\lambda )} {{{v}_{\varphi }}{\mathbf{e}}_{\varphi }^{'}} (\lambda ){\text{ ,}}$

ориентированную так же, как и орт ${\mathbf{e}}(\lambda )$. Здесь ${\mathbf{e}}_{\varphi }^{'}(\lambda )$ – орты взаимных (контравариантных) базисов по отношению к исходным базисам, образованным ортами ${\mathbf{e}}(\lambda )$. Базисный объем дается формулой ${{v}_{\varphi }} = \frac{1}{6}\left| {{\mathbf{e}}({{\lambda }_{1}}) \times {\mathbf{e}}({{\lambda }_{2}})} \right|$ для треугольной ячейки $\Omega $, содержащей базис $\varphi $, и ${{v}_{\varphi }} = \frac{1}{4}\left| {{\mathbf{e}}({{\lambda }_{1}}) \times {\mathbf{e}}({{\lambda }_{2}})} \right|$ для четырехугольной ячейки, если ${{\lambda }_{1}}(\varphi )$ и ${{\lambda }_{2}}(\varphi )$ – ребра, образующие базис $\varphi $. Наконец, $\sum\nolimits_{\varphi (\lambda )} {\kern 1pt} $ – суммирование по всем базисам $\varphi $, в образовании которых приняло участие ребро $\lambda $. Замкнутые вокруг узла $\omega $ поверхности ${\sigma (\lambda (\omega ))}$ образуют узловые домены ${d(\omega )}$.

К узлам отнесем искомую сеточную функцию u. На ребрах выделим положительное направление (фиг. 1) и отнесем к ним сеточную функцию

${{\Delta }_{\lambda }}{\mathbf{u}} = - \sum\limits_{\omega (\lambda )} {{{s}_{\lambda }}} (\omega ){{{\mathbf{u}}}_{\omega }} = {{{\mathbf{u}}}_{{\omega {\kern 1pt} '}}} - {{{\mathbf{u}}}_{\omega }}.$

Сеточные тензорные поля X задаются своими представлениями в базисах ${{{X}_{\varphi }}}$.

Внутреннюю дивергенцию тензорного поля DIN: $(\varphi ) \to (\omega )$ определим, аппроксимируя теорему Гаусса на ${d(\omega )}$:

$\operatorname{DIN} X = \sum\limits_{\lambda (\omega )} {{{s}_{\lambda }}} (\omega ){{\tau }_{X}}(\lambda ),\quad {{\tau }_{X}}(\lambda ) = \sum\limits_{\varphi (\lambda )} {{{v}_{\varphi }}} ({{\mathbf{e}}_{\varphi }^{'}}(\lambda ),{{{X}_{\varphi }}}),$

где $\sum\nolimits_{\lambda (\omega )} {\kern 1pt} $ – суммирование по всем ребрам $\lambda $, имеющим общий узел $\omega $. Обозначая через ${{(\;)}_{\Delta }}$ аппроксимацию соответствующих дифференциальных выражений, имеем

$\begin{gathered} {{{{\left( {\int\limits_O {{\text{tr}}} \left( {\nabla {\mathbf{U}}{{X}^{T}}} \right)d{v}} \right)}}_{\Delta }} = - {{{\left( {\int\limits_O {\mathbf{U}} \operatorname{div} Xdv{\text{ }} - {\text{ }}\int\limits_{\partial O} {(X{\mathbf{U}},{\mathbf{ds}})} } \right)}}_{\Delta }}} = - \sum\limits_\omega {({{{\mathbf{u}}}_{\omega }},\operatorname{DIN} X)} = \\ = \sum\limits_\varphi {{{v}_{\varphi }}} \operatorname{tr} \left( {{{{\left( {\frac{{D{\mathbf{u}}}}{{D{\mathbf{r}}}}} \right)}}_{\varphi }}{{X}_{\varphi }}} \right) = \sum\limits_\varphi {{{v}_{\varphi }}} \operatorname{tr} (\nabla {{{\mathbf{u}}}_{\varphi }}X_{{_{\varphi }}}^{T}) \equiv \left\langle {\nabla {\mathbf{u}},{\text{ }}X} \right\rangle = \left\langle {X,\nabla {\mathbf{u}}} \right\rangle = \sum\limits_\varphi {{{v}_{\varphi }}} \operatorname{tr} ({{t}_{{\Delta u\varphi }}}{{X}_{\varphi }})\left| {_{{{{X}_{\varphi }} = X_{\varphi }^{T}}}} \right. = \left\langle {{{t}_{{\Delta u}}},X} \right\rangle . \\ \end{gathered} $

Здесь $\left\langle {\nabla {\mathbf{u}},{\text{X}}} \right\rangle $ определяется как скалярное произведение сеточных тензорных полей, аппроксимирующее ${{{{\left( {\int_O {\operatorname{tr} } (\nabla {\mathbf{U}}{{X}^{T}})dv} \right)}}_{\Delta }}}$.

Тензорные поля $\frac{{D{\mathbf{u}}}}{{D{\mathbf{r}}}},$ $\nabla {\mathbf{u}}$ и ${{t}_{{\Delta {\mathbf{u}}}}}$, а также тензорное поле напряжений ${{X}_{{\Delta {\mathbf{u}}}}}$ даются своими представлениями в базисах

$\begin{gathered} {{\left( {\frac{{D{\mathbf{u}}}}{{D{\mathbf{r}}}}} \right)}_{\varphi }} = {\text{ }}\sum\limits_{\lambda (\varphi )} {{{\Delta }_{\lambda }}} {\mathbf{u}}\;\centerdot \;{\mathbf{e}}_{\varphi }^{'}(\lambda ),\quad \nabla {{{\mathbf{u}}}_{\varphi }} = \sum\limits_{\lambda (\varphi )} {{\mathbf{e}}_{\varphi }^{'}(\lambda )\;\centerdot \;{{\Delta }_{\lambda }}{\mathbf{u}}} , \\ {{t}_{{\Delta {\mathbf{u}}\varphi }}} = \frac{1}{2}\left( {{{{\left( {\frac{{D{\mathbf{u}}}}{{D{\mathbf{r}}}}} \right)}}_{\varphi }} + \nabla {{{\mathbf{u}}}_{\varphi }}} \right),\quad {{X}_{{\Delta {\mathbf{u}}\varphi }}} = 2\mu {{t}_{{\Delta {\mathbf{u}}\varphi }}} + \upsilon \operatorname{tr} ({{t}_{{\Delta {\mathbf{u}}\varphi }}})\delta . \\ \end{gathered} $

Под $\sum\nolimits_{\lambda (\varphi )} {\kern 1pt} $ понимается суммирование по ребрам $\lambda $, образующим базис $\varphi $.

Произвольному сеточному векторному полю инкрементов ${{\vec {\partial }}_{\lambda }} \in {{H}_{\lambda }}$ на ребрах $\lambda $ (определяющему поле смещений ${{\vec {\partial }}^{\& }}$) сопоставим поле симметризованного тензора смещений ${{t}_{\partial }}$ и поле тензора напряжений ${{X}_{\partial }}$ по формулам

${{t}_{{\partial \phi }}} = \frac{1}{2}\mathop \sum \limits_{\lambda (\varphi )} \left( {{{{\vec {\partial }}}_{\lambda }}\;\centerdot \;~{{e}_{\phi }}(\lambda ) + {\mathbf{e}}_{\phi }^{'}(\lambda )\;\centerdot \;~{{{\vec {\partial }}}_{\lambda }}} \right) = \frac{1}{2}\left( {{{{\left( {\frac{{D{{{\vec {\partial }}}^{\& }}}}{{Dr}}} \right)}}_{\phi }} + \nabla \vec {\partial }_{\phi }^{\& }} \right),$

${{X}_{{\partial \phi }}} = 2\mu ~{{t}_{{\partial \phi }}} + \upsilon \operatorname{tr} ({{t}_{{\partial \phi }}})\delta .$

Для вариации континуального векторного поля смещения $\delta {\mathbf{U}}$ на расстоянии $\delta {\mathbf{r}}$ справедливо $\delta {\mathbf{U}} = \left( {\frac{{d{\mathbf{U}}}}{{d{\mathbf{r}}}},\delta {\mathbf{r}}} \right)$. В силу определения взаимного базиса $({\mathbf{e}}(\lambda ),{\mathbf{e}}_{\varphi }^{'}(\lambda {\kern 1pt} ')) = {{\delta }_{{\lambda \lambda {\kern 1pt} '}}}$ для ребра ${\lambda (\varphi )}$ очевиден сеточный аналог этого тождества:

${{\Delta }_{\lambda }}{\mathbf{u}} = \left( {{{{\left( {\frac{{D{\mathbf{u}}}}{{D{\mathbf{r}}}}} \right)}}_{\varphi }},{\mathbf{e}}(\lambda )} \right).$

Уточним также силу ${{\tau }_{X}}(\lambda )$, действующую на поверхность $\sigma (\lambda )$ в поле напряжений ${{X}_{{\Delta {\mathbf{u}}}}}$:

$\begin{gathered} {{{\mathbf{\tau }}}_{{X\Delta {\mathbf{u}}}}}(\lambda ) = \mathop \sum \limits_{\phi \left( \lambda \right)} \,{{{v}}_{\phi }}({\mathbf{e}}_{\phi }^{'}(\lambda ),~{{X}_{{\Delta {\mathbf{u}}\phi }}}) = \mathop \sum \limits_{\phi \left( \lambda \right)} \,{{{v}}_{\phi }}\{ ~\mu ~\mathop \sum \limits_{\lambda '\left( \phi \right)} \,[({\mathbf{e}}_{\phi }^{'}(\lambda ),~{\mathbf{e}}_{\phi }^{'}(\lambda {\kern 1pt} '))~{{\Delta }_{{\lambda {\kern 1pt} '}}}{\mathbf{u}} + ~ \\ + \;({\mathbf{e}}_{\phi }^{'}(\lambda ),~{{\Delta }_{{\lambda {\kern 1pt} '}}}{\mathbf{u}})~{\mathbf{e}}_{\phi }^{'}(\lambda {\kern 1pt} ')] + \upsilon {\mathbf{e}}_{\phi }^{'}(\lambda )\mathop \sum \limits_{\lambda '(\phi )} \,({\mathbf{e}}_{\phi }^{'}(\lambda {\kern 1pt} '),~{{\Delta }_{{\lambda {\kern 1pt} '}}}u)\} . \\ \end{gathered} $

На поле инкрементов ${{\vec {\partial }}_{\lambda }} \in {{H}_{{\lambda O}}}$ (не вращающем среду твердотельно) эта сила определяет самосопряженный, положительно определенный метрический оператор ${{G}_{O}}$: ${\text{(}}\lambda {\text{)}} \to {\text{(}}\lambda {\text{)}}$, ${{G}_{O}} = G_{O}^{*} > 0$:

${{{\mathbf{\tau }}}_{{X\partial }}}(\lambda ) = \mathop \sum \limits_{\varphi (\lambda )} \,{{{v}}_{\varphi }}({\mathbf{e}}_{\varphi }^{'}(\lambda ),{{X}_{{\partial \varphi }}}),\quad {{{\mathbf{\tau }}}_{{X\partial }}} = {{G}_{O}}\vec {\partial }.$

Для ${{\overrightarrow {\partial 1} }_{\lambda }}$ и ${{\overrightarrow {\partial 2} }_{\lambda }}$ из ${{H}_{{\lambda O}}}$ скалярные произведения ${{(\overrightarrow {\partial 1} {\text{\;,\;}}\overrightarrow {\partial 2} )}_{\lambda }}$ и $\left\langle {\left\langle {\overrightarrow {\partial 1} {\text{\;,\;}}\overrightarrow {\partial 2} } \right\rangle } \right\rangle $ определяются по формулам

${{(\overrightarrow {\partial 1} ,\overrightarrow {\partial 2} )}_{\lambda }} = \mathop \sum \limits_\lambda \,({{\overrightarrow {\partial 1} }_{\lambda }},{{\overrightarrow {\partial 2} }_{\lambda }}),$

$\begin{gathered} {{(\overrightarrow {\partial 1} ,{{G}_{O}}\overrightarrow {\partial 2} )}_{\lambda }} = \mathop \sum \limits_\phi \,{{{v}}_{\phi }}~[2\mu \operatorname{tr} ({{t}_{{\partial 1\phi }}}~{{t}_{{\partial 2\phi }}}) + \upsilon \operatorname{tr} ({{t}_{{\partial 1\phi }}})\operatorname{tr} ({{t}_{{\partial 2\phi }}})] = \\ = {{({{G}_{O}}\overrightarrow {\partial 1} ,\overrightarrow {\partial 2} )}_{\lambda }} = \left\langle {\nabla {{{\overrightarrow {\partial 1} }}^{\& }},{{X}_{{\partial 2}}}} \right\rangle = \left\langle {\nabla {{{\overrightarrow {\partial 2} }}^{\& }},{{X}_{{\partial 1}}}} \right\rangle \equiv \left\langle {\left\langle {\overrightarrow {\partial 1} ,\overrightarrow {\partial 2} } \right\rangle } \right\rangle . \\ \end{gathered} $

Энергия деформации среды, производимой полем инкрементов ${{\vec {\partial }}_{\lambda }}$, есть энергия метрического оператора в этом поле:

${{({{G}_{O}}\vec {\partial },\vec {\partial })}_{\lambda }} = \left\langle {\left\langle {\vec {\partial },\vec {\partial }} \right\rangle } \right\rangle .$

3.2. Поворотно-нейтральные разностные схемы

Определим твердотельное вращение (2D и 3D сетка) как возмущение поля инкрементов ${{\overrightarrow \partial }_{\lambda }} \in {{H}_{\lambda }}$ следующим образом:

${{\vec {\partial }}_{\lambda }} \to {{\vec {\partial }}_{\lambda }} + [{\mathbf{\omega }}~~ \times ~{\mathbf{~e}}(\lambda )],\quad {\mathbf{\omega }} = {\text{const}}{\text{.}}$

Симметризованный тензор смещений ${{t}_{{\partial \varphi }}} = \frac{1}{2}\sum\nolimits_{\lambda (\varphi )}^{} {({{{\vec {\partial }}}_{\lambda }} \bullet {\mathbf{e}}_{\varphi }^{'}(\lambda ) + {\mathbf{e}}_{\varphi }^{'}(\lambda ) \bullet {{{\vec {\partial }}}_{\lambda }})} $ испытает возмущение ${{t}_{{\partial \varphi }}} \to {{t}_{{\partial \varphi }}} + {{t}_{{\omega \varphi }}}$.

Здесь

${{t}_{{\omega \varphi }}} = \frac{1}{2}\sum\limits_{\lambda (\varphi )} {([\omega ~~ \times ~~{\mathbf{e}}(\lambda )] \bullet {\mathbf{e}}_{\varphi }^{'}(\lambda ) + {\mathbf{e}}_{\varphi }^{'}(\lambda ) \bullet [\omega ~~ \times ~~{\mathbf{e}}(\lambda )])} .$

Для произвольного вектора a имеем

$\begin{gathered} {{t}_{{\omega \varphi }}}{\mathbf{a}} = \frac{1}{2}\left\{ {\left[ {{\mathbf{\omega }} \times \sum\limits_{\lambda (\varphi )} {~{\mathbf{e}}(\lambda )({\mathbf{e}}_{\varphi }^{'}(\lambda ),{\mathbf{a}})} } \right] + \sum\limits_{\lambda (\varphi )} ~ {\mathbf{e}}_{\varphi }^{'}(\lambda )([{\mathbf{\omega }}~ \times ~~{\mathbf{e}}(\lambda )],{\mathbf{a}})} \right\} = \\ = \frac{1}{2}\left\{ {[{\mathbf{\omega }}~~ \times ~~{\mathbf{a}}] + \sum\limits_{\lambda (\varphi )} ~ {\mathbf{e}}_{\varphi }^{'}(\lambda )([{\mathbf{a}}~~ \times ~~{\mathbf{\omega }}],{\mathbf{e}}(\lambda ))} \right\} = \frac{1}{2}\left\{ {[{\mathbf{\omega }}~~ \times ~~{\mathbf{a}}] + [{\mathbf{a}}~~ \times ~~{\mathbf{\omega }}]} \right\} = 0, \\ \end{gathered} $

т.е. тензор ${{t}_{{\omega \varphi }}} = 0$. Следовательно, тензор ${{t}_{{\partial \varphi }}}$ инвариантен к твердотельным вращениям $[\omega ~~ \times ~~{\mathbf{e}}(\lambda )],$ $\omega = {\text{const}}$. Такое свойство разностной схемы, смещения в которой входят только через инвариантный к твердотельным вращениям симметризованный тензор ${{t}_{{\partial \varphi }}}$ в базисах сетки, называется поворотной нейтральностью.

Упругие силовые и энергетические характеристики таких схем инвариантны не только к параллельному переносу, но и к твердотельным вращениям. Отсюда же следует, что для строгой положительной определенности оператора $\operatorname{DIN} {{X}_{{\mathbf{u}}}}({{X}_{{\Delta {\mathbf{u}}\varphi }}} = 2\mu ~{{t}_{{\Delta {\mathbf{u}}\varphi }}} + \upsilon \operatorname{tr} {\text{(}}~{{t}_{{\Delta {\mathbf{u}}\varphi }}})~\delta )$ необходимо помимо параллельного переноса сетки $({{\Delta }_{\lambda }}{\mathbf{u}} \equiv 0)$ исключить ее твердотельное вращение (${{\Delta }_{\lambda }}{\mathbf{u}} = [\omega \times {\mathbf{e}}(\lambda )],$ $\omega = {\text{const}}$). В трехмерном случае это достигается принципом “трех гвоздей”, т.е. в трех узлах, образующих плоскость, должны быть заданы нулевые смещения. Соответственно в двумерном случае должна быть закреплена прямая (двумя узлами).

3.3. Скалярно-дивергентные задачи

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

В пространственной области O с границей $\partial O$ рассмотрим скалярно-дивергентную задачу

$\operatorname{div} {{{\mathbf{X}}}_{u}} = f(r),\quad {{{\mathbf{X}}}_{u}} = K\nabla u$

(с некоторыми граничными условиями) наряду с соответствующим интегральным соотношением:

${\int\limits_{\text{O}} {({\mathbf{X}}\nabla u){\kern 1pt} dv} + \int\limits_{\text{O}} {u\operatorname{div} {\mathbf{X}}dv} = \int\limits_{\partial {\text{O}}} {u({\mathbf{X}}{\text{,}}{\mathbf{ds}})} }.$

Здесь u – скаляр (температура, давление и т.п.), X – произвольный вектор, ${{{{\mathbf{X}}}_{u}}}$ – потоки, вызванные градиентом функции u в среде со свойствами проводимости, определяемыми тензором K.

Как и ранее (фиг. 1) в области О вводится семейство нерегулярных разностных сеток, обладающих метрическими свойствами, и соответствующих сеточных функций. Внутреннюю дивергенцию векторного поля DIN: $(\varphi ) \to (\omega )$ определим аналогично п. 3.1, аппроксимируя теорему Гаусса на ${d(\omega )}$. Сеточное векторное поле X задается своими представлениями в базисах ${{{{\mathbf{X}}}_{\varphi }}}$. Обозначая через ${{(\;)}_{\Delta }}$ аппроксимацию соответствующих дифференциальных выражений, имеем

${{{{\left( {\int\limits_O {({\mathbf{X}},\nabla u)d{v}} } \right)}}_{\Delta }} = - {{{\left( {\int\limits_O {u\operatorname{div} {\mathbf{X}}d{v}} {\text{ }} - {\text{ }}\int\limits_{\partial O} {u({\mathbf{X}},{\mathbf{d}}s)} } \right)}}_{\Delta }}} = - \sum\limits_\omega {({{u}_{\omega }},\operatorname{DIN} {\mathbf{X}})} = \sum\limits_\varphi {{{{v}}_{\varphi }}} ({{{\mathbf{X}}}_{\varphi }},\operatorname{GRAD} u).$

Градиентное векторное поле GRAD: ${(\omega ) \to (\varphi )}$ дается своими представлениями в базисах:

$\operatorname{GRAD} u = \sum\limits_{\lambda (\varphi )} {{{\Delta }_{\lambda }}u{\text{ }}{\mathbf{e}}_{\varphi }^{'}(\lambda )} ,\quad {{\Delta }_{\lambda }}u = - \sum\limits_{\omega (\lambda )} {{{s}_{\lambda }}(\omega ){{u}_{\omega }}} = {{u}_{{\omega *}}} - {{u}_{\omega }}.$

Полагая в базисах $\varphi $ в качестве ${{{{\mathbf{X}}}_{\varphi }}}$ векторное поле ${{{{\mathbf{X}}}_{{\nu \varphi }}} = {{K}_{\varphi }}\operatorname{GRAD} \nu }$, получаем самосопряженный неотрицательный оператор ${\operatorname{DIN} {{{\mathbf{X}}}_{\nu }}:{\text{(}}\omega {\text{)}} \to {\text{(}}\omega {\text{)}}}$ или ${\operatorname{DIN} K\operatorname{GRAD} :{\text{(}}\omega {\text{)}} \to {\text{(}}\omega {\text{)}}}$. Здесь потоковое векторное поле ${{{{\mathbf{X}}}_{\nu }}}$ дается своими компонентами в базисах ${{{{\mathbf{X}}}_{{\nu \varphi }}}}$. Оно определяется градиентными свойствами скалярной сеточной функции $\nu $, заданной в узлах $\omega $, и сеточным тензорным полем проводимости $K$, задаваемым своими представлениями в базисах ${{K}_{\varphi }}$. Этот оператор будет строго положительным, если хотя бы в одном граничном узле связной разностной сетки задана первая краевая задача, т.е. в этом граничном узле скалярная сеточная функция обращается в ноль.

3.4. Дифференциально-разностная схема

Введем обозначения для функций сетки метода опорных операторов (п. 3, см. также фиг. 1). К ее узлам $\omega $ будем относить ранее представленные в континуальной модели компоненты вектора скорости ${\mathbf{v}}$ и смешения ${\mathbf{U}}$, а также термодинамические величины $\rho ,\;\eta ,\;P,\;T,\;\xi ,\;\varepsilon $. Аналогичные разностные схемы, где термодинамические функции заданы в ячейках, представлены в [1], [12], [13]. Дифференциально-разностная схема в узлах сетки для уравнений (2.1)–(2.3) с “нулевыми” граничными условиями может быть записана как

$\frac{{dV}}{{dt}} = \operatorname{DIN} {\mathbf{v}},$

$\frac{d}{{dt}}(m{\mathbf{v}}) = \sum\limits_{\varphi (\omega )} {{{V}_{\varphi }}\operatorname{GRAD} P} + \operatorname{DIN} \zeta ,$

$\frac{d}{{dt}}(m\zeta ) = - P\operatorname{DIN} {\mathbf{v}} + \sum\limits_{\varphi (\omega )} {{{V}_{\varphi }}\operatorname{tr} ({{\zeta }_{\varphi }}{{t}_{{\Delta {v}\varphi }}})} - \operatorname{DIN} {\mathbf{W}} + Vq.$

Здесь приузловой объем ${{V}_{\omega }} = \sum\nolimits_{\varphi (\omega )}^{} {{{V}_{\varphi }}} $ определяется базисными объемами ${{V}_{\varphi }}$. Приузловая масса постоянна ${{m}_{\omega }} = {{\rho }_{\omega }}{{V}_{\omega }} = {\text{const}}$. Скорость ${{{\mathbf{v}}}_{\varphi }}$ в базисе $\varphi $ с центральным узлом $\omega $ берется как ${{{\mathbf{v}}}_{\varphi }} = {{{\mathbf{v}}}_{{\omega \,(\varphi )}}}$. Сеточные тензорные поля сдвиговых напряжений задаются представлениями в базисах ${{\zeta }_{\phi }} = 2{{\mu }_{\phi }}\left( {{{t}_{{\Delta {\mathbf{U}}\phi }}} - \frac{1}{3}\operatorname{tr} ({{t}_{{\Delta {\mathbf{U}}\phi }}})\delta } \right)$. Тепловой поток в базисах определяется как ${{{\mathbf{W}}}_{\varphi }} = - {{\kappa }_{\varphi }}\operatorname{GRAD} T$, где ${{\kappa }_{\varphi }}({{P}_{\omega }},{{T}_{\omega }})$ – теплопроводность в базисе $\varphi $ (ячейки $\Omega \supset \varphi $, содержащей этот базис) с центральным узлом $\omega (\varphi )$. Учитывающая адиабатические деформации внутренняя энергия узла $\omega $ с массой m определится как

$m\xi = m\varepsilon (\rho ,T) + \frac{1}{4}\mathop \sum \limits_{\phi (\omega )} \frac{{{{V}_{\phi }}}}{{{{\mu }_{\phi }}}}\operatorname{tr} (\zeta _{\phi }^{2}).$

Здесь модуль сдвига ${{\mu }_{\varphi }}$ на сетке определяется в базисе $\varphi $ аналогично теплопроводности ${{\kappa }_{\varphi }}$.

4. РЕЗУЛЬТАТЫ РАСЧЕТОВ

Рассматривается движение вещества, описываемое уравнениями неразрывности, Эйлера и энергетическими балансами с адиабатическими деформациями твердой фазы в рамках линейной теории упругости (см. [11]).

Для исследования однородной конечно-разностной аппроксимации упругих сил в рамках решения сопряженных задач гидродинамики и упругости (см. [1]) в случае, когда вещество находится в твердом состоянии, проведено моделирование распространения звуковых волн в трехмерной ортогональной пластине вследствие удара, колебания пластины вследствие сдвиговой деформации или деформации изгиба.

Во всех этих случаях для вещества использовались параметры алюминия с плотностью ${{\rho }_{s}} = 2700$ кг/м³, модулем сдвига $\mu = 26.1~$ ГПа и модулем сжатия $k = 72.9~$ ГПа. Соответствующие значения скорости продольного и поперечного звука составляют ${{c}_{p}} = \sqrt {\left( {k + \frac{4}{3}\mu } \right){{\rho }_{s}}} = 6300$ м/с и ${{c}_{s}} = \sqrt {\mu {\text{/}}{{\rho }_{s}}} = 3100~$ м/с. Задача о распространении волны в плоской трехмерной пластине рассматривалась в следующей постановке: $x \in [0;10],$ $y \in [0;5],$ $z \in [0;0.1]$, ${{{v}}_{x}}(t = 0) = 0,$ ${{{v}}_{y}}(t = 0) = {{{v}}_{0}}{{e}^{{ - 4y - {{{(x - 5)}}^{2}}}}},$ ${{{v}}_{z}}(t = 0) = 0$. Здесь координаты выражены в сантиметрах; амплитуда волны ${{{v}}_{0}} = 20$ м/с много меньше скорости звука в алюминии. В начальный момент времени вещество предполагалось недеформированным. На правой границе (фиг. 2) ставилось условие отражения, в то время как на остальных границах задавалось отсутствие внешних сил. Динамика мишени, полученная в расчетах в xyz-геометрии, представлена на фиг. 2. Для анализа вызванных лазерным воздействием гидродинамических и упруго-волновых процессов важна численная информация об обмене импульсом и энергией между испаренной и неиспаренной частями мишени. Для проверки качества методики в этом отношении, аналогично [1], проводился анализ динамики интегральных уравнений внутренней, кинетической и полной энергий, полученных в расчетах на сетке 500 × 200 × 10 (шаг 0.01 см).

Фиг. 2.

Динамика плотности на моменты (сверху вниз, слева направо) 2, 4, 6, 8 и 10 мкс, z = 0.05 см.

5. ЗАКЛЮЧЕНИЕ

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

Список литературы

Цыгвинцев И.П., Круковский А.Ю., Повещенко Ю.А., Гасилов В.А., Бойков Д.С., Попов С.Б. Однородные разностные схемы для сопряженных задач гидродинамики и упругости // Учен. зап. Казан. ун-та. Сер. Физ.-матем. науки. 2019. Т. 161. № 3. С. 377–392.
Самарский A.A., Тишкин В.Ф., Фаворский А.П., Шашков М.Ю. Операторные разностные схемы // Дифференц. ур-ния. 1981. Т. 17. № 7. С. 1317–1327.
Самарский A.A., Тишкин В.Ф., Фаворский А.П., Шашков М.Ю. О представлении разностных схем математической физики в операторной форме // Докл. АН СССР. 1981. Т. 258. № 5. С. 1092–1096.
Коршия Т.К., Тишкин В.Ф., Самарский А.А. и др. Вариационно-операторные разностные схемы для уравнений математической физики. Тбилиси: Изд-во Тбил. ун-та, 1983. 143 с.
Денисов А.А., Колдоба А.В., Повещенко Ю.А. О сходимости разностных схем метода опорных операторов для уравнения Пуассона на обобщенных решениях // Ж. вычисл. матем. и матем. физ. 1989. Т. 29. № 3. С. 371–381.
Денисов А.А., Колдоба А.В., Повещенко Ю.А. О сходимости разностных схем метода опорных операторов для осесимметричного уравнения Пуассона на обобщенных решениях // Ж. вычисл. матем. и матем. физ. 1990. Т. 30. № 10. С. 1477–1486.
Самарский А.А., Колдоба А.В., Повещенко Ю.А., Тишкин В.Ф., Фаворский А.П. Разностные схемы на нерегулярных сетках. Минск: ЗАО “Критерий”, 1996. 275 с.
Колдоба А.В., Повещенко Ю.А., Гасилова И.В., Дорофеева Е.Ю. Разностные схемы метода опорных операторов для уравнений теории упругости // Матем. моделирование. 2012. Т. 24. № 12. С. 86–96.
Повещенко Ю.А., Круковский А.Ю., Подрыга В.О., Головченко Е.Н. Разностные схемы метода опорных операторов для уравнений упругости с азимутальным вращением // Препринты ИПМ им. М.В. Келдыша. 2019. № 10. 36 с.
Новиков В.Г., Никифоров В.Г., Уваров В.Б. Квантово-статистические модели высокотемпературной плазмы и методы расчета росселандовых пробегов и уравнений состояния. М.: Физматлит, 2000. 400 с.
Ландау Л.Д., Лифшиц Е.М. Теоретическая физика. Т. 7. Теория упругости. М.: Наука, 1987. 248 с.
Круковский А.Ю., Повещенко Ю.А., Клочкова Л.В., Сузан Д.В. Оценки сходимости итерационных алгоритмов численного решения трехмерных нестационарных задач магнитной гидродинамики // Препринты ИПМ им. М.В. Келдыша. 2019. № 94. 17 с.
Круковский А.Ю., Новиков В.Г., Цыгвинцев И.П. Численные алгоритмы для решения трехмерных нестационарных задач магнитной гидродинамики // Препринты ИПМ им. М.В. Келдыша. 2014. № 6. 20 с.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Журнал вычислительной математики и математической физики

Архивы выпусков Информация о журнале Отправить рукопись в журнал