Доклады Российской академии наук. Математика, информатика, процессы управления, 2021, T. 501, № 1, стр. 74-78

Математическое моделирование роста материала нео-Гука

Член-корреспондент РАН П. И. Плотников^1, *

¹ Институт гидродинамики им. М.А. Лаврентьева Сибирского отделения Российской академии наук
Новосибирск, Россия

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

DOI: 10.31857/S268695432106014X

Аннотация

Дается вывод математической модели объемного роста несжимаемого материала нео-Гука. Такого рода модели используются для описании эволюции головного мозга под действием внешней нагрузки. В работе доказывается, что множество полей деформаций, соответствующих гомеостатическим состояниям, совпадает с группой Мёбиуса конформных преобразований. Дается анализ линейной краевой задачи, полученной линеаризацией уравнений нелинейной модели на гомеостатическом состоянии. Исследуется поведение решений при неограниченном росте временной переменной. Основной вывод состоит в том, что изменения в материале, вызванные временным повышением давления (гидроцефалия), являются необратимыми.

Ключевые слова: объемный рост, материал нео-Гука, уравнения Стокса, группа Мёбиуса

МАТЕМАТИЧЕСКАЯ МОДЕЛЬ. МАТЕРИАЛ НЕО-ГУКА

В настоящей статье рассматривается краевая задача для системы дифференциальных уравнений, описывающей рост биологического материала под действием приложенной нагрузки. Несжимаемый материал нео-Гука является простейшей моделью нелинейного материала, которая находит широкое применение при математическом моделировании поведения биологических субстанций и, в частности, при моделировании эволюции головного мозга. На протяжении статьи мы будем предполагать, что материал занимает ограниченную область $\Omega \subset {{\mathbb{R}}^{3}}$ с гладкой границей в пространстве расчетных (лагранжевых) координат $x = ({{x}_{1}},{{x}_{2}},{{x}_{3}})$. Состояние материала полностью характеризуется полем деформаций ${\mathbf{u}}{\text{:}}\,\,\Omega \to {{\mathbb{R}}^{3}}$ и распределением давления $p{\text{:}}\,\,\Omega \to \mathbb{R}$. Градиент поля деформаций $\nabla {\mathbf{u}}$ совпадает с матрицей Якоби отображения ${\mathbf{u}}{\text{:}}\,\,\Omega \to {{\mathbb{R}}^{3}}$. Для краткости письма мы будем использовать обозначение ${\mathbf{F}}\, = \,\nabla {\mathbf{u}}$. Мы также будем использовать обозначение ${\text{cof}}{\mathbf{A}}\,: = \,\det {\mathbf{A}}{{{\mathbf{A}}}^{{ - {\rm T}}}}$ для матрицы, сопряженной к присоединенной матрице adj A. Далее через c обозначаются различные постоянные, зависящие только от области Ω. Напомним, что материал является несжимаемым, если в Ω выполняется равенство $\det \nabla {\mathbf{u}} = 1$. Состояние гиперупругого материала полностью характеризуется плотностью упругой запасенной энергии. Для несжимаемого материала нео-Гука плотность запасенной энергии имеет вид [1]

(1)

$\Psi = \frac{1}{2}({\text{|}}\nabla {\mathbf{u}}{{{\text{|}}}^{2}} - 3).$

ОСНОВНЫЕ ПОЛОЖЕНИЯ ТЕОРИИ ОБЪЕМНОГО РОСТА НЕСЖИМАЕМОГО МАТЕРИАЛА НЕО-ГУКА

Основные факты теории объемного роста гиперупругого материала можно найти в статьях [3, 6, 7]. Наш анализ опирается на подход, развитый в работах [4] и [2]. Основным постулатом теории является гипотеза о том, что $\nabla {\mathbf{u}} = {{{\mathbf{F}}}_{e}}{{{\mathbf{F}}}_{g}}$, где матричнозначная функция ${{{\mathbf{F}}}_{g}}$ отвечает за рост материала и называется фактором роста, а ${{{\mathbf{F}}}_{e}}$ отвечает за корректирующую упругую деформацию. Плотность запасенной энергии ${{\Psi }_{g}}$ растущего материала определяется через плотность свободной энергии исходного гиперупругого материала $\Psi $ посредством равенства, см. [4],

$\begin{gathered} {{\Psi }_{g}} = {\text{det}}{{{\mathbf{F}}}_{g}}\Psi ({{{\mathbf{F}}}_{e}}) \equiv {\text{det}}{{{\mathbf{F}}}_{g}}\Psi ({\mathbf{FF}}_{g}^{{ - 1}}) \equiv \\ \equiv {\text{det}}{{{\mathbf{F}}}_{g}}\frac{1}{2}({\text{|}}{\mathbf{FF}}_{g}^{{ - 1}}{{{\text{|}}}^{2}} - 3). \\ \end{gathered} $

Тензор F_e соответствует упругим деформациям исходного несжимаемого материала нео-Гука, что влечет равенство ${\text{det}}{{{\mathbf{F}}}_{e}} = 1$. Отсюда вытекает определяющее соотношение

(2)

$\begin{gathered} {\text{det}}{\mathbf{F}} = {\text{det}}{{{\mathbf{F}}}_{g}}{\text{,}}\quad {\text{или}}\,\,{\text{эквивалентно}} \\ {\text{det}}\nabla {\mathbf{u}} = {\text{det}}{{{\mathbf{F}}}_{g}}. \\ \end{gathered} $

В большинстве случаев теория ограничивается изучением изотропного роста, когда матрица фактора роста имеет простой вид ${{{\mathbf{F}}}_{g}}(x,t)$ = w(x, t)I, $w{\text{:}}\,\,\Omega \times (0,T) \to {{\mathbb{R}}^{ + }}$. Отсутствию роста или атрофии соответствует значение фактора роста $w = 1$. Значения $w > 1$ соответствуют росту, а значения $w < 1$ соответствуют атрофии материала. Напомним, что поле деформаций является положительно ориентируемым, т.е. $\det \nabla {\mathbf{u}} > 0$ и $w > 0$. При сделанных предположениях задача о росте гиперупругого материала сводится к определению поля деформаций и фактора роста. Для этого необходимо вывести систему управляющих дифференциальных уравнений для искомых величин. Мы получим эти уравнения как следствие неравенства Клаузиуса–Дюгема. Далее мы будем предполагать, что механическая система удовлетворяет следующим принципам.

Принцип независимости движений. Определяющие реологические соотношения должны выполняться для всех достаточно гладких полей деформаций ${\mathbf{u}}$.

Второй принцип термодинамики. Для механической системы должно выполняться неравенство Клаузиуса–Дюгема. В изотермическом случае в теории объемного роста гиперупругого материала неравенство Клаузиуса–Дюгема имеет вид [4, 2 ]

(3)

$\int\limits_\Omega ^{} {\left( {\frac{{\partial {{\Psi }_{g}}({\mathbf{F}},w)}}{{\partial {\mathbf{F}}}}:{{\partial }_{t}}{\mathbf{F}} + \frac{{\partial {{\Psi }_{g}}({\mathbf{F}},w)}}{{\partial w}}{{\partial }_{t}}w\, - \,{\mathbf{T}}:{{\partial }_{t}}{\mathbf{F}}} \right)} dx\, \leqslant \,0,$

где T – тензор напряжений,

$\begin{array}{*{20}{l}} \begin{gathered} {{\Psi }_{g}}({\mathbf{F}},w) = \frac{1}{2}(w{\text{|}}{\mathbf{F}}{{{\text{|}}}^{2}} - 3{{w}^{3}}),\quad \frac{{\partial {{\Psi }_{g}}({\mathbf{F}},w)}}{{\partial {\mathbf{F}}}} = w{\mathbf{F}}, \\ \frac{{\partial {{\Psi }_{g}}({\mathbf{F}},w)}}{{\partial w}} = \frac{1}{2}{\text{|}}{\mathbf{F}}{{{\text{|}}}^{2}} - \frac{9}{2}{{w}^{2}}. \\ \end{gathered} \end{array}$

Эволюционный принцип. Фактор роста должен удовлетворять эволюционному уравнению. Это означает, что производная ${{\partial }_{t}}w$ должна полностью определяться полем деформаций u.

Из эволюционного принципа и тождества ${{\partial }_{t}}(\det \nabla {\mathbf{u}}) \equiv ({\text{cof}}\,\nabla {\mathbf{u}}):{{\partial }_{t}}\nabla {\mathbf{u}}$ вытекает, что при произвольно фиксированном поле деформаций величина $\sigma $, определенная соотношениями

$\begin{gathered} \sigma = ({\text{cof}}\nabla {\mathbf{u}}):{{\partial }_{t}}\nabla {\mathbf{u}}, \\ {{\partial }_{t}}w = \frac{1}{2}{{(\det \nabla {\mathbf{u}})}^{{\frac{{ - 2}}{3}}}}\sigma = \frac{1}{3}{{w}^{{ - 2}}}\sigma , \\ \end{gathered} $

также фиксирована. Отсюда и из принципа независимости движений следует, что в неравенстве Клаузиуса–Дюгема в качестве ${{\partial }_{t}}\nabla u$ можно выбрать произвольное векторное поле с фиксированной величиной $\sigma $, которая, в свою очередь, полностью определяется полем деформаций. Легко видеть, что требуемым условиям удовлетворяет каждое поле вида

(4)

$\begin{gathered} {{\partial }_{t}}\nabla {\mathbf{u}} = {\mathbf{\xi }} + {\mathbf{\eta }}, \\ {\mathbf{\xi }} = ({\text{cof}}\,\nabla {\mathbf{u}}{{)}^{{ - {\rm T}}}}{\text{rot}}{\mathbf{A}},\quad {\mathbf{\eta }} = ({\text{cof}}\,\nabla {\mathbf{u}}{{)}^{{ - {\rm T}}}}\nabla \varphi , \\ \end{gathered} $

где ${\mathbf{A}} \in C_{0}^{\infty }(\Omega )$ – произвольное векторное поле и φ – решение краевой задачи

$\Delta \varphi = \sigma \,\,{\text{в}}\,\,\Omega ,\quad \varphi = 0\,\,{\text{на}}\,\,\partial \Omega .$

Подстановка соотношения (4) в (3) приводит к неравенству

(5)

$\begin{gathered} \int\limits_\Omega ^{} {\left( {\frac{{\partial {{\Psi }_{g}}({\mathbf{F}},w)}}{{\partial {\mathbf{F}}}} - {\mathbf{T}}} \right):\nabla {\mathbf{\xi }}dx + } \\ + \,\int\limits_\Omega ^{} {\left\{ {\left( {\frac{{\partial {{\Psi }_{g}}({\mathbf{F}},w)}}{{\partial {\mathbf{F}}}} - {\mathbf{T}}} \right):\nabla {\mathbf{\eta }} + \frac{{\partial {{\Psi }_{g}}({\mathbf{F}},w)}}{{\partial w}}{{\partial }_{t}}w} \right\}dx \leqslant 0.} \\ \end{gathered} $

Так как поле A в представлении для ξ произвольно, то первый интеграл в левой части (5) равен нулю, что приводит к равенству

$\int\limits_\Omega ^{} {{{{({\text{cof}}\,\nabla {\mathbf{u}})}}^{{ - 1}}}{\text{div}}\left( {\frac{{\partial {{\Psi }_{g}}({\mathbf{F}},w)}}{{\partial {\mathbf{F}}}} - {\mathbf{T}}} \right) \cdot {\text{rot}}{\mathbf{A}}dx = 0} .$

Отсюда и из теоремы Вейля вытекает соотношение

$\begin{gathered} {\text{div}}\left( {\frac{{\partial {{\Psi }_{g}}({\mathbf{F}},w)}}{{\partial {\mathbf{F}}}} - {\mathbf{T}}} \right) = \\ = \,({\text{cof}}\,\nabla {\mathbf{u}})\nabla p \equiv {\text{div}}(p{\text{cof}}\,\nabla {\mathbf{u}})\,\,{\text{в}}\,\,\Omega , \\ \end{gathered} $

в котором $p{\text{:}}\,\,\Omega \to \mathbb{R}$ – функция давления. Отсюда вытекает следующее выражение для тензора напряжений:

(6)

${\mathbf{T}} = \frac{{\partial {{\Psi }_{g}}({\mathbf{F}},w)}}{{\partial {\mathbf{F}}}} + p{\text{cof}}\,\nabla {\mathbf{u}}.$

Подстановка T и η во второй интеграл в (5) приводит к неравенству

(7)

$\begin{gathered} \begin{array}{*{20}{l}} {\int\limits_\Omega ^{} {{{{\mathbf{I}}}_{1}}{{w}^{{ - 1}}}{{\partial }_{t}}wdx \leqslant 0,} } \end{array} \\ {{{\mathbf{I}}}_{1}} = \frac{1}{2}w{\text{|}}\nabla {\mathbf{u}}{{{\text{|}}}^{2}} - 3\left( {\frac{3}{2} + p} \right){{w}^{3}}. \\ \end{gathered} $

Величина ${{{\mathbf{I}}}_{1}}$ равна следу материального тензора Эшелби, ср. [4, 2 ].

УРАВНЕНИЯ ДИНАМИКИ РАСТУЩЕГО МАТЕРИАЛА НЕО-ГУКА

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

(8)

$\begin{gathered} {\text{div}}{\mathbf{T}} = {\mathbf{f}}\quad {\text{в}}\quad \Omega \times (0,T), \\ {\mathbf{Tn}} + {\mathbf{g}} = 0\quad {\text{на}}\quad \partial \Omega \times (0,T), \\ \end{gathered} $

где ${\mathbf{n}}$ – вектор внешней нормали к границе $\Omega $. Не существует общепринятых уравнений для определения фактора роста $w$. Как правило, используются различные варианты неконсервативной модели, основанной на предположении, что временная производная фактора роста является функцией фактора роста и градиента поля деформаций. Дополнительной гипотезой является предположение о том, что правая часть этого уравнения служит тензорной изотропной функцией от тензора Эшелби, см. [4, 2 ]. Для модели изотропного роста, основанной на неравенстве Клаузиуса–Дюгема, уравнение эволюции фактора роста определяется однозначно с точностью до одной скалярной функции следа тензора Эшелби ${{{\mathbf{I}}}_{1}}$. Действительно, из неравенства (7) и принципа независимости движений вытекает, что

$\begin{gathered} {{w}^{{ - 1}}}{{\partial }_{t}}w = - {\mathbf{a}}({{{\mathbf{I}}}_{1}}), \\ {\text{где}}\quad {\mathbf{a}}(s)s \geqslant 0\quad {\text{для}}\,\,{\text{всех}}\quad s \in \mathbb{R}. \\ \end{gathered} $

При выбранном a это соотношение дает уравнение для определения фактора роста. Комбинируя эти результаты с (2), (6) и (8), мы приходим к следующей системе уравнений и граничных условий для поля деформаций, давления и фактора роста:

(9)

$\begin{gathered} {\text{div}}(w\nabla {\mathbf{u}} + p{\text{cof}}\,\nabla {\mathbf{u}}) = {\mathbf{f}}, \\ {\text{det}}\nabla {\mathbf{u}} - {{w}^{3}} = 0\quad {\text{в}}\quad \Omega \times (0,T), \\ {{w}^{{ - 1}}}{{\partial }_{t}}w + {\mathbf{a}}({{{\mathbf{I}}}_{1}}) = 0\quad {\text{в}}\quad \Omega \times (0,T), \\ (w\nabla {\mathbf{u}} + p{\text{cof}}\,\nabla {\mathbf{u}}){\mathbf{n}} = - {\mathbf{g}}\quad {\text{на}}\quad \partial \Omega \times (0,T), \\ w(x,0) = {{w}_{0}}(x)\quad {\text{в}}\quad \Omega . \\ \end{gathered} $

ГОМЕОСТАТИЧЕСКИЕ СОСТОЯНИЯ

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

${\text{|}}\nabla {\mathbf{u}}{{{\text{|}}}^{2}} = 3{{w}^{2}},\quad \det \nabla {\mathbf{u}} = {{w}^{3}},\quad {{{\mathbf{I}}}_{1}} = 0.$

Первые два равенства равносильны тому, что отображение ${\mathbf{u}}{\text{:}}\,\,\Omega \to {{\mathbb{R}}^{3}}$ является конформным. Кроме того, из приведенных соотношений следует, что

(10)

$w = (\det {\mathbf{u}}{{)}^{{1/3}}}{{ = 3}^{{ - 1/2}}}{\text{|}}\nabla {\mathbf{u}}{\text{|}},\quad p = - 1.$

Тем самым мы приходим к следующему утверждению.

Лемма 1. Тройка $({\mathbf{u}},p,w)$ определяет гомеостатическое состояние тогда и только тогда, когда отображение ${\mathbf{u}}{\text{:}}\,\,\Omega \to {{\mathbb{R}}^{3}}$ является конформным, фактор роста и давление определяются равенствами (10). Каждое гомеостатическое состояние является стационарным решением краевой задачи (9).

Согласно теореме Лиувилля, каждое конформное отображение u в ${{\mathbb{R}}^{3}}$ является суперпозицией преобразований сдвига, растяжения, вращения и инверсии. Физический смысл имеют только такие преобразования инверсии, у которых полюс инверсии находится вне области $\Omega $. Совокупность всех конформных отображений в ${{\mathbb{R}}^{3}}$ образует группу Мёбиуса. Группа Мёбиуса в ${{\mathbb{R}}^{3}}$ является 10-мерным многообразием. Касательное пространство к этому многообразию в единице (инфинитезимальное пространство группы) состоит из всех отображений вида

(11)

${\mathbf{\zeta }}(x) = {\mathbf{b}} + {\mathbf{S}}x + \lambda x + (A \cdot x)x - \frac{1}{2}{\text{|}}x{{{\text{|}}}^{2}}A,$

где ${\mathbf{b}},A \in {{\mathbb{R}}^{3}}$ – произвольные векторы, $\lambda $ – произвольная постоянная и S – произвольная антисимметрическая матрица. Инфинитезимальное пространство совпадает с пространством решений обобщенных уравнений Киллинга группы конформных преобразований ${{\mathbb{R}}^{3}}$,

$\nabla {\mathbf{\zeta }} + \nabla {{{\mathbf{\zeta }}}^{{\rm T}}} - \frac{2}{3}{\text{div}}{\mathbf{\zeta I}} = 0.$

Определение 1. Через $T\mathcal{M}$ обозначим линейное пространство всех отображений ζ, допускающих представление (11). Через $T{{\mathcal{M}}_{{{\text{div}}}}}$ обозначим линейное пространство дивергенций всех векторных полей, принадлежащих $T\mathcal{M}$.

ЛИНЕАРИЗОВАННЫЕ УРАВНЕНИЯ

Исследование системы (9) целесообразно начать с анализа линеаризованных уравнений. Особый интерес представляет исследование линеаризации на гомеостатических состояниях. Из групповой инвариантности множества гомеостатических состояний следует, что можно ограничиться изучением линеаризации на тривиальном равновесном состоянии ${\mathbf{u}} = {\mathbf{id}}$, $p = - 1$, $w = 1$. Обозначим через ${\mathbf{v}}$, π, $\varphi $ отклонения поля деформаций, давления и фактора роста от равновесного состояния, вызванные действием распределенной силы f и поверхностной силы g. Разложим поле деформаций ${\mathbf{v}}$ и давление π на две составляющие

${\mathbf{v}} = {{{\mathbf{v}}}_{f}} + {\mathbf{w}},\quad {\mathbf{\pi }} = {{{\mathbf{\pi }}}_{f}} + {\mathbf{\mu }},$

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

(12)

$\begin{gathered} {\text{div}}(\nabla {{{\mathbf{v}}}_{f}} + \nabla {\mathbf{v}}_{f}^{ \top } + {{{\mathbf{\pi }}}_{f}}{\mathbf{I}}) = {\mathbf{f}},\quad {\text{div}}{{{\mathbf{v}}}_{f}} = 0\quad {\text{в}}\quad \Omega , \\ (\nabla {{{\mathbf{v}}}_{f}} + \nabla {\mathbf{v}}_{f}^{ \top } + {{{\mathbf{\pi }}}_{f}}{\mathbf{I}}) \cdot {\mathbf{n}} + {\mathbf{g}} = 0\quad {\text{на}}\quad \partial \Omega , \\ \end{gathered} $

(13)

$\begin{gathered} {\text{div}}\left( {\nabla {\mathbf{w}} + \nabla {{{\mathbf{w}}}^{ \top }} - \frac{2}{3}{\text{div}}\,{\mathbf{wI}} + {\mathbf{\mu I}}} \right) = 0\quad {\text{в}}\quad \Omega , \\ {\text{div}}{\mathbf{w}} = 3\varphi \quad {\text{в}}\quad \Omega , \\ \left( {\nabla {\mathbf{w}} + \nabla {{{\mathbf{w}}}^{{\rm T}}} - \frac{2}{3}{\text{div}}\,{\mathbf{wI}} + {\mathbf{\mu I}}} \right){\mathbf{n}} = 0\quad {\text{на}}\quad \partial \Omega , \\ \end{gathered} $

(14)

$\begin{gathered} {{\partial }_{t}}\varphi = \alpha {\mathbf{\mu }} + \alpha {{{\mathbf{\pi }}}_{f}}\quad {\text{в}}\quad \Omega \times (0,T), \\ \varphi (x,0) = {{\varphi }_{0}}(x)\quad {\text{в}}\quad \Omega , \\ \end{gathered} $

где $\alpha = 3{\mathbf{a}}{\kern 1pt} '(0) > 0$, ${{\varphi }_{0}}$ – заданное начальное распределение фактора роста. Нашей целью является исследование корректности краевой задачи (12)–(14). Для этого рассмотрим следующую конструкцию. Разложим инфинитезимальное пространство $T{{\mathcal{M}}_{{{\text{div}}}}}$ в прямую сумму подпространств

$\begin{array}{*{20}{l}} \begin{gathered} {{\mathcal{E}}_{{{\text{el}}}}} = \{ {\mathbf{c}} + {\mathbf{S}}x{\text{:}}\,\,{\mathbf{c}} \in {{\mathbb{R}}^{3}},{\mathbf{S}} = - {{{\mathbf{S}}}^{{\rm T}}}\} , \\ {{\mathcal{E}}_{{{\text{gr}}}}} = \left\{ {\lambda x + ({\mathbf{A}} \cdot x)x - \frac{{{\text{|}}x{{{\text{|}}}^{2}}}}{2}{\mathbf{A}},\lambda \in \mathbb{R},{\mathbf{A}} \in {{\mathbb{R}}^{3}}} \right\}. \\ \end{gathered} \end{array}$

Обозначим через $\mathcal{V}$ и $\mathcal{W}$ подпространства ${{W}^{{1,2}}}(\Omega )$, а через $\mathcal{H}$ – замкнутое подпространство ${{L}^{2}}(\Omega )$, определенные следующими равенствами:

$\begin{gathered} \mathcal{V}\, = \,\left\{ {{\mathbf{v}}\, \in \,{{W}^{{1,2}}}(\Omega ){\text{:}}\,\,\int\limits_\Omega ^{} {{\mathbf{v}} \cdot \zeta dx = 0\,\,{\text{для}}\,\,{\text{всех}}\,\,\zeta \in {{\mathcal{E}}_{{{\text{el}}}}}} } \right\}, \\ \mathcal{W}\, = \,\left\{ {{\mathbf{w}}\, \in \,{{W}^{{1,2}}}(\Omega ){\text{:}}\,\,\int\limits_\Omega ^{} {{\mathbf{w}} \cdot \zeta dx = 0\,\,{\text{для}}\,\,{\text{всех}}\,\,\zeta \in T\mathcal{M}} } \right\}, \\ \mathcal{H}\, = \,\left\{ {\varphi \, \in \,{{L}^{2}}(\Omega ){\text{:}}\,\,\int\limits_\Omega ^{} {\varphi \cdot \zeta dx = 0\,\,{\text{для}}\,\,{\text{всех}}\,\,\zeta \in T{{\mathcal{M}}_{{{\text{div}}}}}} } \right\}. \\ \end{gathered} $

Через $\Pi {\text{:}}\,\,{{L}^{2}}(\Omega ) \to \mathcal{H}$ мы будем обозначать оператор ортогонального проектирования пространства ${{L}^{2}}(\Omega )$ на подпространство $\mathcal{H}$.

Подпространство ${{\mathcal{E}}_{{{\text{el}}}}}$ состоит из полей деформаций, которые порождены твердотельными движениями упругого материала и не влияют на процесс роста. Поэтому целесообразно факторизовать пространства решений задач (12) и (13) по подпространству ${{\mathcal{E}}_{{el}}}$. Далее мы будем предполагать, что ${{{\mathbf{v}}}_{f}} \in \mathcal{V}$ и ${\mathbf{w}} \in {{\mathcal{E}}_{{{\text{gr}}}}} + \mathcal{W}.$ Следующие два утверждения, гарантируют корректность задач (12) и (13).

Предложение 1. Пусть пара $({\mathbf{f}},{\mathbf{g}}) \in {{L}^{2}}(\Omega )$ × × L²$(\partial \Omega )$, уравновешена, т.е. удовлетворяет условиям

$\int\limits_\Omega ^{} {{\mathbf{f}} \cdot \zeta dx + \int\limits_{\partial \Omega }^{} {{\mathbf{g}} \cdot \zeta ds = 0\,\,для\,\,всех\,\,\zeta \in {{\mathcal{E}}_{{{\text{el}}}}}.} } $

Тогда задача (12) имеет единственное слабое решение ${{{\mathbf{v}}}_{f}} \in \mathcal{V}$, ${{{\mathbf{\pi }}}_{f}} \in {{L}^{2}}(\Omega )$, допускающее оценку

${\text{||}}{{{\mathbf{v}}}_{f}}{\text{|}}{{{\text{|}}}_{{{{W}^{{1,2}}}(\Omega )}}} + \,{\text{||}}{{{\mathbf{\pi }}}_{f}}{\text{|}}{{{\text{|}}}_{{{{L}^{2}}(\Omega )}}} \leqslant c\left( {{\text{||}}{\mathbf{f}}{\text{|}}{{{\text{|}}}_{{{{L}^{2}}(\Omega )}}} + \,{\text{||}}{\mathbf{g}}{\text{|}}{{{\text{|}}}_{{{{L}^{2}}(\partial \Omega )}}}} \right).$

Предложение 2. Для любого $\varphi \in {{L}^{2}}(\Omega )$ задача (13) имеет единственное слабое решение ${\mathbf{w}} \in {{\mathcal{E}}_{{{\text{gr}}}}} + \mathcal{W}$, ${\mathbf{\mu }} \in \mathcal{H}$, допускающее оценку

${\text{||}}{\mathbf{w}}{\text{|}}{{{\text{|}}}_{{{{W}^{{1/2}}}(\Omega )}}} + \,{\text{||}}{\mathbf{\mu }}{\text{|}}{{{\text{|}}}_{{{{L}^{2}}(\Omega )}}} \leqslant c{\text{||}}\varphi {\text{|}}{{{\text{|}}}_{{{{L}^{2}}(\Omega )}}}.$

Функция μ обращается в нуль тогда и только тогда, когда $\varphi \in T{{\mathcal{M}}_{{{\text{div}}}}}$. Более того, справедлива оценка

$\int\limits_\Omega ^{} {{\mathbf{\mu }}\varphi dx \leqslant - {{c}^{{ - 1}}}{\text{||}}\Pi \varphi {\text{||}}_{{{{L}^{2}}(\Omega )}}^{2}.} $

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

$\begin{array}{*{20}{l}} {{\mathbf{\Phi }}({\mathbf{w}}*,{\mathbf{\mu }}*) = \mathop {\max }\limits_{{\mathbf{\mu }} \in \mathcal{H}} \mathop {\min }\limits_{{\mathbf{w}} \in \mathcal{W}} {\mathbf{\Phi }}({\mathbf{w}},{\mathbf{\mu }})} \end{array}$

для функционала

$\begin{gathered} \begin{array}{*{20}{l}} {{\mathbf{\Phi }}({\mathbf{w}},{\mathbf{\mu }}) = \frac{1}{4}\int\limits_\Omega ^{} {{{{\left| {\nabla {\mathbf{w}} + \nabla {{{\mathbf{w}}}^{{\rm T}}} - \frac{2}{3}{\text{div}}\,{\mathbf{wI}}} \right|}}^{2}}dx + } } \end{array} \\ + \,\int\limits_\Omega ^{} {{\mathbf{\mu }}({\text{div}}\,{\mathbf{w}} - 3\varphi )dx.} \\ \end{gathered} $

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

Лемма 2. Пусть $({\mathbf{w}},{\mathbf{\mu }})$ – решение задачи (13) с произвольной функцией $\varphi \in {{L}^{2}}(\Omega )$. Тогда отображение ${\mathbf{B}}{\text{:}}\,\,\varphi \to {\mathbf{\mu }}$ определяет ограниченный симметричный оператор в ${{L}^{2}}(\Omega )$. Ортогональные подпространства $T{{\mathcal{M}}_{{{\text{div}}}}}$ и $\mathcal{H}$ являются инвариантными подпространствами B. Оператор B обращается в нуль на $T{{\mathcal{M}}_{{{\text{div}}}}}$ и строго отрицателен на $\mathcal{H}$, т.е.

$\begin{gathered} {\mathbf{B}}({\mathbf{I}} - \Pi ) = 0, \\ {{\langle {\mathbf{B}}\varphi ,\varphi \rangle }_{{{{L}^{2}}(\Omega )}}} \leqslant - {{c}^{{ - 1}}}{\text{||}}\Pi \varphi {\text{||}}_{{{{L}^{2}}(\Omega )}}^{2}\,\,для\,\,всех\,\,\varphi \in {{L}^{2}}(\Omega ). \\ \end{gathered} $

Далее, из предложения 1 вытекает, что π_f полностью определяется силами f, g. Поэтому уравнение (14) может быть записано в форме операторного уравнения для фактора роста $\varphi $ :

(15)

$\begin{gathered} {{\partial }_{t}}\varphi = \alpha {\mathbf{B}}\varphi + \alpha {{{\mathbf{\pi }}}_{f}}\,\,{\text{на}}\,\,{\text{интервале}}\,\,(0,\infty ), \\ \varphi (0) = {{\varphi }_{0}}. \\ \end{gathered} $

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

Теорема 1. Пусть заданные силы ${\mathbf{f}} \in {{L}^{1}}(0,\infty $; L²(Ω)), ${\mathbf{g}} \in {{L}^{1}}(0,\infty ;{{L}^{2}}(\partial \Omega ))$ уравновешены в каждый момент времени. Тогда задача (15) имеет единственное решение $\varphi \in C(0,\infty ;{{L}^{2}}(\Omega ))$, которое при неограниченном возрастании временной переменной стремится к предельному значению

$\begin{gathered} \varphi (t) \to ({\mathbf{I}} - \Pi ){{\varphi }_{0}} + \alpha \int\limits_0^\infty {({\mathbf{I}} - \Pi ){{{\mathbf{\pi }}}_{f}}(s)ds} \\ в\quad {{L}^{2}}(\Omega )\quad при\quad t \to \infty . \\ \end{gathered} $

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

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

Ciarlet P. Mathematical elasticity 1: Three-dimensional Elasticity. Basel: Elsevier Science Publishers, 1988.
Ciarletta P., Ambrosi D., Maugin G.A. Mass transport in morphogenetic processes: a second gradient theory for volumetric growth and material remodeling // J. Mech. Phys. Solids. 2012. V. 60. P. 432–450.
Cowin S.C. Tissue growth and remodeling // Annu. Rev. Biomed. Eng. 2004. V. 6. P. 77–107.
Epstein M., Maugin G.A. Thermomechanics of volumetric growth in uniform bodies//Int. J. Plasticity. 2000. V. 16. P. 951–978.
Решетняк Ю.Г. Оценки для некоторых дифференциальных операторов с конечномерным ядром// Сиб. матем. журн. 1970. Т. 11. № 2. С. 414–428.
Rodriguez E., Hoger A., McCulloch A. Stress-dependent finite growth law in soft elastic tissue // J. Biomech. 1994. V. 27. P. 455–467.
Skalak R., Dasgupta G., Moss M., Otten E., Dullemeijer P., Vilmann H. Analytical description of growth.// J. Theor. Biol. 1982. V. 94. P. 555–577.

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

Инструменты

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

Доклады Российской академии наук. Математика, информатика, процессы управления

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