Журнал вычислительной математики и математической физики, 2020, T. 60, № 4, стр. 738-751

О термодинамической согласованности и математической корректности в теории упругопластических, сыпучих и пористых сред

В. М. Садовский^*

ИВМ СО РАН
660036 Красноярск, Академгородок, 50/44, Россия

Поступила в редакцию 14.11.2019
После доработки 14.11.2019
Принята к публикации 16.12.2019

DOI: 10.31857/S0044466920040158

Аннотация

Математические модели динамики упругопластических, сыпучих и пористых сред приводятся к вариационным неравенствам для гиперболических дифференциальных операторов, термодинамически согласованных по Годунову. На этой основе формулируется понятие обобщенных решений с диссипативными ударными волнами, строятся априорные оценки гладких решений в характеристических коноидах операторов, дающие представление о корректности постановки задачи Коши и краевых задач с диссипативными граничными условиями, конструируются эффективные вычислительные методы сквозного счета, адаптированные к разрывам решений. Библ. 23. Фиг. 2.

Ключевые слова: динамика, ударная волна, упругость, пластичность, сыпучая среда, пористая среда, термодинамическая согласованность, вариационное неравенство, метод сквозного счета.

1. ВВЕДЕНИЕ

Специальные формулировки нестационарных уравнений механики сплошных сред в виде термодинамически согласованных систем законов сохранения изучались ранее в работах С.К. Годунова, его коллег и учеников [1]–[4]. Они были получены для уравнений газовой динамики, теории упругости, магнитной гидродинамики, электродинамики и др., и оказались весьма полезными при исследовании вопросов корректности постановки краевых задач с начальными данными и граничными условиями, при анализе разрывных решений с ударными волнами, при построении конечно-разностных и конечно-объемных схем для численного решения краевых задач.

Формулировка определяющих уравнений для $m$-мерного вектора $U(t,x)$ функций состояния сплошной среды в виде термодинамически согласованной системы законов сохранения предполагает наличие так называемых производящих потенциалов $\Phi (U)$ и ${{\Psi }_{i}}(U)$, ($i = 1,...,n$; $n$ – пространственная размерность модели). Уравнения записываются в следующей форме:

(1.1)

$\mathcal{D}\left\langle U \right\rangle = g(U),$

где $g$ – заданный $m$-мерный вектор правой части, а $\mathcal{D}$ – квазилинейный дифференциальный оператор первого порядка (угловые скобки указывают на функциональную зависимость):

$\mathcal{D}\left\langle U \right\rangle = \frac{{\partial \varphi (U)}}{{\partial t}} - \sum\limits_{i = 1}^n \,\frac{{\partial {{\psi }_{i}}(U)}}{{\partial {{x}_{i}}}},\quad \varphi (U) = \frac{{\partial \Phi }}{{\partial U}},\quad {{\psi }_{i}}(U) = \frac{{\partial {{\Psi }_{i}}}}{{\partial U}}.$

Вектор $g(U)$, вообще говоря, зависит от переменных $t,\;x$. Производящие потенциалы могут также зависеть от $t$ и $x$ как от параметров, если среда неоднородна или ее свойства меняются со временем. Однако далее для упрощения формул эти аргументы не указываются.

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

(1.2)

$\frac{\partial }{{\partial t}}(U \cdot \varphi (U) - \Phi (U)) = \sum\limits_{i = 1}^n \,\frac{\partial }{{\partial {{x}_{i}}}}(U \cdot {{\psi }_{i}}(U) - {{\Psi }_{i}}(U)) + U \cdot g(U),$

и соответствующим ему интегральным законом сохранения, который согласуется с законами сохранения основной системы (1.1) на гладких решениях, но может противоречить им при наличии сильных разрывов – ударных волн. В моделях механики деформируемых сред уравнение (1.2) представляет собой дифференциальную форму закона сохранения энергии. Обычно производящий потенциал $\Phi (U)$ зависит от дополнительного параметра (температуры или энтропии), который остается постоянным в области гладкости решения и меняется скачкообразно, обеспечивая выполнение этого закона, при переходе через поверхности разрывов. Но если в частных моделях влияние дополнительного параметра не учитывается (если рассматриваются изотермические или изэнтропические процессы), то невыполнение закона на разрывах вполне обосновано, так как должен происходить отток энергии за фронтом ударной волны для обеспечения изотермичности или изэнтропичности.

При условии непрерывной дифференцируемости производящих потенциалов уравнения основной системы преобразуются к недивергентной матричной форме:

(1.3)

$A(U)\frac{{\partial U}}{{\partial t}} = \sum\limits_{i = 1}^n \,{{B}^{i}}(U)\frac{{\partial U}}{{\partial {{x}_{i}}}} + g(U),\quad A(U) = \frac{{\partial \varphi }}{{\partial U}},\quad {{B}^{i}}(U) = \frac{{\partial {{\psi }_{i}}}}{{\partial U}},$

с симметричными матрицами–коэффициентами $A$ и ${{B}^{i}}$. Если матрица $A(U)$ положительно определена, что выполняется в случае сильно выпуклого производящего потенциала $\Phi (U)$, то рассматриваемая система относится к классу систем уравнений $t$-гиперболических по Фридрихсу [5]. Гиперболичность в известном смысле гарантирует корректность математической модели и позволяет применить к ее анализу хорошо разработанные вычислительные методы и технологии.

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

2. ВАРИАЦИОННЫЕ НЕРАВЕНСТВА

В качестве обобщения системы уравнений (1.1) будем рассматривать вариационное неравенство

(2.1)

$(\tilde {U} - U) \cdot (\mathcal{D}\left\langle U \right\rangle - g(U)) \geqslant 0,\quad U,\tilde {U} \in F.$

Здесь $F$ – выпуклое и замкнутое в $m$-мерном арифметическом пространстве множество допустимых вариаций вектора $U$, которое, вообще говоря, может параметрически зависеть от $t$ и $x$; $\tilde {U}$ – произвольный элемент $F$. В механике упругопластических сред вариационное неравенство представляет собой формулировку принципа Мизеса, в соответствии с которым при заданной скорости пластической деформации мощность диссипации принимает максимальное значение на действительных напряжениях. Граница множества $F$ в пространстве напряжений является поверхностью текучести материала. Если вектор $U$ лежит внутри $F$, то в силу произвольности вариации из (2.1) следует система уравнений (1.1), описывающая упругий процесс. Если $U$ – граничная точка, то выполняется ассоциированный закон пластического течения, в соответствии с которым вектор пластической скорости деформации, равный разности $g(U) - \mathcal{D}\left\langle U \right\rangle $, направлен по внешней нормали к границе $F$.

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

Обобщение принципа максимума Мизеса на произвольные термодинамически необратимые процессы, известное как фундаментальный принцип Циглера, играет исключительно важную роль в физике, химии и биологии (см. [6]). Его можно также использовать при моделировании условий на разрывах решений в задачах неравновесной термодинамики.

Приведем некоторые примеры моделей механики деформируемых сред, допускающих формулировку (2.1).

2.1. Теория упругопластического течения Прандтля–Рейсса

Вариационное неравенство в этой теории записывается относительно вектора скорости $v$ и симметричного тензора напряжений $\sigma $ [7]:

$(\tilde {v} - v) \cdot (\rho \dot {v} - \nabla \cdot \sigma - \rho g) + (\tilde {\sigma } - \sigma ):(a:\dot {\sigma } - \nabla v) \geqslant 0,\quad \sigma ,\tilde {\sigma } \in F.$

Здесь $\rho $ – плотность среды, $a$ – положительно-определенный тензор модулей упругой податливости четвертого ранга, обладающий специальной симметрией. Используются общепринятые обозначения и операции тензорного анализа.

В силу произвольности вариации вектора скорости, из этого неравенства следует система уравнений движения. Условие неотрицательности слагаемого с вариацией тензора напряжений, которое вытекает из него при $\tilde {v} = v$, в точности соответствует формулировке принципа Мизеса. Компонентами вектора $U$ при переходе к вариационному неравенству (2.1) являются компоненты вектора скорости и тензора напряжений в декартовой системе координат, записанные в столбец. В качестве производящих потенциалов берутся квадратичные функции (нижний индекс означает проекцию вектора на координатную ось):

(2.2)

$\Phi (U) = \frac{1}{2}(\rho {{v}^{2}} + \sigma :a:\sigma ),\quad {{\Psi }_{i}}(U) = {{(\sigma \cdot v)}_{i}}.$

В этом примере $\mathcal{D}\left\langle U \right\rangle $ является линейным дифференциальным оператором с постоянными коэффициентами.

Аналогичное неравенство описывает динамические процессы в структурно неоднородных материалах в рамках теории моментного континуума Коссера, в которой учитываются вращательные степени свободы частиц микроструктуры. Соответствующий вектор $U$ включает в себя компоненты вектора линейной скорости, несимметричного тензора напряжений, вектора угловой скорости и тензора моментных напряжений, также несимметричного, см. [8]. Более общее вариационное неравенство для нелинейного оператора возникает в теории упрочняющейся упругопластической среды. В этом случае вектор $U$, кроме компонент вектора скорости и тензора напряжений, содержит компоненты тензорного параметра трансляционного упрочнения и скалярный параметр изотропного упрочнения [7].

2.2. Модель упругопластической сыпучей среды

Определяющие уравнения для описания деформации сыпучей среды, по-разному сопротивляющейся растяжению и сжатию, строятся с помощью обобщенного реологического метода [9]. Математическая модель идеальной среды с упругопластическими частицами приводится к следующему вариационному неравенству:

$\begin{gathered} (\tilde {v} - v) \cdot (\rho \dot {v} - \nabla \cdot \sigma - \rho g) + (\tilde {\sigma } - \sigma ):(a:\dot {s} - \nabla {v}) \geqslant 0, \\ \sigma ,\tilde {\sigma } \in F,\quad \sigma = {{\pi }_{K}}(s). \\ \end{gathered} $

Здесь $s$ – симметричный тензор условных напряжений; ${{\pi }_{K}}$ – проектор на выпуклый и замкнутый конус $K$ с вершиной в нуле, который включает в себя всевозможные допустимые тензоры напряжений, отвечающие состоянию сжатия. Проекция на конус вычисляется по евклидовой норме ${\text{|}}\sigma {\kern 1pt} {{{\text{|}}}_{a}} = \sqrt {\sigma :a:\sigma } $.

В идеальной сыпучей среде отсутствуют связи между частицами, поэтому напряженные состояния, кроме состояний сжатия, в ней не реализуются. Такое свойство описывает замыкающее уравнение для определения тензора $\sigma $ через $s$.

Обозначим через $U$ и $V$ векторы с компонентами $v$, $s$ и с компонентами $v$, $\sigma $ соответственно. В этих обозначениях вариационное неравенство модели преобразуется к матричной форме:

(2.3)

$(\tilde {V} - V) \cdot \left( {A\frac{{\partial U}}{{\partial t}} - \sum\limits_{i = 1}^n \,{{B}^{i}}\frac{{\partial V}}{{\partial {{x}_{i}}}}} \right) \geqslant 0,\quad V,\tilde {V} \in F,\quad V = {{\pi }_{K}}(U),$

где $A$ и ${{B}^{i}}$ – симметричные матрицы с постоянными коэффициентами, которые могут быть получены двукратным дифференцированием по $U$ квадратичных потенциалов (2.2) теории упругопластического течения; ${{\pi }_{K}}$ – проектор на конус $K$ по норме, ассоциированной с положительно-определенной матрицей $A$.

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

Регуляризация состоит в замене входящего в (2.3) замыкающего уравнения уравнением с малым параметром $\varepsilon > 0$ (коэффициентом упругости связей):

$V = \varepsilon U + (1 - \varepsilon ){{\pi }_{K}}(U).$

По свойству проекции на выпуклый конус с вершиной в нуле (см., например, [9, гл. 2]) имеем ${{\pi }_{K}}(V) = {{\pi }_{K}}(U)$. Поэтому последнее уравнение можно представить в обращенной форме:

$U = \frac{1}{\varepsilon }V - \frac{{1 - \varepsilon }}{\varepsilon }{{\pi }_{K}}(V).$

Далее, полагая

(2.4)

$\Phi (V) = \frac{1}{{2\varepsilon }}(V \cdot AV - (1 - \varepsilon ){{\pi }_{K}}(V) \cdot A{{\pi }_{K}}(V)),\quad {{\Psi }_{i}}(V) = \frac{1}{2}V \cdot {{B}^{i}}V,$

приводим модель к вариационному неравенству (2.1) относительно вектор-функции $V(t,x)$.

Действительно, установленные в [9] свойства проекции на конус позволяют доказать, что производящий потенциал $\Phi $ является непрерывно-дифференцируемой функцией и что

$\frac{{\partial \Phi }}{{\partial V}} = \frac{1}{\varepsilon }AV - \frac{{1 - \varepsilon }}{\varepsilon }A{{\pi }_{K}}(V) = AU.$

Кроме того, так как проекция на конус с вершиной в нуле по норме ${\text{|}}V{\kern 1pt} {{{\text{|}}}_{A}} = \sqrt {V \cdot AV} $ удовлетворяет уравнению ${{\pi }_{K}}(V) \cdot A(V - {{\pi }_{K}}(V)) = 0$, то справедлива цепочка равенств:

$\begin{gathered} V \cdot AV + \frac{{1 - \varepsilon }}{\varepsilon }\left( {V \cdot AV - {{\pi }_{K}}(V) \cdot A{{\pi }_{K}}(V)} \right)\;{\text{ = }}\;{\text{|}}V{\kern 1pt} {\text{|}}_{A}^{2} + \frac{{1 - \varepsilon }}{\varepsilon }(V + {{\pi }_{K}}(V)) \cdot A(V - {{\pi }_{K}}(V)) = \\ \, = {\text{|}}V{\kern 1pt} {\text{|}}_{A}^{2} + \frac{{1 - \varepsilon }}{\varepsilon }\left( {{\text{|}}V - {{\pi }_{K}}(V){\kern 1pt} {\text{|}}_{A}^{2} + 2{{\pi }_{K}}(V) \cdot A(V - {{\pi }_{K}}(V))} \right)\;{\text{ = }}\;{\text{|}}V{\kern 1pt} {\text{|}}_{A}^{2} + \frac{{1 - \varepsilon }}{\varepsilon }{\text{|}}V - {{\pi }_{K}}(V){\kern 1pt} {\text{|}}_{A}^{2}, \\ \end{gathered} $

из которой следует, что $\Phi (V)$ – сильно выпуклая функция при $\varepsilon \leqslant 1$. Это гарантирует гиперболичность дифференциального оператора $\mathcal{D}\left\langle V \right\rangle $.

2.3. Модель пористой среды

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

$\begin{gathered} (\tilde {v} - v) \cdot (\rho \dot {v} - \nabla \cdot \sigma - \rho g) + (\tilde {s} - s):(a:\dot {s} - \nabla v) + (\tilde {q} - {{\pi }_{K}}(q)):(b:\dot {q} - \nabla v) \geqslant 0, \\ s,\tilde {s} \in F,\quad \sigma = s + {{\pi }_{K}}(q). \\ \end{gathered} $

Здесь $a$ – симметричный и положительно-определенный тензор четвертого ранга модулей упругой податливости пористого скелета; $b$ – тензор с такими же свойствами, характеризующий повышение жесткости материала при схлопывании пор. Симметричные тензоры напряжений $s$ и $q$, через которые определяется тензор действительных напряжений $\sigma $, описывают напряженное состояние пористого скелета и изменение этого состояния за счет контактного взаимодействия перегородок скелета после схлопывания пор.

В зависимости от пористости среды, пластичность может наступать до схлопывания пор (как правило, это происходит в высокопористых материалах) или после него. Как и в предыдущих моделях, переход в пластическое состояние описывается с помощью выпуклого и замкнутого множества $F$ в пространстве напряжений, ограниченного поверхностью текучести материала. Начальная пористость, которая в общем случае неоднородно распределена по объему среды, учитывается при постановке начальных данных.

Составим вектор $U$ из компонент вектора скорости и тензоров напряжений $s$ и $q$, а в векторе $V$ вместо тензора $q$ будем использовать его проекцию ${{\pi }_{K}}(q)$ на конус $K$. Тогда эти векторы свяжутся уравнением $V = {{\pi }_{K}}(U)$, а вариационное неравенство модели запишется в виде (2.3) с матрицами–коэффициентами соответствующей размерности. Таким образом, в общих обозначениях рассматриваемая модель пористой среды полностью повторяет модель идеальной сыпучей среды. Следовательно, переходя к регуляризованному замыкающему уравнению с малым параметром $\varepsilon > 0$:

$\sigma = s + \varepsilon q + (1 - \varepsilon ){{\pi }_{K}}(q),$

можно привести ее к вариационному неравенству (2.1) с производящими потенциалами (2.4).

Заметим, что если в приведенных выше моделях учесть вращательное движение частиц, то получатся новые, уточненные математические модели структурно неоднородных сыпучих и пористых материалов, которые формулируются в виде (2.1).

3. ИНТЕГРАЛЬНЫЕ ОЦЕНКИ

Рассмотрим некоторые общие свойства решений вариационного неравенства (2.1). Пусть $U(t,x)$ и $U{\kern 1pt} {\text{'}}(t,x)$ – непрерывно-дифференцируемые решения, определенные в замкнутой области $G$. Ниже будет предполагаться, что производящие потенциалы $\Phi (U)$ и ${{\Psi }_{i}}(U)$ являются трижды непрерывно-дифференцируемыми функциями, а вектор-функция $g(U)$ удовлетворяет условию Липшица по $U$.

Полагая $\tilde {U} = U{\kern 1pt} {\text{'}} \in F$ в (2.1), $\tilde {U} = U \in F$ в вариационном неравенстве, характеризующем решение $U{\kern 1pt} {\text{'}}$, и суммируя результаты, получаем

(3.1)

$(U{\kern 1pt} {\text{'}} - U) \cdot \left( {\mathcal{D}\left\langle {U{\kern 1pt} {\text{'}}} \right\rangle - \mathcal{D}\left\langle U \right\rangle } \right) \leqslant (U{\kern 1pt} {\text{'}} - U) \cdot (g{\text{'}}(U{\kern 1pt} {\text{'}}) - g(U)).$

Чтобы оценить влияние вектора правой части на решение, рассмотрим общий случай, когда вектор-функция $g{\text{'}}$, входящая в вариационное неравенство для $U{\kern 1pt} {\text{'}}$, отличается от вектор-функции $g$ в (2.1), и когда обе они зависят от $t$ и $x$.

Неравенство, аналогичное (3.1), является отправной точкой при выводе априорных оценок интегральной нормы разности решений гиперболической системы квазилинейных уравнений (1.3) в характеристических коноидах дифференциального оператора. Оценки получены, например, в [2], [11], [12] и в точности воспроизводятся для решений вариационного неравенства.

Предположим, что $G$ – усеченный коноид в пространстве переменных $(t,x)$, основаниями которого служат две гиперплоскости $t = {{t}_{0}}$ и $t = {{t}_{1}}$ (фиг. 1), а боковая поверхность ${{S}_{c}}$, уравнение которой имеет вид $h(t,x) = 0$, удовлетворяет неравенству Гамильтона–Якоби для оператора $\mathcal{D}\left\langle U \right\rangle $:

$\dot {h} + H(\nabla h) \geqslant 0.$

Здесь $H(\nu )$ – наименьший среди $m$ действительных корней $c = {{H}_{k}}(\nu )$ $(k = 1,...,m)$ характеристического уравнения

$detD(c,\nu ) = 0,\quad D(c,\nu ) = cA(U) + \sum\limits_{i = 1}^n \,{{\nu }_{i}}{{B}^{i}}(U).$

Важнейшее свойство неравенства Гамильтона–Якоби состоит в том, что в точках построенной таким образом конической поверхности ${{S}_{c}}$ выполняется условие неотрицательной определенности матрицы $D(\dot {h}, - \nabla h)$. Так как матрица $A(U)$ положительно определена, то для выполнения последнего условия необходимо и достаточно, чтобы были неотрицательны все корни следующего многочлена от $\lambda $:

$\det (D(\dot {h}, - \nabla h) - \lambda A(U)) = detD(\dot {h} - \lambda , - \nabla h).$

Корни этого многочлена легко вычисляются с учетом принятых обозначений:

${{\lambda }_{k}} = \dot {h} - {{H}_{k}}( - \nabla h) = \dot {h} + {{H}_{k}}(\nabla h),\quad k = 1,...,m,$

а условие неотрицательности принимает вид неравенства Гамильтона–Якоби.

Фиг. 1.

Усеченный коноид для оценки решений.

Способ построения конической области $G$ с такими свойствами внутри произвольной окрестности, принадлежащей области определения решений, на основе метода характеристик для уравнения Гамильтона–Якоби детально описан в [2, с. 162–180].

Определим норму разности решений через интеграл от квадратичной формы по сечению $\Omega (t)$ области $G$ гиперплоскостью $t = {\text{const}}$:

${{\left\| {U{\kern 1pt} {\text{'}} - U} \right\|}^{2}}(t) = \int\limits_{\Omega (t)} \,(U{\kern 1pt} {\text{'}} - U) \cdot A(U)(U{\kern 1pt} {\text{'}} - U)d\Omega .$

Для получения оценки преобразуем левую часть (3.1) к виду

$\frac{1}{2}\frac{\partial }{{\partial t}}((U{\kern 1pt} {\text{'}} - U) \cdot A(U)(U{\kern 1pt} {\text{'}} - U)) - \frac{1}{2}\frac{\partial }{{\partial {{x}_{i}}}}\sum\limits_{i = 1}^n \,((U{\kern 1pt} {\text{'}} - U) \cdot {{B}^{i}}(U)(U{\kern 1pt} {\text{'}} - U)) - (U{\kern 1pt} {\text{'}} - U) \cdot \mathcal{N}\left\langle {U{\kern 1pt} {\text{'}},U} \right\rangle ,$

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

$\mathcal{N}\left\langle {U{\kern 1pt} {\text{'}},U} \right\rangle = \frac{1}{2}\left( {\frac{\partial }{{\partial t}}A(U) - \sum\limits_{i = 1}^n \,\frac{\partial }{{\partial {{x}_{i}}}}{{B}^{i}}(U)} \right)(U{\kern 1pt} {\text{'}} - U) - (A(U{\kern 1pt} {\text{'}}) - A(U))\frac{{\partial U{\kern 1pt} {\text{'}}}}{{\partial t}} + \sum\limits_{i = 1}^n \,({{B}^{i}}(U{\kern 1pt} {\text{'}}) - {{B}^{i}}(U))\frac{{\partial U{\kern 1pt} {\text{'}}}}{{\partial {{x}_{i}}}}.$

После таких преобразований проинтегрируем обе части неравенства по области $G$, применяя формулу Грина к дивергентным слагаемым. Это дает

$\begin{gathered} {{\left\| {U{\kern 1pt} {\text{'}} - U} \right\|}^{2}}({{t}_{1}}) - {{\left\| {U{\kern 1pt} {\text{'}} - U} \right\|}^{2}}({{t}_{0}}) - \int\limits_{{{S}_{c}}} \,(U{\kern 1pt} {\text{'}} - U) \cdot D(c,\nu )(U{\kern 1pt} {\text{'}} - U)\frac{{dS}}{{\sqrt {1 + {{c}^{2}}} }} \leqslant \\ \, \leqslant 2\iint\limits_G {(U{\kern 1pt} {\text{'}} - U) \cdot (\mathcal{N}\left\langle {U{\kern 1pt} {\text{'}},U} \right\rangle + g{\text{'}}(U{\kern 1pt} {\text{'}}) - g(U))d\Omega dt.} \\ \end{gathered} $

Здесь $c = - \dot {h}{\text{/|}}\nabla h{\kern 1pt} {\text{|}}$ – скорость движения границы области $\Omega (t)$ в направлении нормали $\nu = \nabla h{\text{/|}}\nabla h{\kern 1pt} {\text{|}}$. Единичный вектор с компонентами $( - c,\nu ){\text{/}}\sqrt {1 + {{c}^{2}}} $ является внешней по отношению к $G$ нормалью к конической поверхности ${{S}_{c}}$ в пространстве переменных $t,\;x$.

Правую часть полученного неравенства разложим на слагаемые, полагая в ней

$g{\text{'}}(U{\kern 1pt} {\text{'}}) - g(U) = (g{\text{'}}(U{\kern 1pt} {\text{'}}) - g{\text{'}}(U)) + (g{\text{'}}(U) - g(U)).$

После этого разность в первой скобке оценивается по условию Липшица, вторая скобка не преобразуется, так как именно она характеризует изменение $g$. Для интеграла в правой части в силу неравенства Коши–Буняковского справедлива оценка

$\int\limits_{{{t}_{0}}}^{{{t}_{1}}} \,\left( {a{{{\left\| {U{\kern 1pt} {\text{'}} - U} \right\|}}^{2}} + b\left\| {g{\text{'}}(U{\kern 1pt} {\text{'}}) - g(U{\text{'}})} \right\|\left\| {U{\kern 1pt} {\text{'}} - U} \right\|} \right)(t)dt$

с положительными постоянными $a$ и $b$, зависящими, вообще говоря, от обоих решений, их первых производных по $t$ и ${{x}_{i}}$, а также от константы Липшица для вектор-функции $g(U)$. Как результат, с учетом неположительности подынтегрального выражения в интеграле по ${{S}_{c}}$ получим неравенство:

${{\left\| {U{\kern 1pt} {\text{'}} - U} \right\|}^{2}}({{t}_{1}}) \leqslant {{\left\| {U{\kern 1pt} {\text{'}} - U} \right\|}^{2}}({{t}_{0}}) + 2a\int\limits_{{{t}_{0}}}^{{{t}_{1}}} \,{{\left\| {U{\kern 1pt} {\text{'}} - U} \right\|}^{2}}(t)dt + 2b\int\limits_{{{t}_{0}}}^{{{t}_{1}}} \,\left( {\left\| {g{\text{'}}(U{\kern 1pt} {\text{'}}) - g(U{\kern 1pt} {\text{'}})} \right\|\left\| {U{\kern 1pt} {\text{'}} - U} \right\|} \right)(t)dt.$

Это же неравенство остается справедливым, если заменить в нем ${{t}_{0}}$ и ${{t}_{1}}$ на произвольные моменты времени $t{\text{'}} < t$, принадлежащие интервалу $({{t}_{0}},{{t}_{1}})$. При стремлении $t{\text{'}}$ к $t$ полученное неравенство преобразуется к более простому виду:

$\frac{d}{{dt}}\left\| {U{\kern 1pt} {\text{'}} - U} \right\|(t) \leqslant a\left\| {U{\kern 1pt} {\text{'}} - U} \right\|(t) + b\left\| {g{\text{'}}(U{\kern 1pt} {\text{'}}) - g(U{\kern 1pt} {\text{'}})} \right\|(t)$

или

${{e}^{{at}}}\frac{d}{{dt}}\left( {{{e}^{{ - at}}}\left\| {U{\kern 1pt} {\text{'}} - U} \right\|(t)} \right) \leqslant b\left\| {g{\text{'}}(U{\kern 1pt} {\text{'}}) - g(U{\kern 1pt} {\text{'}})} \right\|(t),$

откуда после интегрирования по $t$ от ${{t}_{0}}$ до ${{t}_{1}}$ следует окончательная оценка:

(3.2)

$\left\| {U{\kern 1pt} {\text{'}} - U} \right\|({{t}_{1}}) \leqslant \left\| {U{\kern 1pt} {\text{'}} - U} \right\|({{t}_{0}}){{{\text{e}}}^{{a({{t}_{1}} - {{t}_{0}})}}} + b\int\limits_{{{t}_{0}}}^{{{t}_{1}}} \,\left\| {g{\text{'}}(U{\kern 1pt} {\text{'}}) - g(U{\kern 1pt} {\text{'}})} \right\|(t){{{\text{e}}}^{{a({{t}_{1}} - t)}}}dt.$

Из оценки (3.2) вытекает утверждение о единственности и непрерывной зависимости по начальным данным решения задачи Коши:

$U{{{\text{|}}}_{{t = 0}}} = {{U}_{0}}(x) \in F.$

Действительно, если одно из решений задачи Коши есть $U(t,x)$, то в силу (3.2) всякое другое решение этой же задачи $U{\kern 1pt} {\text{'}}(t,x)$ совпадает с ним в любом усеченном коноиде $G$, боковая поверхность которого удовлетворяет неравенству Гамильтона–Якоби. Кроме того, малое изменение вектор-функции ${{U}_{0}}(x)$, как и малое изменение вектор-функции $g(U)$, приводит к малому изменению решения в любом сечении $\Omega (t)$ области $G$ гиперплоскостью $t = {\text{const}}$.

На основе (3.2) доказывается также ограниченность скорости распространения возмущений (конечность областей зависимости и влияния решений).

Аналогичную оценку можно получить в окрестности неподвижной гиперповерхности ${{S}_{b}}$, на которой заданы диссипативные граничные условия, т.е. условия общего вида, выполнение которых для вектор-функций $U$ и $U{\kern 1pt} {\text{'}}$ обеспечивает выполнение в точках ${{S}_{b}}$ неравенства

$(U{\kern 1pt} {\text{'}} - U) \cdot \sum\limits_{i = 1}^n \,{{\nu }_{i}}{{B}^{i}}(U)(U{\kern 1pt} {\text{'}} - U) \leqslant 0,$

где $\nu $ – вектор внешней по отношению к области решения задачи нормали к сечению ${{S}_{b}}$ гиперплоскостью $t = {\text{const}}$. При выводе оценки в качестве области $G$ нужно взять часть усеченного коноида, опирающуюся на ${{S}_{b}}$ (см. фиг. 2).

Фиг. 2.

Усеченный коноид, опирающийся на границу области.

4. РАЗРЫВНЫЕ РЕШЕНИЯ

Вариационное неравенство (2.1) записывается в дивергентной форме, на основе которой можно строить обобщенные (разрывные) решения:

$\tilde {U} \cdot \mathcal{D}\left\langle U \right\rangle - (\tilde {U} - U) \cdot g(U) \geqslant \frac{\partial }{{\partial t}}(U \cdot \varphi (U) - \Phi (U)) - \sum\limits_{i = 1}^n \,\frac{\partial }{{\partial {{x}_{i}}}}(U \cdot {{\psi }_{i}}(U) - {{\Psi }_{i}}(U)),$

и имеет, таким образом, интегральное обобщение, эквивалентное ему на гладких решениях:

(4.1)

$\begin{gathered} \iint\limits_G {\left( { - \frac{{\partial (\chi \tilde {U})}}{{\partial t}} \cdot \varphi (U) + \sum\limits_{i = 1}^n \,\frac{{\partial (\chi \tilde {U})}}{{\partial {{x}_{i}}}} \cdot {{\psi }_{i}}(U) - \chi (\tilde {U} - U) \cdot g(U)} \right)d\Omega dt} \geqslant \\ \, \geqslant \iint\limits_G {\left( { - \frac{{\partial \chi }}{{\partial t}}(U \cdot \varphi (U) - \Phi (U)) + \sum\limits_{i = 1}^n \,\frac{{\partial \chi }}{{\partial {{x}_{i}}}}(U \cdot {{\psi }_{i}}(U) - {{\Psi }_{i}}(U))} \right)}{\kern 1pt} d\Omega dt. \\ \end{gathered} $

Интегральное обобщение получается после умножения обеих частей дивергентного неравенства на произвольную финитную в $G$ неотрицательную функцию $\chi \in {{C}^{\infty }}(G)$ и последующего интегрирования по области $G$ с применением формулы Грина ко всем слагаемым, содержащим производные от решения по времени и пространственным переменным. Предполагается, что “пробная” вектор-функция $\tilde {U}(t,x)$ в неравенстве (4.1) непрерывно-дифференцируема в области $G$ и поточечно удовлетворяет ограничению $\tilde {U} \in F$.

Интегральное неравенство естественным образом определяет класс обобщенных решений. К этому классу относятся такие вектор-функции $U \in {{L}_{1}}(G)$, удовлетворяющие (4.1) при всех допустимых $\tilde {U}$, для которых входящие в это неравенство интегралы определены в смысле Лебега. Сюда, в частности, подпадают решения с сильным разрывом, имеющие разрыв I рода на некоторой гиперповерхности ${{S}_{0}} \subset G$ и непрерывно-дифференцируемые в оставшейся части области $G$.

В результате применения формулы Грина к интегралам по подобластям непрерывности такого решения, из (4.1) получается система соотношений, выполненных в точках гиперповерхности ${{S}_{0}}$. При этом входящие в (4.1) интегралы, содержащие производные по времени и по пространственным переменным

$\iint\limits_{{{G}^{ \pm }}} {\frac{{\partial (\chi \tilde {U})}}{{\partial t}} \cdot \varphi (U)d\Omega dt},\quad \iint\limits_{{{G}^{ \pm }}} {\frac{{\partial (\chi \tilde {U})}}{{\partial {{x}_{i}}}} \cdot {{\psi }_{i}}(U)d\Omega dt},$

преобразуются как

$ \mp \int\limits_{{{S}_{0}}} \,c\tilde {U} \cdot \varphi (U)\frac{{\chi dS}}{{\sqrt {1 + {{c}^{2}}} }} - \iint\limits_{{{G}^{ \pm }}} {\tilde {U} \cdot \frac{{\partial \varphi }}{{\partial t}}\chi d\Omega dt},\quad \pm {\kern 1pt} \int\limits_{{{S}_{0}}} \,{{\nu }_{i}}\tilde {U} \cdot {{\psi }_{i}}(U)\frac{{\chi dS}}{{\sqrt {1 + {{c}^{2}}} }} - \iint\limits_{{{G}^{ \pm }}} {\tilde {U}} \cdot \frac{{\partial {{\psi }_{i}}(U)}}{{\partial {{x}_{i}}}}\chi d\Omega dt.$

Остальные интегралы преобразуются аналогично. После сложения интегралов по подобластям ${{G}^{ + }}$ и ${{G}^{ - }}$ получается неравенство

$\begin{array}{*{20}{c}} {\int\limits_{{{S}_{0}}} \,\tilde {U} \cdot (c\varphi (U) + \sum\limits_{i = 1}^n \,{{\nu }_{i}}{{\psi }_{i}}(U))\frac{{\chi dS}}{{\sqrt {1 + {{c}^{2}}} }} + \iint\limits_{{{G}^{ + }} \cup {{G}^{ - }}} {\tilde {U} \cdot \left( {\frac{{\partial \varphi (U)}}{{\partial t}} - \sum\limits_{i = 1}^n \,\frac{{\partial {{\psi }_{i}}(U)}}{{\partial {{x}_{i}}}}} \right)}{\kern 1pt} \chi d\Omega dt - } \\ { - \;\iint\limits_G {(\tilde {U} - U) \cdot g(U)\chi d\Omega dt \geqslant }\int\limits_{{{S}_{0}}} \,\left( {c[U \cdot \varphi (U) - \Phi (U)] + \sum\limits_{i = 1}^n \,{{\nu }_{i}}[U \cdot {{\psi }_{i}}(U) - {{\Psi }_{i}}(U)]} \right)\frac{{\chi dS}}{{\sqrt {1 + {{c}^{2}}} }} + } \\ {\, + \iint\limits_{{{G}^{ + }} \cup {{G}^{ - }}} {\left( {\frac{{\partial (U \cdot \varphi (U) - \Phi (U))}}{{\partial t}} - \sum\limits_{i = 1}^n \,\frac{{\partial (U \cdot {{\psi }_{i}}(U) - {{\Psi }_{i}}(U))}}{{\partial {{x}_{i}}}}} \right)\chi d\Omega dt}.} \end{array}$

Здесь квадратные скобки означают скачок функции на разрыве: $[a] = {{a}^{ + }} - {{a}^{ - }}$, где ${{a}^{ \pm }}$ – односторонние пределы на ${{S}_{0}}$ со стороны ${{G}^{ + }}$ и ${{G}^{ - }}$ соответственно.

Пусть $({{t}_{0}},{{x}^{0}}) \in G$ – точка непрерывности решения $U(t,x)$. Тогда существует окрестность этой точки, не пересекающая гиперповерхность ${{S}_{0}}$. Для всякой допустимой “пробной” функции $\chi (t,x)$, носитель которой сосредоточен в данной окрестности, интегралы по ${{S}_{0}}$ в полученном неравенстве равны нулю, следовательно,

$\iint\limits_G {(\tilde {U} - U)} \cdot (\mathcal{D}\left\langle U \right\rangle - g(U))\chi d\Omega dt \geqslant 0.$

Если же $({{t}_{0}},{{x}^{0}}) \in {{S}_{0}}$, то для функции $\chi (t,x)$, носитель которой сосредоточен в окрестности этой точки, имеем

$\begin{array}{*{20}{c}} {\int\limits_{{{S}_{0}}} \,(\tilde {U}(c\varphi (U) + \sum\limits_{i = 1}^n \,{{\nu }_{i}}{{\psi }_{i}}(U)) - c[U \cdot \varphi (U) - \Phi (U)] - \sum\limits_{i = 1}^n \,{{\nu }_{i}}[U \cdot {{\psi }_{s}}(U) - {{\Psi }_{i}}(U)])\frac{{\chi dS}}{{\sqrt {1 + {{c}^{2}}} }} \geqslant } \\ { \geqslant \; - {\kern 1pt} \iint\limits_{{{G}^{ + }} \cup {{G}^{ - }}} {(\tilde {U} - U)} \cdot (\mathcal{D}\left\langle U \right\rangle - g(U))\chi d\Omega dt.} \end{array}$

В силу произвольности функции $\chi \geqslant 0$ и с учетом возможности неограниченного уменьшения рассматриваемой окрестности из этих двух неравенств следует, что в открытых подобластях ${{G}^{ \pm }}$, где решение достаточно гладкое, выполняется вариационное неравенство (2.1), а на ${{S}_{0}}$ справедливо неравенство

$\tilde {U} \cdot [r(U)] \geqslant c[U \cdot \varphi (U) - \Phi (U)] + \sum\limits_{i = 1}^n \,{{\nu }_{i}}[U \cdot {{\psi }_{i}}(U) - {{\Psi }_{i}}(U)],\quad r(U) = c\varphi (U) + \sum\limits_{i = 1}^n \,{{\nu }_{i}}{{\psi }_{i}}(U),$

в котором $c > 0$ – скорость фронта разрыва, представляющего собой сечение ${{S}_{0}}$ гиперплоскостью $t = {\text{const}}$, в направлении нормали $\nu $, внешней по отношению к ${{G}^{ + }}$.

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

${{a}^{ + }}{{b}^{ + }} - {{a}^{ - }}{{b}^{ - }} = ({{a}^{ + }} - {{a}^{ - }})\frac{{{{b}^{ + }} + {{b}^{ - }}}}{2} + \frac{{{{a}^{ + }} + {{a}^{ - }}}}{2}({{b}^{ + }} - {{b}^{ - }}),$

неравенство в точках гиперповерхности разрыва решения переписывается в эквивалентной форме:

(4.2)

$(\tilde {U} - {{U}^{0}})[r(U)] \geqslant - d(U) = c([U] \cdot {{\varphi }^{0}}(U) - [\Phi (U)]) + \sum\limits_{i = 1}^n \,{{\nu }_{i}}([U] \cdot \psi _{i}^{0}(U) - [{{\Psi }_{i}}(U)]),$

где ${{U}^{0}} = ({{U}^{ + }} + {{U}^{ - }}){\text{/}}2$, величины ${{\varphi }^{0}}(U)$ и $\psi _{i}^{0}(U)$ определяются аналогично.

Полагая в (4.2) $\tilde {U} = {{U}^{0}} \in F$, можно получить условие реализуемости разрыва: $d(U) \geqslant 0$. Это же условие возникает при анализе разрывных решений системы уравнений (1.1), являющейся частным случаем вариационного неравенства (2.1), когда множество $F$ совпадает со всем пространством. В моделях механики деформируемых сред величина $d(U)$ представляет собой часть энергии, которая выделяется на фронте ударной волны за счет упругой деформации. Левая часть (4.2) характеризует энергию, выделяемую за счет пластичности. На основе гипотезы о независимости термодинамически обратимых и необратимых процессов, будем считать, что вместо (4.2) на гиперповерхности разрыва выполняется более строгое неравенство:

(4.3)

$(\tilde {U} - {{U}^{0}})[r(U)] \geqslant 0,\quad {{U}^{ \pm }},\tilde {U} \in F.$

На самом деле это неравенство соответствует формулировке принципа Циглера для термодинамически необратимого процесса перехода среды из состояния перед фронтом ударной волны в состояние за ее фронтом, в котором вектор ${{U}^{0}}$ играет роль вектора термодинамических потоков, а $[r(U)]$ – роль вектора термодинамических сил. С этой точки зрения вариационное неравенство (4.2) представляет собой математическую модель сильного разрыва как самостоятельного явления.

Справедливо следующее утверждение: при выполнении (4.3) на фронте разрыва решения для всякого $\alpha :{\text{|}}\alpha {\kern 1pt} {\text{|}} \leqslant 1$ выполняется более общее вариационное неравенство:

(4.4)

$(\tilde {U} - {{U}^{\alpha }})[r(U)] \geqslant 0,\quad {{U}^{\alpha }} = \frac{{1 + \alpha }}{2}{{U}^{ + }} + \frac{{1 - \alpha }}{2}{{U}^{ - }},\quad \tilde {U} \in F.$

Действительно, в силу (4.3) это неравенство выполняется при $\alpha = 0$. Принимая в качестве $\tilde {U} \in F$ сначала вектор ${{U}^{ + }}$, а затем ${{U}^{ - }}$, и суммируя полученные из него неравенства, можно установить, что $[U] \cdot [r(U)] = 0$. Отсюда с учетом формулы ${{U}^{\alpha }} = {{U}^{0}} + \alpha [U]{\text{/}}2$ вытекает неравенство (4.4) при любом значении $\alpha $.

Геометрическая интерпретация доказанного утверждения состоит в следующем. Если какая-либо точка ${{U}^{\alpha }}$ отрезка в $m$-мерном пространстве с концами ${{U}^{ - }}$, ${{U}^{ + }}$ лежит строго внутри множества $F$, то в силу произвольности входящей в (4.4) вариации $\tilde {U} - {{U}^{\alpha }}$ вектор $[r(U)]$ равен нулю. Если же весь отрезок принадлежит границе выпуклого множества $F$, то этот вектор направлен по внутренней нормали к границе, причем направление нормали не меняется с переходом к другой точке отрезка.

Таким образом, обобщенные решения в механике упругопластических сред могут содержать разрывы двух видов: нейтральные (упругие) волны, определяемые системой уравнений $[r(U)] = 0$, и диссипативные (пластические) волны, отвечающие условию градиентности по отношению к поверхности текучести материала, для которых отрезок с концами ${{U}^{ + }}$, ${{U}^{ - }}$ принадлежит поверхности текучести.

Заметим, что в геометрически линейной теории упруго-идеальнопластической среды производящие потенциалы $\Phi (U)$ и ${{\Psi }_{i}}(U)$ являются квадратичными функциями, поэтому $d(U)$ тождественно равно нулю. Следовательно, неравенство (4.3) вытекает из (4.2) без каких-либо дополнительных гипотез и предположений.

С помощью описанного метода была получена полная система соотношений сильного разрыва в теории упруго-идеальнопластического течения [13]. Ранее было показано [14], что система уравнений этой теории, основанная на ассоциированном законе пластичности, неприводима к дивергентной форме. Поэтому ее невозможно обобщить в виде полной системы интегральных законов сохранения. Скорости пластических ударных волн впервые были получены [15] с применением вспомогательной гипотезы о максимальной диссипации энергии на разрыве. В модели линейного изотропного и трансляционного упрочнения, в которой производящие потенциалы также квадратичные, описанный метод был применен в работе [16]. Качественный анализ разрывных решений в рамках этой модели выполнен в [17]. В работах [18], [19] исследовались разрывные решения в теории сыпучей среды с учетом разного сопротивления материала растяжению и сжатию.

5. ВЫЧИСЛИТЕЛЬНЫЕ АЛГОРИТМЫ

Общий способ построения численных методов решения задач, допускающих постановку в виде вариационных неравенств, сводится к двум этапам: аппроксимации дифференциального оператора задачи и аппроксимации ограничения, определяющего множество допустимых вариаций решения. Как обычно, аппроксимация оператора предполагает его замену дискретным аналогом ${{\mathcal{D}}_{h}}\left\langle {{{U}_{h}}} \right\rangle $ ($h$ – характерный параметр пространственно-временнóй сетки). При аппроксимации ограничения фактически указывается, в каком смысле понимается включение искомого сеточного решения ${{U}_{h}}$ в множество $F$. После этого неравенству (2.1) приводится в соответствие дискретное вариационное неравенство

(5.1)

$({{\tilde {U}}_{h}} - {{\hat {U}}_{h}}) \cdot ({{\mathcal{D}}_{h}}\left\langle {{{U}_{h}}} \right\rangle - g({{U}_{h}})) \geqslant 0,\quad {{\hat {U}}_{h}} = {{\mathcal{I}}_{h}}\left\langle {{{U}_{h}}} \right\rangle ,\quad {{\tilde {U}}_{h}} \in F,$

где ${{\mathcal{I}}_{h}}\left\langle {{{U}_{h}}} \right\rangle $ – разностный оператор, осуществляющий аппроксимацию ограничения.

Пусть ${{\mathcal{D}}_{h}}\left\langle {{{U}_{h}}} \right\rangle $ – неконсервативный разностный оператор следующего вида:

${{\mathcal{D}}_{h}}\left\langle {{{U}_{h}}} \right\rangle = A({{U}^{{k - 1}}})\frac{{{{U}^{k}} - {{U}^{{k - 1}}}}}{{\Delta t}} - \sum\limits_{i = 1}^n \,{{B}^{i}}({{U}^{{k - 1}}}){{\mathcal{L}}_{i}}\left\langle {{{U}^{{k - 1}}}} \right\rangle - g({{U}^{{k - 1}}}).$

Здесь ${{\mathcal{L}}_{i}}\left\langle {{{U}^{{k - 1}}}} \right\rangle $ – явная по времени разностная аппроксимация производной по пространственной переменной; верхние индексы указывают на номер слоя по времени, для сокращения обозначений индекс $h$ в развернутой записи формул опускается.

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

Введем вспомогательную сеточную вектор-функцию $\mathop {\bar {U}}\nolimits_h $ как решение разностной схемы для системы квазилинейных уравнений $\mathcal{D}\left\langle U \right\rangle = g(U)$:

(5.2)

$A({{U}^{{k - 1}}})\frac{{{{{\bar {U}}}^{k}} - {{U}^{{k - 1}}}}}{{\Delta t}} = \sum\limits_{i = 1}^n \,{{B}^{i}}({{U}^{{k - 1}}}){{\mathcal{L}}_{i}}\left\langle {{{U}^{{k - 1}}}} \right\rangle + g({{U}^{{k - 1}}}).$

Тогда дискретное вариационное неравенство (5.1) в каждом узле сеточной области запишется в следующем виде:

(5.3)

$(\tilde {U} - {{\hat {U}}^{k}}) \cdot A({{U}^{{k - 1}}})({{U}^{k}} - {{\bar {U}}^{k}}) \geqslant 0,\quad {{\hat {U}}^{k}},\tilde {U} \in F.$

Используя произвол в выборе связи между узловыми значениями ${{\hat {U}}_{h}}$ и ${{U}_{h}}$, можно строить различные корректирующие алгоритмы для вычисления решения ${{U}^{k}}$ вариационного неравенства (5.1) через решение ${{\bar {U}}^{k}}$ системы уравнений (5.2). В простейшем случае, когда ${{\hat {U}}^{k}} = {{U}^{k}}$, решение (5.3) эквивалентно поиску минимума на множестве $F$ квадратичной функции $({{U}^{k}} - {{\bar {U}}^{k}}) \cdot A({{U}^{{k - 1}}})({{U}^{k}} - {{\bar {U}}^{k}})$, определяющей квадрат нормы, ассоциированной с симметричной положительно-определенной матрицей $A({{U}^{{k - 1}}})$, или, что то же самое, нахождению вектора ${{U}^{k}}$ как проекции ${{\pi }_{{k - 1}}}({{\bar {U}}^{k}})$ на множество $F$ по этой норме. Такой алгоритм при численной реализации модели изотропной упруго-идеальнопластической среды с условием пластичности Мизеса известен как алгоритм корректировки напряжений Уилкинса [20].

Задавая эту связь по формуле $2{{\hat {U}}^{k}} = (1 + \alpha ){{U}^{k}} + (1 - \alpha ){{\bar {U}}^{k}}$, можно из (5.3) получить более общий вариант корректировки, зависящий от параметра $\alpha $. В этом случае

(5.4)

${{U}^{k}} = \frac{2}{{1 + \alpha }}{{\hat {U}}^{k}} - \frac{{1 - \alpha }}{{1 + \alpha }}{{\bar {U}}^{k}},\quad {{\hat {U}}^{k}} = {{\pi }_{{k - 1}}}({{\bar {U}}^{k}}).$

Вариационное неравенство (5.3) позволяет исследовать вычислительную устойчивость данного алгоритма. Пусть $\delta {{\bar {U}}^{k}}$ – возмущение вспомогательного решения ${{\bar {U}}^{k}}$ из-за ошибок округления. Полагая $\tilde {U} = {{\hat {U}}^{k}} + \delta {{\hat {U}}^{k}}$ в (5.3) и $\tilde {U} = {{\hat {U}}^{k}}$ в аналогичном неравенстве, записанном для возмущенного решения ${{U}^{k}} + \delta {{U}^{k}}$, и суммируя результаты, можно установить, что

$\delta {{\hat {U}}^{k}} \cdot A({{U}^{{k - 1}}})(\delta {{U}^{k}} - \delta {{\bar {U}}^{k}}) \leqslant 0.$

Отсюда после вычисления вариации $\mathop {\hat {U}}\nolimits^k $ через вариации ${{U}^{k}}$ и $\mathop {\bar {U}}\nolimits^k $ следует неравенство:

${\text{|}}\delta {{U}^{k}}{\text{|}}_{{k - 1}}^{2} + \alpha {\text{|}}\delta {{U}^{k}} - \delta {{\bar {U}}^{k}}{\text{|}}_{{k - 1}}^{2} \leqslant {\text{|}}\delta {{\bar {U}}^{k}}{\text{|}}_{{k - 1}}^{2}$

(прямые скобки обозначают ассоциированную норму вектора). При $\alpha \geqslant 0$ это неравенство приводит к условию невозрастания нормы возмущения на этапе корректировки решения. Значение $\alpha = 0$ отвечает равной нулю дополнительной схемной диссипации энергии, которой соответствует второе слагаемое в левой части неравенства. Такой выбор параметра наиболее предпочтителен. С другой стороны, регулируя значение $\alpha $, можно подавлять нежелательные эффекты численного решения типа осцилляций на фронтах диссипативных ударных волн.

Если при выбранном соотношении шагов пространственно-временнóй сетки разностная схема (5.2) в совокупности с необходимыми граничными условиями задачи удовлетворяет условию устойчивости при переходе на один слой по времени, т.е. если

${{\left\| {\delta {{{\bar {U}}}^{k}}} \right\|}_{{k - 1}}} \leqslant {{\left\| {\delta {{U}^{{k - 1}}}} \right\|}_{{k - 2}}}(1 + C\Delta t),\quad C = {\text{const}} \geqslant 0,$

причем сеточная норма согласована с ассоциированной нормой вектора, то выполняется неравенство

${{\left\| {\delta {{U}^{k}}} \right\|}_{{k - 1}}} \leqslant {{\left\| {\delta {{U}^{1}}} \right\|}_{0}}exp(Ck\Delta t),$

гарантирующее устойчивость решения дискретного вариационного неравенства по начальным данным. Таким образом, применение корректировки (5.4) как этапа решения задачи при $\alpha \geqslant 0$ не приводит к развитию неустойчивости вычислений.

Можно усовершенствовать описанный метод, применив двухшаговую процедуру построения решения, первый шаг которой в точности совпадает с алгоритмом (5.2), (5.3), а второй шаг повторяет этот алгоритм после замены матриц $A({{U}^{{k - 1}}})$, ${{B}^{i}}({{U}^{{k - 1}}})$ и вектора $g({{U}^{{k - 1}}})$ матрицами и вектором:

$A\left( {\frac{{{{U}^{k}} + {{U}^{{k - 1}}}}}{2}} \right),\quad {{B}^{i}}\left( {\frac{{{{U}^{k}} + {{U}^{{k - 1}}}}}{2}} \right),\quad g\left( {{{t}^{{k - 1}}} + \frac{{\Delta t}}{2},x,\frac{{{{U}^{k}} + {{U}^{{k - 1}}}}}{2}} \right).$

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

Более сложная процедура численной реализации вариационного неравенства возникает в случае консервативной аппроксимации дифференциального оператора, которую целесообразно применять, например, при расчете обобщенных решений с сильными разрывами. Пусть теперь ${{\bar {U}}_{h}}$ – решение консервативной разностной схемы для системы квазилинейных уравнений $\mathcal{D}\left\langle U \right\rangle = g(U)$:

(5.5)

$\frac{{\varphi ({{{\bar {U}}}^{k}}) - \varphi ({{U}^{{k - 1}}})}}{{\Delta t}} = \sum\limits_{i = 1}^n \,{{\mathcal{L}}_{i}}\left\langle {{{\psi }_{i}}({{U}^{{k - 1}}})} \right\rangle + g({{U}^{{k - 1}}}).$

Тогда соответствующее дискретное неравенство может быть записано в следующей форме:

(5.6)

$(\tilde {U} - {{\hat {U}}^{k}}) \cdot (\varphi ({{U}^{k}}) - \varphi ({{\bar {U}}^{k}})) \geqslant 0,\quad {{\hat {U}}^{k}},\tilde {U} \in F.$

Принимая здесь ${{\hat {U}}^{k}} = {{U}^{k}}$, с учетом критерия выпуклости производящего потенциала $\Phi (U)$, согласно которому

$\Phi (\tilde {U}) - \Phi ({{U}^{k}}) \geqslant (\tilde {U} - {{U}^{k}}) \cdot \varphi ({{U}^{k}}),$

получаем, что искомый вектор ${{U}^{k}}$ является решением задачи минимизации при ограничениях:

(5.7)

$min\{ \Phi (\tilde {U}) - \tilde {U} \cdot \varphi ({{\bar {U}}^{k}})\,{\text{|}}\,\tilde {U} \in F\} ,$

обобщающей задачу нахождения проекции на множество $F$ при численной реализации вариационного неравенства (5.3). Справедливо также обратное утверждение: единственная точка минимума определенной здесь сильно выпуклой функции при ограничении $\tilde {U} \in F$ представляет собой решение неравенства (5.6).

Применение рассмотренной выше корректировки решения с параметром $\alpha $ в случае консервативной аппроксимации оператора нецелесообразно. Устойчивость такой процедуры и эффект снижения схемной диссипации энергии с уменьшением $\alpha $ установить, вообще говоря, не удается. Более того, численные эксперименты для модельных задач показывают, что получаемая при этом разностная схема теряет точность при наличии поверхностей сильного разрыва решения.

В качестве обобщения рассмотрим вариант корректировки решения [21], для которого

(5.8)

$\varphi ({{U}^{k}}) = \frac{2}{{1 + \alpha }}\varphi ({{\hat {U}}^{k}}) - \frac{{1 - \alpha }}{{1 + \alpha }}\varphi ({{\bar {U}}^{k}}).$

При таком выборе связи между сеточными функциями ${{\hat {U}}_{h}}$ и ${{U}_{h}}$ вариационное неравенство (5.6) приводится к виду

$(\tilde {U} - {{\hat {U}}^{k}}) \cdot (\varphi ({{\hat {U}}^{k}}) - \varphi ({{\bar {U}}^{k}})) \geqslant 0.$

Следовательно, вектор ${{\hat {U}}^{k}}$ также является решением задачи выпуклого программирования (5.7), а вектор ${{U}^{k}}$ определяется как решение системы нелинейных уравнений (5.8).

Для малых возмущений решения в данном случае справедливы неравенство и конечное соотношение

$\delta {{\hat {U}}^{k}} \cdot (\delta \varphi ({{U}^{k}}) - \delta \varphi ({{\hat {U}}^{k}})) \leqslant 0,\quad 2A({{\hat {U}}^{k}})\delta {{\hat {U}}^{k}} = (1 + \alpha )\delta \varphi ({{U}^{k}}) + (1 - \alpha )\delta \varphi ({{\bar {U}}^{k}}).$

После исключения вектора $\delta {{\hat {U}}^{k}}$ отсюда следует оценка

${\text{|}}\delta \varphi ({{U}^{k}}){\kern 1pt} {\text{|}}_{k}^{2} + \alpha {\text{|}}\delta \varphi ({{U}^{k}}) - \delta \varphi ({{\bar {U}}^{k}}){\kern 1pt} {\text{|}}_{k}^{2} \leqslant {\text{|}}\delta \varphi ({{\bar {U}}^{k}}){\kern 1pt} {\text{|}}_{k}^{2}.$

Здесь ${\text{|}}U{{{\text{|}}}_{k}} = \sqrt {U \cdot A{{{({{{\hat {U}}}^{k}})}}^{{ - 1}}}U} $ – норма, ассоциированная с матрицей, обратной к $A({{\hat {U}}^{k}})$. Таким образом, при $\alpha \geqslant 0$ корректирующая процедура обеспечивает устойчивость вычислений, если схема (5.5) устойчива по соответствующей сеточной норме.

Остановимся на вопросах численной реализации корректировки в форме (5.3) и в более общей форме, связанной с задачей (5.6). Так как при решении вариационного неравенства этап корректировки применяется многократно, на каждом слое по времени и в каждом узле пространственной сетки, то эффективность алгоритма в целом определяется скоростью выполнения этого этапа. В простейшей ситуации, когда структура ограничений, задаваемых множеством допустимых вариаций $F$, и матрица $A({{U}^{{k - 1}}})$ таковы, что проектор ${{\pi }_{{k - 1}}}(U)$ выписывается в замкнутом виде, процедура корректировки реализуется непосредственно по формулам (5.3) после вычисления ${{\bar {U}}^{k}}$ из системы уравнений (5.2). Такая ситуация встречается в большинстве приложений, касающихся изотропных упругопластических материалов и классических условий пластичности Мизеса и Треска–Сен-Венана [7].

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

(5.9)

$(\tilde {U} - \hat {U}{{{\text{'}}}^{k}}) \cdot (A({{\hat {U}}^{k}})(\hat {U}{{{\text{'}}}^{k}} - {{\hat {U}}^{k}}) + \omega (\varphi ({{\hat {U}}^{k}}) - \varphi ({{\bar {U}}^{k}}))) \geqslant 0,\quad \hat {U}{{{\text{'}}}^{k}},\tilde {U} \in F.$

Здесь $A({{\hat {U}}^{k}})$ – матрица, вычисляемая на первой итерации; $\omega > 0$ – итерационный параметр, начальное значение которого $\omega = 1$ уменьшается вдвое всякий раз при нарушении релаксационности последовательности приближений. Алгоритм нахождения нового приближения $\hat {U}{{{\text{'}}}^{k}}$ из этого неравенства через предыдущее ${{\hat {U}}^{k}}$ повторяет алгоритм (5.4) решения вариационного неравенства (5.3).

Скорость сходимости итераций (5.9) зависит, вообще говоря, от начального приближения, которое рекомендуется выбирать как ${{\hat {U}}^{k}} = {{U}^{{k - 1}}}$. Для ускорения можно через определенное число шагов производить пересчет матрицы $A({{\hat {U}}^{k}})$ с восстановлением значения $\omega = 1$. Вопросы сходимости аналогичного итерационного процесса в несколько иной интерпретации изучены в [22, с. 209].

При реализации алгоритма с параметром $\alpha = 1$ нет необходимости в последующем решении системы уравнений (5.8), так как ${{U}_{h}} = {{\hat {U}}_{h}}$. В общем случае решение можно вычислить с помощью метода Ньютона–Рафсона.

При исследовании деформирования анизотропных материалов, подчиняющихся специальным условиям пластичности, когда оператор проекции в замкнутой форме выписать не удается, достаточно эффективным может оказаться метод решения задачи выпуклого программирования на основе модифицированной функции Лагранжа (см. [23]). Если множество $F$ допускает параметризацию

$F = \{ U\,{\text{|}}\,{{f}_{j}}(U) \leqslant 0,\;j = 1,...,l\} $

через систему выпуклых функций ${{f}_{j}}(U)$, то модифицированная функция Лагранжа равна

${{\Lambda }^{k}}(U,\gamma ) = \Phi (U) - U \cdot \varphi (\mathop {\bar {U}}\nolimits^k ) + \frac{1}{{2\kappa }}\sum\limits_{j = 1}^l \,(ma{{x}^{2}}\{ 0,{{\gamma }_{j}} + \kappa {{f}_{j}}(U)\} - \gamma _{j}^{2}).$

Здесь $\kappa > 0$ – параметр метода, $\gamma $ – двойственный вектор. Эквивалентная (5.7) задача поиска седловой точки функции Лагранжа

${{\Lambda }^{k}}({{U}^{k}},{{\gamma }^{k}}) = \mathop {min}\limits_{\tilde {U}} \mathop {max}\limits_{\tilde {\gamma }} {{\Lambda }^{k}}(\tilde {U},\tilde {\gamma })$

решается при помощи итерационного процесса

${{A}^{k}}({{U}^{k}})(U{\kern 1pt} {{{\text{'}}}^{k}} - {{U}^{k}}) = - \omega \frac{{\partial {{\Lambda }^{k}}({{U}^{k}},{{\gamma }^{k}})}}{{\partial U}},\quad \gamma {{{\text{'}}}^{k}} - {{\gamma }^{k}} = \lambda \frac{{\partial {{\Lambda }^{k}}({{U}^{k}},{{\gamma }^{k}})}}{{\partial \gamma }}.$

Итерационные параметры $\omega > 0$ и $\lambda > 0$ первоначально полагаются равными единице и подбираются при помощи процедуры деления пополам из соображения релаксационности последовательности приближений для ${{U}^{k}}$ и ${{\gamma }^{k}}$.

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

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

Годунов С.К., Забродин А.В., Иванов М.Я., Крайко А.Н., Прокопов Г.П. Численное решение многомерных задач газовой динамики. М.: Наука, 1976.
Годунов С.К. Уравнения математической физики. М.: Наука, 1979.
Годунов С.К., Роменский Е.И. Элементы механики сплошных сред и законы сохранения. М.: Научная книга, 1998.
Годунов С.К., Пешков И.М. Симметрические гиперболические уравнения нелинейной теории упругости // Ж. вычисл. матем. и матем. физ. 2008. Т. 48. № 6. С. 1034–1055.
Friedrichs K.O. Symmetric hyperbolic linear differential equations // Commun. Pure Appl. Math. 1954. V. 7. № 2. P. 345–392.
Мартюшев Л.М., Селезнев В.Д. Принцип максимальности производства энтропии в физике и смежных областях. Екатеринбург: ГОУ ВПО УГТУ-УПИ, 2006.
Садовский В.М. Разрывные решения в задачах динамики упругопластических сред. М.: Наука, 1997.
Садовский В.М., Гузев М.А., Садовская О.В., Qi Ch. Моделирование пластической деформации на основе теории ортотропного континуума Коссера // Физ. мезомех. 2019. Т. 22. № 2. С. 59–66.
Садовская О.В., Садовский В.М. Математическое моделирование в задачах механики сыпучих сред. М.: Физматлит, 2008.
Садовский В.М., Садовская О.В. Анализ деформации пористой среды с учетом схлопывания пор // Прикл. механ. и техн. физ. 2016. Т. 57. № 5. С. 53–65.
Смирнов В.И. Курс высшей математики. Т. 4. Ч. 2. М.: Наука, 1981.
Овсянников Л.В. Лекции по основам газовой динамики. Москва–Ижевск: Институт компьютерных исследований, 2003.
Садовский В.М. Гиперболические вариационные неравенства в задачах динамики упругопластических тел // Прикл. матем. и механ. 1991. Т. 55. Вып. 6. С. 1041–1048.
Кукуджанов В.Н. К исследованию уравнений динамики упругопластических сред при конечных деформациях. Нелинейные волны деформаций. Таллин, 1977. Т. 2. С. 102–105.
Быковцев Г.И., Кретова Л.Д. О распространении ударных волн в упругопластических средах // Прикл. матем. и механ. 1972. Т. 36. Вып. 1. С. 106–116.
Садовский В.М. Упругопластические волны сильного разрыва в линейно упрочняющихся средах // Изв. РАН. МТТ. 1997. V. 32. № 6. С. 104–111.
Куликовский А.Г., Чугайнова А.П. Исследование разрывов в решениях уравнений упругопластической среды Прандтля–Рейсса // Ж. вычисл. матем. и матем. физ. 2016. Т. 56. № 4. С. 650–663.
Маслов В.П., Мясников В.П., Данилов В.Г. Математическое моделирование аварийного блока Чернобыльской АЭС. М.: Наука, 1988.
Садовский В.М. К теории распространения упругопластических волн в сыпучих средах // Докл. АН. 2002. Т. 386. № 4. С. 487–489.
Wilkins M.L. Calculation of elastic-plastic flow. Methods in Computational Physics: Fundamental Methods in Hydrodynamics. V. 3. P. 211–263. New York: Academic Press, 1964.
Аннин Б.Д., Садовский В.М. О численной реализации вариационного неравенства в задачах динамики упругопластических тел // Ж. вычисл. матем. и матем. физ. 1996. Т. 36. № 9. С. 177–191.
Пшеничный Б.Н., Данилин Ю.М. Численные методы в экстремальных задачах. М.: Наука, 1975.
Жильцов А.В., Намм Р.В. Метод множителей Лагранжа в задаче конечномерного выпуклого программирования // Дальневосточный матем. журнал. 2015. Т. 15. № 1. С. 53–60.

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

Инструменты

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

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

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