Evolution of the field of distributed defects in a crystal during contact interaction with a system of rigid stamps

Full Text

Введение

В настоящее время в микросистемной технике широко применяются микроэлектромеханические устройства (MEMS). Характерный размер таких устройств – от нескольких десятков нанометров до микрона [1; 2]. В отличие от микроэлектронных компонент элементы MEMS подвержены воздействию не только электромагнитных полей, но также действию объемных и поверхностных сил, возникающих в областях механических контактов подвижных частей [3]. Вызванное этими усилиями напряженно-деформированное состояние микроэлектромеханической системы оказывает существенное влияние на ее эксплуатационные свойства и должно учитываться при расчетах и проектировании. Поскольку элементы таких устройств имеют характерные размеры порядка микрона и менее, моделирование их напряженного состояния обладает масштабной спецификой: существенное влияние на него оказывают силы поверхностного натяжения и внутренние источники напряжений, возникающие из-за присутствия дефектов в кристаллической структуре образующего эти элементы материала. В качестве последнего, как правило, используют кристаллический или аморфный кремний. Распределение дефектов в кристалле зависит от напряженного состояния, и его изменение приводит к перераспределению дефектов в объеме кристалла — дефекты мигрируют. Подобная миграция с макроскопической точки зрения может быть представлена как эволюция неупругих (пластических) деформаций. В этой связи актуальным как в теоретическом, так и в прикладном плане, является математическое моделирование эволюции напряженного состояния и связанной с ней структуры дефектов в микромасштабных кристаллических телах при их контактном взаимодействии. В настоящей работе развиваются общие способы постановок подобных контактных задач и рассматриваются модельные примеры о квазистатическом взаимодействии упругого кристалла с системой жестких штампов. При этом в кристалле задается некоторое начальное распределение дефектов, изменение которого в процессе контактного взаимодействия определяется из решения эволюционной задачи.

1. Ход исследования

1°. Континуальная теория дефектов давно является классическим разделом нелинейной механики деформируемого твердого тела, однако её основные положения и методы продолжают интенсивно развиваться и определять перспективные на сегодняшний день общетеоретические и прикладные задачи. Близким в смысле объекта исследования, но существенно отличающимся по методам и моделям является раздел физики твердого тела, изучающий дефекты как нарушения регулярной структуры кристалла. Основное различие состоит в способе математической формализации влияния совокупности дефектов на макроскопические свойства содержащего их тела. Континуальная теория, что соответствует ее названию, представляет всю совокупность дефектов как непрерывное поле локальных деформаций, не удовлетворяющих условиям совместности и вызывающих дополнительную деформацию и соответствующее напряженно-деформированное состояние, подобно тому, как неравномерный нагрев вызывает термоупругие напряжения. Физические теории дефектов, в частности теория дислокаций, рассматривают дискретные структуры кристаллов с нарушенным порядком и, образно говоря, следят за отдельными нарушителями и изменением энергетического состояния кристалла при их миграциях между узлами кристаллической решетки. Несмотря на впечатляющие вычислительные возможности, предоставляемые современными компьютерами, прямое моделирование дискретных структур с дефектами, например, методами молекулярной динамики можно выполнить за разумное время только в очень небольших объемах и на интервалах времени, не превышающих нескольких миллисекунд. В этой связи синтез двух указанных подходов, при котором моделирование влияния большого числа дефектов на напряженно деформированное тело осуществляется в рамках континуального подхода, в то время как уравнения эволюции полей дефектов определяются физическими моделями взаимодействия дефекта с соседними атомами основного материала, представляет актуальную задачу современной механики континуума. Настоящая работа представляет некоторые результаты, полученные в этом направлении.

Исторически континуальная и атомистическая теории дефектов возникли независимо друг от друга, но почти одновременно. Начала континуальной теории были заложены в работах Фойгта, Вайнгартена, Вольтерра более 100 лет назад [4; 5]. Основным объектом в этих исследованиях были уравнения совместности Сен-Венана и возможные способы модификации тела, которые приводили к тому, что уравнения Сен-Венана переставали выполняться (для чего изменялся индекс связности тела за счет разрезов, склеек и пр.). В некотором смысле эти исследования были чисто теоретическими: они позволили математически строго сформулировать понятие несовместной деформации и ее внутренних источников — дислокаций Вольтерра. Следует отметить, что это понятие, несмотря на название, не имело отношения к физике атомистических структур, утративших по тем или иным причинам регулярность своего строения. Но как раз этот вопрос оказался в фокусе внимания в начале XX века, когда теория атомистического строения вещества получила признание (Эйнштейн-броуновское движение, Резерфорд планетарная модель атома, Бор – старая квантовая теория) и стала развиваться семимильными шагами. С аппаратом квантовой механики, казалось, легко можно было бы реализовать грандиозный замысел Коши — получить все уравнения механики континуума из анализа атомистического строения вещества [6; 7]. Однако математическое моделирование макроскопических свойств привело к столь разительному отличию от эксперимента (Френкель [8]), что стала ясна необходимость поиска какой-то совершенно новой идеи. Вскоре она была высказана Поляни [9], Орованом [10] и Тейлором [11; 12]: причина неудовлетворительного соотношения теоретических предсказаний и эксперимента кроется том, что кристаллическая структура материала предполагалась идеальной, в то время как реальные структуры насыщены дефектами. Они-то и являются внутренними источниками напряжений, подобно разрезам и склейкам в теории Вольтерра, а это ” собственное” напряженное состояние кардинально меняет физики-механические свойства макроскопического тела, которые в эксперименте определяются как отклик на внешние поля. Действительно, математическая теория идеального кристалла предполагала, что эти поля являются единственной причиной напряжений, а в эксперименте они накладывались на уже присутствующее в кристалле, самоуравновешенное, но отнюдь не малое напряженное состояние.

Построение континуальной теории дислокаций было начато в работах [11; 12], и достаточно быстро было обнаружено большое сходство между уравнениями, возникающими при описании несовместных деформаций, вызываемых дислокациями, и уравнениями общей теории относительности (Экхарт [13], Бильби [14] и Кондо [15]). Это привело к идее использовать методы неевклидовой геометрии для описания дополнительных полевых величин, характеризующих распределенные дефекты как источники внутренних напряжений (Кренер [16]). Образно говоря, тело, деформации которого несовместны, не может быть целиком преобразовано в ненапряженное состояние посредством какой-либо непрерывной деформации в объемлющем физическом пространстве, однако можно формально перенести его в пространство большего числа измерений, в котором оно найдет свою ненапряженную форму. Эта форма будет существовать на некоторой гиперповерхности расширенного пространства, геометрические свойства которой, такие как кривизна или кручение, не будут отвечать условиям евклидового пространства, и именно эти дополнительные геометрические параметры будут служить математическим описанием отличительной особенности тела, вызванной наличием в нем дефектов. Конечно, эта иллюстрация очень груба, и для более детального обсуждения методов неевклидовой геометрии в приложении к обсуждаемой задаче потребуется более детальное описание элементов теории. Это будет сделано дальше.

Существенное развитие математическая теория распределенных дефектов получила в работах Нолла [17], Вана [18] и Де Вита [19]. В них было сформулировано математически строго понятие материального многообразия как структуры, полностью описывающей совокупность материальных точек, составляющих тело и служащей основой для построения различных геометрий, евклидовых и неевклидовых, посредством задания метрики и аффинной связности на них. Кроме того, были выделены источники внутренних напряжений особого вида — дисклинации, которые возникали при определенных условиях в дополнение к распределенным дислокациям, но оказывали влияние на самонапряженное состояние образом, отличным от последних. Замечательное достижение математической теории дефектов состояло в том, что дефектам различных типов нашлись принципиально различные геометрические образы: дислокациям соответствует кручение материальной связности, дисклинациям — кривизна, объемным дефектам — неметричность.

2°. В работе используются методы анализа на гладких многообразиях, с которыми можно ознакомиться например, в руководствах Спивака [20] и Ли [21]. Приложение этих методов к механике континуума подробно изложено в монографиях Марсдена & Хьюгса [22], Ракотомананы [23], Френкеля [25], Эпстейна [24], Штайнмана [26], Лычева & Койфмана [27].

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

Итак, тело $\mathfrak{B}$ может быть математически формализовано посредством следующей структуры, задаваемой упорядоченной четверкой

$\mathfrak{B}:=\left(\mathcal{B}\mathrm{,}\mathcal{F}\mathrm{,}\mathcal{A}\mathrm{,}\mu \right)\mathrm{.}$

Первый элемент четверки, $\mathcal{B}$, задает некоторое непустое множество. Его элементы уникальным образом идентифицируют материальные точки, составляющие тело. Второй элемент, $\mathcal{F}$, определяет топологию на этом множестве, превращая его в топологическое пространство. Третий элемент, $\mathcal{A}$, представляет собой максимальный атлас, элементы которого — карты — позволяют установить соответствие между элементами множества $\mathcal{B}$ и упорядоченными n-ками действительных чисел, т. е. элементами ${ℝ}^n$. Такие соответствия в областях картрирования $U_{\alpha }$ задаются картрирующими отображеними - гомеоморфизмами ${\phi }_{\alpha }:\mathcal{B}\supset U_{\alpha }\to {ℝ}^n$, где индекс $\alpha$ пробегает множество $\mathcal{I}$, индексирующее все карты атласа. При этом требуется выполнение двух свойств: полноты и $C^k$ согласованности.

$\underset{\alpha \in \mathcal{I}}{\cup }{U}_{\alpha }\mathfrak{B}\mathrm{,} \forall \alpha \mathrm{,}\beta \in \mathcal{I} {\phi }_{\alpha }\circ {\overline{){\phi }_{\beta }^{-1}}}_{{\phi }_{\beta }(U_{\alpha }\cap U_{\beta })}\in C^k\left({ℝ}^n\right)$

Последний элемент определяет меру на $\mathfrak{B}$, с помощью которой можно выполнять интегрирование на $\mathfrak{B}$ и его подмногообразиях. Пока что все перечисленные элементы никак не конкретизированы и, по существу, представляют стандартные конструкции дифференциальной геометрии [21]. Это, с одной стороны, подчеркивает тот факт, что тело не является непосредственно наблюдаемым, а представляет себя посредством форм-образов вложения в физическое пространство, т. е. более конкретное определение может быть дано только после формализации соответствующих структур для физического пространства. С другой стороны, это обстоятельство предоставляет определенную свободу для определения на теле структур более общих, чем те, что можно получить из физического пространства. Эта возможность существенно используется в настоящей работе для построения неевклидовой отсчетной формы [28; 29].

3°. Тело как совокупность материальных точек (либо некоторых абстрактных меток этих точек, что эквивалентно приведенному выше определению, поскольку здесь материальные точки не несут дополнительных физических атрибутов, таких как масса, заряд и пр.) не является непосредственно наблюдаемым и в некотором смысле находится в абстрактном платоновом мире идей. Наблюдаемыми будут лишь формы тела-образы его вложения в физическое пространство, которые можно математически формализовать в том же ключе, как подмногообразия некоторого объемлющего гладкого многообразия, формально характеризующего физическое пространство. Но, в отличие от тела, на которое при его формализации не накладывалось никаких ограничений на выбор геометрии (связности и метрики), формы оказываются не такими свободными. Поскольку они по своему физическому смыслу должны представлять некоторые регионы физического пространства, то геометрия на них должна быть индуцирована последним. Это означает, что многообразия, представляющие формы, наделяются вполне определенной связностью и метрикой, индуцированной (образно говоря, ”скопированной”) объемлющим физическим пространством. Таким образом, для полного задания форм требуется представить физическое пространство тоже как многообразие, но наделенное специальной метрикой и связностью. Оставаясь в рамках классической физики, определим эти свойства с помощью классической вспомогательной структуры, которая, однако, оказывается весьма полезной сама по себе в дальнейших построениях.

В рамках классического (нерелятивистского) подхода физическое пространство $\mathfrak{P}$ может быть математически формализовано как аффинное пространство, определяемое структурой

$\mathfrak{P}:=\left(\mathcal{P}\mathrm{,}\mathcal{V}\mathrm{,}\text{v}ec\mathrm{,}o\right)\mathrm{.}$ (1)

Первый элемент этой структуры, $\mathcal{P}$, представляет собой непустое множество. Элементы этого множества с физической точки зрения указывают на места, занимаемые материальными точками в составе формы. Вторым элементом является векторное трансляционное пространство $\mathcal{V}$, которое само является структурой над полем действительных чисел $ℝ$:

$\mathcal{V}=\left(\mathfrak{V}\mathrm{,}ℝ,\oplus ,\odot \right)\mathrm{,}$

где $\mathfrak{V}$ — множество, элементами которого определяют трансляционные векторы, $\oplus$ — операция внутренней композиции $\mathfrak{V}\times \mathfrak{V}\to \mathfrak{V}$, $\odot$ — операция внешней композиции $ℝ\times \mathfrak{V}\to \mathfrak{V}$. Векторное пространство обладает характеристикой — размерностью. В работе она полагается равной трем, и эта величина интерпретируется как размерность аффинного пространства $\mathfrak{P}$. Третий элемент структуры $\text{v}ec$ задает операцию, согласующую элементы множества $\mathcal{P}$ и трансляционные векторы из $\mathcal{V}$, она определяется как отображение $\mathcal{P}\times \mathcal{P}\to \mathcal{V}$, удовлетворяющее условиям

 $\forall O\in \mathcal{P} V_O:\mathcal{P}∍A\mapsto vec\left(O\mathrm{,}A\right)\in \mathcal{V}$ биективна,

$\forall a\mathrm{,}b\mathrm{,}c\in \mathcal{P} vec\left(a\mathrm{,}b\right)+vec\left(b\mathrm{,}c\right)+vec\left(c\mathrm{,}a\right)=0\in \mathcal{V}\mathrm{.}$

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

Выбор некоторой точки $O\in \mathcal{P}$ и некоторого базиса в $\mathcal{V}*$, скажем ${\left({\mathrm{e}}^{\mathrm{\alpha }}\right)}_{\alpha =1}^3$, позволяет установить взаимно однозначное соответствие между точками $\mathcal{P}$ и упорядоченными n-ками чисел. Это соответствие задается отображением, которое далее будем называть арифметизацией пространства $\mathcal{P}$:

$C_{O;{\left({\mathrm{e}}^{\mathrm{\alpha }}\right)}_{\alpha =1}^3}:={\left({\mathrm{e}}^{\mathrm{\alpha }}\mathrm{\circ }\mathrm{vec}\mathrm{(}\mathrm{O}\mathrm{,}\mathrm{\cdot }\mathrm{)}\right)}_{\alpha =1}^3\mathrm{.}$  (2)

Алгебраическая структура физического пространства, построенная к настоящему моменту изложения, достаточно бесформенна — в ней не определены метрические свойства, такие как длины, углы, не задано свойство ортогональности. Все это определяется заданием метрики, наделяющей аффинное пространство евклидовой структурой. Метрика может быть задана как билинейное отображение, ставящее в соответствие каждой паре векторов $(u,v)$ число

$g(u,v)=\alpha \in ℝ\mathrm{,}$

удовлетворяющее условиям

$\forall u\mathrm{,}v\in \mathfrak{P} g\left(u\mathrm{,}v\right)=g\left(v\mathrm{,}u\right), \forall u\in \mathfrak{P} g\left(u\mathrm{,}u\right)\ge \mathrm{0,} \left(g\right(u\mathrm{,}u)=0)\iff (u=0)$

В координатной форме это отображение может быть представлено тензором типа (0,2):

$g=g_{\alpha \beta }{e}^{\alpha }\otimes e^{\beta }\mathrm{,}$

где $g_{\alpha \beta }$ – элементы симметричной положительно определенной матрицы, $e^{\mathrm{\alpha }}$ – элементы дуального базиса в $\mathfrak{V}*$. Выбор конкретных значений метрической матрицы $g_{\alpha \beta }$ осуществляется относительно некоторого выбранного базиса, выбор которого апеллирует к некоторой процедуре экспериментальной идентификации, которую, в свою очередь, можно осуществить только после полного определения метрической матрицы. Эта ситуация похожа на классическую проблему определения шахматной фигуры, скажем коня: его абстрактное определение основано на правилах игры, которые, в свою очередь, требуют для полного описания предварительное определение всех фигур, в том числе и коня. Для выхода из порочного круга следует сослаться на известную теорему алгебры, согласно которой для любой положительно определенной симметричной матрицы существует базис, в котором она приводится к единичной. В этой связи можно полагать, что выбор метрики и базиса, в котором она представлена единичной матрицей, осуществляется в рамках единого акта идентификации, а выбранный базис приобретает особый статус декартового базиса $c^{\mathbf{\alpha }}$. В нем

$g={\delta }_{\alpha \beta }{c}^{\alpha }\otimes c^{\beta }\mathrm{.}$

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

$e_{\mathrm{\alpha }}={\mathrm{\Omega }}_{\mathrm{\alpha }}^{\mathrm{\beta }}{\mathrm{c}}_{\mathrm{\beta }}\mathrm{,} {\mathrm{e}}^{\mathrm{\alpha }}={\mathrm{\mho }}_{\mathrm{\beta }}^{\mathrm{\alpha }}{\mathrm{c}}^{\mathrm{\beta }}\mathrm{,} \mathrm{g}={\mathrm{g}}_{\mathrm{\alpha \beta }}{\mathrm{e}}^{\mathrm{\alpha }}\otimes {\mathrm{e}}^{\mathrm{\beta }}\mathrm{,} {\mathrm{\Omega }}_{\mathrm{\alpha }}^{\mathrm{\beta }}{\mathrm{\mho }}_{\mathrm{\gamma }}^{\mathrm{\alpha }}={\mathrm{\delta }}_{\mathrm{\gamma }}^{\mathrm{\beta }}\mathrm{,} {\mathrm{g}}_{\mathrm{\alpha \beta }}={\mathrm{\delta }}_{\mathrm{\mu \nu }}{\mathrm{\Omega }}_{\mathrm{\alpha }}^{\mathrm{\mu }}{\mathrm{\Omega }}_{\mathrm{\beta }}^{\mathrm{\nu }}\mathrm{.}$

Заметим, что поскольку при определении $\mathfrak{P}$ была выбрана одна их возможных ориентаций, то $\det {\Omega }_{\alpha }^{\beta }>0$.

Метрический тензор определяет согласованные с ним меры длины, площади объема посредством внешних форм

$e=c^1\wedge c^2\wedge c^3=J{e}^1\wedge e^2\wedge e^3\mathrm{,} s_1=c^2\wedge {\mathrm{c}}^3={\mathrm{J}}_1{\mathrm{e}}^2\wedge {\mathrm{e}}^3\mathrm{,} {\mathrm{s}}_2={\mathrm{c}}^3\wedge {\mathrm{c}}^1={\mathrm{J}}_2{\mathrm{e}}^3\wedge {\mathrm{e}}^1\mathrm{,} {\mathrm{s}}_3={\mathrm{c}}^1\wedge {\mathrm{c}}^2={\mathrm{J}}_3{\mathrm{e}}^1\wedge {\mathrm{e}}^2\mathrm{,}$ (3)

где $J=det\left[{\Omega }_{\alpha }^{\beta }\right]$, а $J_1\mathrm{,}{J}_2\mathrm{,}{J}_3$ — обобщенные миноры матрицы $\left[{\Omega }_{\alpha }^{\beta }\right]$.

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

$\mathfrak{M}=\left(\mathcal{P}\mathrm{,}\mathcal{T}\mathrm{,}\mathcal{L}\mathrm{,}\nu \right)\mathrm{.}$ (4)

На первом месте стоит множество $\mathcal{P}$ - то же, что и в алгебраической структуре (1). Второй элемент, $\mathcal{T}$, определяет топологию, которая, в отличие от топологии тела, не определяется произвольно, а индуцируется метрической топологией открытых шаров в аффинном пространстве

$\mathcal{T}=\left\{\overline{)O}{C}_{o,{\left(e_{\alpha }\right)}_{\alpha =1}^3}\left(O\right) открыт в {ℝ}^3\right\}\mathrm{.}$

Третий элемент определяет максимальный атлас, который содержит декартово картрирование, порождаемое арифметизацией (2), и всевозможные непрерывные отображения открытых подмножеств ${ℝ}^3$ в ${ℝ}^3$. Последние определяют различные криволинейные координаты, используемые для описания мест в физическом пространстве. Последний элемент структуры (4) задает меру на многообразии таким образом, что мера линейных объектов (векторов, бивекторов, поливекторов) совпадает со значениями форм (3).

5°. Поскольку в рамках настоящей работы релятивистские эффекты не учитываются, изменения напряженно-деформированного состояния и связанных с ним полей дефектов можно рассматривать как зависимости от абсолютного времени. Оно количественно характеризуется координатами отдельного одномерного (хронометрического) многообразия $\mathfrak{T}$, топологически эквивалентного $ℝ$:

$\mathfrak{T}=(\left(a\mathrm{,}b)T\mathrm{,}R\mathrm{,}M\right)\mathrm{,} (a\mathrm{,}b)\subset ℝ\mathrm{.}$

Здесь $T$ – сужение стандартной топологии интервалов числовой оси на $\left(a\mathrm{,}b\right)$, $R$ – атлас, состоящий из всех непрерывных отображений $\left(a\mathrm{,}b\right)\to ℝ$, $M$ – мера, которая порождается мерой интервалов

$M\left(a\mathrm{,}b\right):=\phi \left(b\right)-\phi \left(a\right), \phi \in R\mathrm{.}$

Пространство-время $\mathfrak{S}$ как четырехмерное многообразие представляет собой глобально тривиальное расслоение над базой $\mathfrak{T}$ и типовыми слоями $\mathfrak{M}$, топология которого представляет собой топологию прямого произведения пространств $\mathfrak{T}$ и $\mathfrak{M}$, атлас — прямое произведение атласов этих многообразий, а мера задается как произведение мер.

6°. Тело $\mathfrak{B}$ в физическом пространстве представлено своей формой ${\mathcal{S}}_t$, которая может изменяться во времени, т. е. зависит от значения координаты $t$ на базе расслоения $\mathfrak{T}$. Форма определяется как образ многообразия $\mathfrak{B}$ относительно функции $\varkappa$, осуществляющей вложение многообразия $\mathfrak{B}$ в многообразие $\mathfrak{P}$. Требование к функции $\varkappa$ быть вложением исключает самопересечения образа и нарушение его непрерывности, наделяя тем самым форму свойствами подмногообразия физического многообразия $\mathfrak{P}$. В этой связи форму можно интерпретировать как ”непрерывное” множество мест материальных точек, в совокупности составляющих тело. Функции ${\varkappa }_t$ будем называть конфигурациями. Композиция двух конфигураций

$\gamma ={\varkappa }_1\circ {\varkappa }_2^{-1}$

определяет непрерывное отображение одной формы на другую:

$\gamma :\mathfrak{P}\supset {\mathcal{S}}_1\to {\mathcal{S}}_2\subset \mathfrak{P}\mathrm{.}$

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

Семейство форм в $\mathfrak{S}$ образует мировую трубку тела,

$\mathfrak{S}\supset \mathcal{U}=\widehat{\varkappa }=\left(\mathfrak{B}\mathrm{,}\mathfrak{T}\right),\widehat{\varkappa }\left(x\mathrm{,}t\right):={\varkappa }_t\left(x\right)$

которая представляет все формы как единый геометрический объект. Семейство деформаций ${\left\{{\gamma }_t\right\}}_{t\in \mathfrak{T}}$, послойно преобразующее сечения двух мировых трубок, определяемых одним и тем же временным слоем, позволяет однозначно восстановить одну из мировых трубок по другой. При этом одна из мировых трубок может быть ассоциирована с отсчетным состоянием, а вторая — с актуальным. Формы, мировые трубки, конфигурации и деформации тела, обладающего ненапряженной формой, иллюстрируются рис. 1. На этом рисунке сечение отсчетной трубки не изменяется во времени. В общем случае изменение сечения отсчетной трубки характеризует эволюцию отсчетного описания. Изменение сечения актуальной трубки характеризует наблюдаемую в физическом пространстве деформацию. Рис. 2 иллюстрирует эту ситуацию.

 

Рис. 1: Формы, мировые трубки, конфигурации и деформации тела, обладающего ненапряженной формой

   

Рис. 2: Формы, мировые трубки, конфигурации и деформации тела в общем случае эволюции отсчетной и актуальной форм

 

7°. Физическое состояние тела определяется из условия стационарности действия

$\delta \mathcal{A}=0,$

Действие $\mathcal{A}$ задается как скалярнозначная функция, аргументами которой служат две мировые трубки тела, первая — отсчетная трубка, вторая — актуальная:

$\mathcal{A}:\left({\mathcal{U}}_R\mathrm{,}\mathcal{U}\right)\to ℝ\mathrm{.}$

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

$\mathcal{A}:\left(\mathfrak{B}\mathrm{,}\mathfrak{T}\mathrm{,}{\widehat{\varkappa }}_R\mathrm{,}\widehat{\varkappa }\right)\to ℝ\mathrm{.}$

Число аргументов можно сократить до трех, если вместо тела использовать одну из его форм, ${\mathcal{S}}_i$, которую далее будем называть промежуточной, и двух обобщенных деформаций:

$\mathcal{A}:\left({\mathcal{S}}_i\mathrm{,}{\widehat{\gamma }}_R\mathrm{,}\widehat{\gamma }\right)\to ℝ\mathrm{.}$

На рис. 1, 2 промежуточная форма обозначена символом ${\mathcal{S}}_I$ и изображена синим цветом.

Существенного упрощения формулировки можно достичь, если ввести следующие предположения.

1. Действие может быть представлено в аддитивной форме

$\mathcal{A}\underset{\mathcal{U}}{=\int }\mathcal{L}\left(x\mathrm{,}{\widehat{\gamma }}_R\mathrm{,}\widehat{\gamma }\right)dv\left(x\right)$

Здесь, однако, подынтегральное выражение – плотность действия $\mathcal{L}$ – все еще очень сложно, поскольку второй и третий его аргументы представляют собой функции, определяющие соответствующие трубки, т. е. в такой формулировке плотность действия представляет собой функционал, действующий на функциональном пространстве обобщенных деформаций.

2. В силу принципа локализации обобщенные деформации могут быть заменены их формулами Тейлора первого порядка:

$\widehat{\gamma }(x+h)=\gamma \left(x\right)+{\hat{F}}_{x}h+o(\parallel h\parallel ), {\hat{F}}_x=D_x\widehat{\gamma }\mathrm{.}$

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

$\mathcal{L}\left(x\mathrm{,}{\widehat{\gamma }}_R\mathrm{,}\widehat{\gamma }\right)\approx \widehat{\mathcal{L}}\left(x\mathrm{,}{\widehat{\gamma }}_R\right(x),\widehat{\gamma }(x),{\hat{F}}_{xR}\mathrm{,}{\hat{F}}_x)$

3. В силу (нерелятивистского) принципа объективности плотность действия должна быть инвариантна относительно преобразований Галилея, т. е. трансляций и поворотов в слоях и сдвигов базы $\mathfrak{T}$. В этой связи $\mathcal{L}$ не может явно зависеть от $\widehat{\gamma }\left(x\right)$, ${\hat{F}}_{\mathrm{xR}}$ и, кроме того, зависимость ${\hat{F}}_{\mathrm{xR}}$, ${\hat{F}}_x$ усекается до зависимости только от симметричных частей этих тензоров, т. е.

$\mathcal{L}\left(x\mathrm{,}{\widehat{\gamma }}_R\mathrm{,}\widehat{\gamma }\right)\approx \widehat{\mathcal{L}(}x\mathrm{,}{F}_{xR}^{\text{s}ym}\mathrm{,}{F}_x^{\text{s}ym})$

Эквивалентную зависимость в более лаконичном (и более привычном) виде можно получить, если ввести тензоры деформаций

$C_{\mathrm{R}}:=F_{\mathrm{xR}}^{\text{T}}\circ {\mathrm{F}}_{\mathrm{xR}}\mathrm{,} \mathrm{C}:={\mathrm{F}}_{\mathrm{x}}^{\text{T}}\circ {\mathrm{F}}_{\mathrm{x}}\mathbf{.}$

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

$\mathcal{L}(x\mathrm{,}{I}_1^R\mathrm{,}\dots \mathrm{,}{I}_n^R;I_1\mathrm{,}\dots \mathrm{,}{I}_n),$

где $I_1^R\mathrm{,}\dots \mathrm{,}{I}_n^R$, $I_1\mathrm{,}\dots \mathrm{,}{I}_n$ – инварианты тензоров $C_{\mathrm{R}}$, $C$, причем в случае анизотропии материала n=3 и $I_1^R=C_R:I, I_2^R=\frac{1}{2}\left({(C_R:I)}^2-C_R^2:I\right)\mathrm{,} I_3^R=\frac{1}{3}\left({\left(C{}^{\mathrm{R}}\right)}^3:I-({I_1^R)}^3+3I_1^R{I}_2^R\right);$

$I_1=C:I, I_2=\frac{1}{2}\left({(C:I)}^2-C^2:I\right)\mathrm{,} I_3=\frac{1}{3}\left({\mathrm{C}}^3:\mathrm{I}-{\mathrm{I}}_1^3+3{\mathrm{I}}_1{\mathrm{I}}_2\right)\mathrm{,}$

а для анизотропного материала количество инвариантов возрастает. Так, для ортотропного материала, который, в основном используется в MEMS, добавляются два инварианта, совместных с вектором ориентации $A$ (или $A_R$) плоскости анизотропии:

$I_4^R=A_R\cdot C_R\cdot A_R\mathrm{,} I_5^R=A_R\cdot C_R\cdot C_R\cdot A_R; I_4=A\cdot C\cdot A\mathrm{,} I_5=A\cdot C\cdot C_R\cdot A\mathrm{.}$

Явные представления плотности действия могут быть получены как минимум двумя способами. Первый предполагает вывод соотношений из более детализированных физических теорий, опираясь на моделирование атомистических и квантомеханических представлений. Этот подход, интенсивно развиваемый в настоящее время, получил название ” первых принципов” или ab initio. Второй предполагает непосредственную аппроксимацию результатов испытаний некоторого представительного объема материала, который, как правило, в начале испытания приводится в некоторое однородное состояние. Оба подхода предполагают построение соотношений для достаточно малой части тела, причем внешнее воздействие, реакция на которое непосредственно описывается этими соотношениями, прикладывается к части, заранее приведенной в некоторое стандартное состояние. Задача моделирования деформирования тела — воспользоваться этими соотношениями и получить зависимости для физических полей во всем теле, представляя его как соединение достаточно большого числа достаточно малых частей (строго говоря, инфинитезимальных). Если удается найти такое соединение, в котором все части окажутся в стандартном состоянии, то задача существенно упрощается, а моделирование осуществляется классическими методами механики континуума. Такое соединение будет отвечать стандартной (свободной от напряжений) отсчетной форме. Проблема лишь в том, что в реальности таких форм не существует, поскольку во всех телах присутствуют внутренние источники напряжений — дефекты. В этом случае вместо одной отсчетной формы необходимо рассматривать континуальное семейство форм (рис. 3):

$\mathfrak{B}∍x\mapsto {\mathcal{U}}_x\mathrm{.}$ (5)

   

Рис. 3: Семейство отсчетных форм и семейство отсчетных трубок

 

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

 

${\gamma }_x{\overline{){\mathcal{S}}_R}}_{x}typ\mathrm{,}$

где символом $\text{t}yp$ обозначен одинаковый для всех точек элементарный объем-прототип. Градиент этих деформаций, вычисленный для каждой точки тела, определяет поле импланта K:

$\mathfrak{x}\mapsto \text{g}rad{\gamma }_{\mathfrak{x}}=:K^{-1}\mathrm{,}$

которое может быть использовано для задания метрики на $\mathfrak{B}$ посредством следующих соотношений:

$G(u,v):=g\left(K^{-1}u\mathrm{,}{K}^{-1}v\right)$ (6)

Поскольку поле импланта K в общем случае не является градиентом векторного поля, то построенная таким образом метрика определяет риманово пространство с ненулевой кривизной. Этой метрике отвечает (материальная) связность

${\nabla }_{\mathrm{u}}v=u^i\left({\partial }_i{v}^k+v^j{{\Gamma }_i^k}_j\right){\partial }_k=u^i{\nabla }_i{v}^k{\partial }_k\mathrm{,}$

${{\Gamma }_i^k}_j=\frac{G^{kl}}{2}\left({\partial }_i{G}_{jl}+{\partial }_j{G}_{il}-{\partial }_l{G}_{ij}\right)\mathrm{.}$

Кривизна построенной таким образом связности:

${R^i}_{jkl}={\partial }_k{\Gamma }_{lj}^i-{\partial }_l{\Gamma }_{kj}^i+{\Gamma }_{kr}^i{\Gamma }_{lj}^r-{\Gamma }_{lr}^i{\Gamma }_{kj}^r\mathrm{,}$

характеризует меру несовместности деформаций, или, иначе, плотность распределения дефектов, или же плотность распределения собственных источников напряжений. Все три определения с макроскопической точки зрения эквивалентны. Дифференциально-геометрическое описание является наиболее строгой математической характеристикой феномена собственного напряженного состояния, порождаемого распределенными в теле дефектами, на что мы уже указывали в начале статьи, ссылаясь на классические работы [13; 15–17]. Более подробно эти вопросы изложены, в частности, в [30–32].

   

Рис. 4: Неевклидова отсчетная форма и соответствующая трубка

 

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

Конечноэлементуную модель определим с помощью конформной тетраэдральной сетки в области ${\mathcal{S}}_I$, занимаемой промежуточной формой. При этом перемещения из промежуточной формы ${\mathcal{S}}_I$ в актуальную ${\mathcal{S}}_A$ (т. е. разность значений x деформации ${\gamma }_A$ и мест соответствующих точек в промежуточной форме X) аппроксимируется непрерывными функциями, линейными на каждом тетраэдре.

Пусть один из таких тетраэдров ${\mathcal{T}}_R$ с вершинами $P_1$, $P_2$, $P_3$, $P_4$ и ориентированным объемом

${\mathcal{V}}_R=\frac{1}{6}(P_4-P_1)\cdot \left(\left({\mathrm{P}}_{\mathrm{2}}\mathrm{-}{\mathrm{P}}_{\mathrm{1}}\mathrm{)}\mathrm{\times }\mathrm{(}{\mathrm{P}}_{\mathrm{3}}\mathrm{-}{\mathrm{P}}_{\mathrm{1}}\right)\right)$

трансформируется в результате деформации $x\left(X\mathrm{,}t\right)$ в тетраэдр $\mathcal{T}$ с вершинами $S_1$, $S_2$, $S_3$, $S_4$ и ориентированным объемом

$\mathcal{V}=\frac{1}{6}(S_4-S_1)\cdot (\left({\mathrm{S}}_2-{\mathrm{S}}_1)\times ({\mathrm{S}}_3-{\mathrm{S}}_1)\right)\mathrm{.}$

Символом $J$ обозначим якобиан преобразования, $J=\mathcal{V}/{\mathcal{V}}_R>0$.

Пусть ${\lambda }_1\left(X\right)$, ${\lambda }_2\left(X\right)$, ${\lambda }_3\left(X\right)$, ${\lambda }_4\left(X\right)$ — барицентрические координаты точки X:

$\left(\begin{array}{c}{\lambda }_1\left(X\right)\\ {\lambda }_2\left(X\right)\\ {\lambda }_3\left(X\right)\end{array}\right)=T^{-1}\cdot \left(X-P_1\right)\mathrm{,} T=\left(P_2-P_1 P_3-P_1 P_4-P_1\right)\mathrm{,}$

${\lambda }_4=1-{\lambda }_1-{\lambda }_2-{\lambda }_3\mathrm{,}$

где тензор T представлен ”декомпозицией по столбцам”, роль которых играют векторы $P_2-P_1,...,P_4-P_1$. Тогда координаты каждой точки X внутри тетраэдра перед деформацией и координаты ее положения x после деформирования могут быть представлены в форме

$X=\sum _{i=1}^4{\lambda }_i{P}_i\mathrm{,} x=\sum _{i=4}^4{\lambda }_i\left(X\right)S_i\mathrm{.}$

Используя формулу для обращения матрицы $3\times 3$:

$A^{-1}=\frac{1}{\det A}\left[\begin{array}{c}{({\overrightarrow{x}}_2\times {\overrightarrow{x}}_3)}^{\text{T}}\\ {({\overrightarrow{x}}_3\times {\overrightarrow{x}}_1)}^{\text{T}}\\ {({\overrightarrow{x}}_1\times {\overrightarrow{x}}_2)}^{\text{T}}\end{array}\right]\mathrm{,} A=\left[\begin{array}{ccc}{\overrightarrow{x}}_1& {\overrightarrow{x}}_2& {\overrightarrow{x}}_3\end{array}\right]\mathrm{,}$

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

$F=\sum _{i=1}^4{S}_i\otimes D_i\mathrm{,} D_i=\frac{{(-1)}^i}{6{\mathcal{V}}_R}\left(\left({\mathrm{P}}_{\mathrm{\{}\mathrm{i}\mathrm{+}\mathrm{2}\mathrm{\}}}\mathrm{-}{\mathrm{P}}_{\mathrm{\{}\mathrm{i}\mathrm{+}\mathrm{1}\mathrm{\}}}\mathrm{)}\mathrm{\times }\mathrm{(}{\mathrm{P}}_{\mathrm{\{}\mathrm{i}\mathrm{+}\mathrm{3}\mathrm{\}}}\mathrm{-}{\mathrm{P}}_{\mathrm{\{}\mathrm{i}\mathrm{+}\mathrm{1}\mathrm{\}}}\right)\right)\mathrm{,}$

где $\left\{i\right\}:=i-4⌊\frac{i-1}{4}⌋$ обозначает оператор индексирования узлов в диапазоне (1...4), $⌊...⌋$ — целая часть числа.

Локально искажения при деформации из промежуточной формы в актуальную характеризуются правым тензором деформаций Коши — Грина:

${\mathrm{F}}^{\mathrm{T}}\cdot \mathrm{F}=\sum _{\mathrm{i}\mathrm{,}\mathrm{j}=1}^4{\mathrm{S}}_{\mathrm{i}}\cdot {\mathrm{S}}_{\mathrm{j}}{\mathrm{D}}_{\mathrm{i}}\otimes {\mathrm{D}}_{\mathrm{j}}\mathrm{,}$

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

 

$C=K^T·F^T·F·K=\mathrm{\sum }_{\mathrm{i}\mathrm{,}\mathrm{j}\mathrm{=}\mathrm{1}}^{\mathrm{4}}{S}_i·S_j·D_i·K\otimes {\mathrm{D}}_{\mathrm{j}}·\mathrm{K}$

Инварианты тензора C находятся из соотношений:

null

В силу линейности функций формы упругий потенциал $W\left(C\right)$ постоянен внутри тетраэдра. При этом энергия упругого деформирования, запасаемая в тетраэдре, определяется равенством

${\mathcal{U}}_R={\mathcal{V}}_{R}W\left(C\right)$

Сила упругого взаимодействия, действующая на i-ю вершину тетраэдра ${\mathcal{T}}_R$, находится из соотношения (дискретного аналога формулы Дойля — Эриксена):

$f_{\mathrm{i}}^{\left({\mathcal{T}}_{\mathrm{R}}\right)}=-\frac{\partial {\mathcal{U}}_{\mathrm{R}}}{\partial {\mathrm{S}}_{\mathrm{i}}}=-\left|{\mathcal{V}}_{\mathrm{R}}\right|\sum _{\mathrm{s}=1}^5\frac{\partial \mathrm{W}}{\partial {\mathrm{I}}_{\mathrm{s}}}\frac{\partial {\mathrm{I}}_{\mathrm{s}}}{\partial {\mathrm{S}}_{\mathrm{i}}}\mathrm{.}$

Входящие в это соотношение производные инвариантов по узловым смещениям вычисляются по формулам:

$$\begin{gathered} \frac{\partial I_1}{\partial S_i}=2\sum _{n=1}^4{D}_n\cdot K\cdot K^{\text{T}}\cdot D_i{S}_n \\ \frac{{\displaystyle \partial I_2}}{{\displaystyle \partial S_i}}=\sum _{n=1}^4 D_n\cdot K\cdot \left(I_{1}I-C\right)\cdot K^{\text{T}}\cdot D_i{S}_n \\ \frac{{\displaystyle \partial I_3}}{{\displaystyle \partial S_i}}=2J|K{|}^2\frac{\partial J}{\partial S_i}, \frac{{\displaystyle \partial J}}{{\displaystyle \partial S_i}}=\frac{{(-1)}^i}{6{\mathcal{V}}_R}(S_{\{i+1\}}-S_{\{i+3\}})\times (S_{\{i+2\}}-S_{\{i+3\}}), \\ \frac{{\displaystyle \partial I_4}}{{\displaystyle \partial S_i}}=2D_i\cdot K\cdot A_1\sum _{n=1}^4 D_n\cdot K\cdot A_1{S}_n\mathrm{,} \\ \frac{{\displaystyle \partial I_5}}{{\displaystyle \partial S_i}}=2\sum _{n=1}^4 \{D_i K\cdot A_1{D}_{n}K\cdot C\cdot A_1+D_i\cdot K·C\cdot A_1{D}_n\cdot K\cdot A_1\}{S}_n\mathrm{.} \end{gathered}$$

Если положить K=I, то эти формулы преобразуются к известным соотношениям [33].

Производные плотности упругой энергии $\frac{\partial W}{\partial I_s}$ зависят от выбранного гиперупругого потенциала $W$. Четыре варианта формулировки $W$, наиболее часто используемые при анализе конечных деформаций гиперупругих тел, приведены в таб. 1. Выражения для производных $\frac{\partial W}{\partial I_s}$ в явной форме приведены в таб. 2. Нумерация строк таб. 1, указанная в первой колонке, соответствует нумерации потенциалов в таб. 2.

 

Таблица 1

Table 1

 Гиперупругие потенциалы

Наименование	Потенциал
Муни — Ривлина (1)	$\frac{\mu }{4}\left[(1-\beta )\left(\frac{I_1}{I_3^{1/3}}-3\right)+(\beta +1)\left(\frac{I_2}{I_3^{2/3}}-3\right)\right]+\frac{3\lambda +2\mu }{6}{\left(\sqrt{I_3}-1\right)}^2$
Синьорини (2)	$\sqrt{I_3}\left(\frac{p}{2}\left(\frac{I_1}{I_3}+1\right)+\frac{\lambda +\mu -p}{8}\left(\frac{I_2^2}{I_3^2}+3\right)+\frac{p-3\lambda -\mu }{4}\left(\frac{I_2}{I_3}-1\right)\right)-\mu -p$
Мурнагана (3)	$\begin{array}{l}\frac{\lambda +2\mu }{8}({I_1-3)}^2-\frac{\mu }{2}(3-2I_1+I_2)+a({I_1-3)}^2+b(I_1-3)(3-2I_1+I_2)+c(I_1-I_2+I_3-1)\\ \end{array}$
Блейтца и Ко (4)	$\frac{\mu }{4}(\beta +1)\left(I_1+\frac{1}{\nu }\left(I_3^{-\nu }-1\right)-3\right)+\frac{\mu }{4}(1-\beta )\left(\frac{I_2}{I_3}+\frac{1}{\nu }\left(I_3^{\nu }-1\right)-3\right)$

 

Таблица 2

Table 2

Производные гиперупругих потенциалов

	$\frac{\partial W}{\partial I_1}$	$\frac{\partial W}{\partial I_2}$	$\frac{\partial W}{\partial I_3}$
1	$\frac{(1-\beta )\mu }{4\sqrt[3]{I_3}}$	$\frac{(\beta +1)\mu }{4I_3^{2/3}}$	$\frac{(\beta -1)I_1\sqrt[3]{I_3}\mu -2(\beta +1)I_2\mu +\left(I_3^{5/3}-I_3^{7/6}\right)(6\lambda +4\mu )}{12I_3^{5/3}}$
2	$\frac{p}{2\sqrt{I_3}}$	$\frac{I_2(\lambda +\mu -p)+I_3(p-3\lambda -\mu )}{4I_3^{3/2}}$	$\frac{I_3\left(I_3\right(9\lambda +5\mu -p)-4I_{1}p)+3I_2^2(p-\lambda -\mu )+2I_2{I}_3(3\lambda +\mu -p)}{16I_3^{5/2}}$
3	$2\left(a-b+\frac{\lambda +2\mu }{4}\right)(I_1-3)+b(3-2 I_1+I_2)+\mu +c$	$b(I_1-3)-c-\frac{\mu }{2}$	c
4	$\frac{1}{4}(\beta +1)\mu$	$-\frac{(\beta -1)\mu }{4I_3}$	$\begin{array}{r}\frac{\mu }{4}{I}_3^{-\nu -2}\left(\right(\beta -1)\left(I_2{I}_3^{\nu }-I_3^{2\nu +1}\right)-(\beta +1)I_3\end{array})$

 

Таким образом, силы в вершинах тетраэдров могут быть найдены в виде аналитических (алгебраических) зависимостей от позиций вершин, задаваемых векторами $S_i$, $i\in (1,\dots ,4)$.

Сила, действующая на i-ю вершину тетраэдра со стороны внешних силовых полей (объемной плотностью b и поверхностной, t), определяется следующим интегралом:

 

$f_{\mathrm{i}\mathrm{,}\text{e}\mathrm{xt}}^{\left({\mathcal{T}}_{\mathrm{R}}\right)}=\underset{{\mathrm{\Gamma }}_{\mathrm{\sigma }}\left(\mathrm{t}\right)\cap \mathcal{T}}{\int }t{\mathrm{\lambda }}_{\mathrm{i}}\mathrm{ds}+\underset{\mathcal{T}}{\int }{\mathrm{\lambda }}_{\mathrm{i}}\mathrm{d\Omega }\mathrm{.}$

Здесь $\mathcal{T}$ — регион, занимаемый тетраэдром в составе актуальной формы.

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

$\sum _{\mathcal{T}\in {\Sigma }_i}\left({\mathrm{f}}_{\mathrm{i}}^{\left({\mathcal{T}}_{\mathrm{R}}\right)}+{\mathrm{f}}_{\mathrm{i}\mathrm{,}\text{e}\mathrm{xt}}^{\left({\mathcal{T}}_{\mathrm{R}}\right)}\right)=\mathrm{0,}$ (7)

где ${\Sigma }_i$ обозначает множество тетраэдров с вершиной в узле номер i.

Эта система содержит значения тензорного поля импланта K в качестве параметра. Если положить K=I, то промежуточная форма окажется свободной от напряжений, и решение системы (7) будет определять напряженно-деформированное состояние актуальной формы в смысле классической гиперупругости. Отличие K от единичного тензора характеризует внутренние источники деформаций, роль которых играют распределенные дефекты. В этом случае с каждым абстрактным тетраэдром конечноэлементной модели связаны три его метрических представления: ненапряженный тетраэдр, отсчетный тетраэдр в составе промежуточной формы и актуальный тетраэдр в составе актуальной формы (рис. 5). Все три тетраэдра в отдельности являются регионами физического пространства, и потому для их описания не требуется неевклидовых представлений. Вместе с тем совокупности тетраэдров, образующих конечноэлементные сети, могут быть вложены в физическое (евклидово) пространство только для промежуточной и актуальной форм. Отсчетные тетраэдры в общем случае могут быть вложены в физическое пространство в составе конечноэлементной сети, индивидуальной для каждого из них. В этом состоит аналогия дискретизированного состояния с семейством отсчетных форм (5). Попытка погрузить совокупность отсчетных тетраэдров, не нарушая совпадения узлов, в физическое пространство окажется безуспешной, поскольку длины ребер соседних тетраэдров не будут совпадать в евклидовой метрике. Надлежащая ”подстройка” метрики (6) над многообразием, получаемым из аффинного пространства за счет ” стирания” метрики ассоциированного с ним трансляционного евклидова пространства, приводит к построению риманового пространства, в которое все ненапряженные тетраэдры вкладываются в совокупности без разрывов конечноэлементной сети.

 

Рис. 5: Неевклидова отсчетная форма, семейство отсчетных форм, промежуточная и актуальная формы в конечноэлементном представлении

 

Несколько слов следует сказать о вычислении производных второго и третьего порядков, требуемых для вычисления связности и кривизны по полю импланта, определяемого из конечноэлементного расчета. Непосредственное дифференцирование кусочно-линейной аппроксимации ничего, кроме нулевых (либо несуществующих) значений, не даст, и для вычисления производных требуется иной подход. При построении алгоритма производные вычислялись как проекции на специальные базисные функции, т. е. процедура дифференцирования заменялась интегрированием по объему тела. Базисные функции определялись как элементы системы, биортогональной к элементарным мономам ($\mathrm{1,}x\mathrm{,}{x}^2\mathrm{,}xy\mathrm{,}xz\mathrm{,}\dots$). Биортогональность понималась в смысле квадратичной метрики пространства функций, задаваемых над областью, занимаемой промежуточной формой.

9°. Если поле импланта K полагать заданным, то в результате решения алгебраической системы (7) можно определить напряженно-деформированное состояние отсчетной и актуальной форм, вызванное действием внешних полей и внутренних источников напряжений, количественно полностью характеризуемых полем K. Такая задача в нелинейной теории упругости обычно называется задачей о (несовместной) дисторсии, и ей посвящены многочисленные работы, начиная с работ Вольтерры. Более сложной является задача о переменной дисторсии, которая эволюционирует в процессе деформирования тела. В этом случае поле K изменяется во времени, и это изменение следует находить из дополнительного уравнения, характеризующего физические процессы типа миграции дефектов или какого-либо иного физического механизма, с макроскопической точки зрения описываемого как вязко-пластичность. Для того чтобы конкретизировать этот механизм, следует выбрать некоторое уравнение эволюции для поля . В рамках настоящей статьи, не погружаясь в физические детали теории дислокаций, которые можно найти в работах [34–36], будем полагать, что эволюция отсчетной формы определяется в рамках дислокационной модели Александера — Хаасена — Сумино. Поле этом K определяется из эволюционного уравнения, связанного с эволюционным уравнением для плотности дефектов N:

$\begin{array}{l}K=AN{⟨\sqrt{J_2}-B\sqrt{N}⟩}^p\frac{\mathrm{S}}{\sqrt{J_2}}\mathrm{,}\\ \dot{N}=C{⟨\sqrt{J_2}-A\sqrt{N}⟩}^{p+\lambda }\mathrm{.}\end{array}$ (8)

 Здесь A, B, C, p, λ — эмпирические константы (более точно, функции температуры), причем A, C характеризует подвижность дефектов, а B — упрочнение материала, связанное с их перераспределением; S – девиаторная часть тензора напряжений Коши, $J_2=S:S^{\text{T}}/2$, а угловые скобки определяют ”одностороннюю связь”:

$⟨x⟩=\{x, x>00, x\le 0.$

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

 

Рис. 6: Плотность распределения дефектов в начале и в конце процесса контактного взаимодействия

 

Рис. 7: Распределения характеристик несовместимости деформаций

 

10°. При действии на тело и его границу заданных силовых полей на каждом временном слое решения нелинейной системы уравнений (7) полностью определяют напряжения и деформации. При этом на части узлов может быть задано смещение, в результате чего эти узлы исключаются из множества неизвестных и число уравнений сокращается. В этой связи учет контактного взаимодействия с жестким штампом при фиксированной области контакта выполняется по описанному выше алгоритму, который лишь упрощается за счет уменьшения числа неизвестных. Конечно, при этом задается смещение штампа, а результирующие сила и момент вычисляются по ним. Если из условий расчета известны усилия, приложенные к штампу, то для моделирования контактного взаимодействия требуется выполнить итерационный процесс, в ходе которого задаются пробные смещения, которые затем корректируются из условия минимизации фактических усилий на штампе и их заданных значений. Однако если область контакта заранее неизвестна, то, даже при кинематическом задании внедрения штампа указать смещения узлов сразу не удается. Здесь также требуется итерационный вычислительный процесс, на первом шаге которого задаются смещения узлов (по профилю штампа), актуальные позиции которых попадают в пробную область контакта. При этом позиции части узлов, не попавших в пробную область, могут оказаться внутри области пространства, занимаемой штампом. Такое ”проникновение” означает, что пробная область слишком мала и её требуется увеличить. Добавляя соседние узлы к пробной области, повторяем расчет до тех пор, пока не исключим ”проникновение”. Если внедрение штампа задано не кинематически, а посредством результирующей силы и момента, то итерационный процесс происходит по силовым и контактным параметрам — результирующим усилиям и форме пробной области. В случае осесимметричной формы штампа итерационный процесс существенно упрощается, поскольку пробная область задается скалярным параметром — ее диаметром. Вычисления, проводимые в настоящей работе, предполагают такую симметрию штампов.

Заключение

С помощью изложенного выше алгоритма были произведены тестовое моделирование эволюции напряженно-деформированного состояния и распределения дефектов в параллелепипеде при внедрении системы трех жестких штампов со сферической формой поверхности контакта. Материальные характеристики, используемые при моделировании, соответствуют кремнию (в рамках потенциала Муни — Ривлина при $\beta =0$ и не учитывая фактора анизотропии), а размеры параллелепипеда $2\mu m\times 1\mu m\times 1\mu m$. Диаметр каждого штампа составил 1/10 от наибольшего размера тела, т. е. $0.1\mu$, расстояние между штампами $0.2\mu$. Параметры эволюционных уравнений для дефектов взяты из [37]. Программа нагружения штампов – линейная по времени, до достижения максимального напряжения, равного пределу пластичности. Длительность процесса нагружения – 100 мин, интервал одного временного слоя – 1 мин. Начальное распределение дефектов полагалось квазипериодическим, графически пространственное распределение N в среднем сечении параллелепипеда показано на рис. 6а. Финальное распределение — на рис. 6b. Поле К в начале расчета полагалось единичным, что соответствует отсутствию собственных (несовместных) деформаций в начале нагружения. Для характеристики изменения К в финале процесса на рис. 7а показано распределение в среднем сечении скалярной величины – второго инварианта его девиаторной части. (такой выбор продиктован ассоциацией с определяющей характеристикой поля напряжений в классической пластичности). Рис. 7b иллюстрирует распределение скалярной кривизны материальной связности в среднем сечении параллелепипеда.

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

Работа выполнена при финансовой поддержке гранта РФФИ 20-01-00400 А.

About the authors

Tatiana N. Lycheva

Ishlinsky Institute for Problems in Mechanics RAS

Author for correspondence.
Email: taniko81@mail.ru
ORCID iD: 0000-0002-8802-022X

tribology laboratory engineer

Russian Federation, Moscow

Sergey A. Lychev

Ishlinsky Institute for Problems in Mechanics RAS

Email: lychevsa@mail.ru
ORCID iD: 0000-0001-7590-1389

Doctor of Physical and Mathematical Sciences, associate professor, leading researcher, Laboratory of Mechanics of Technological Processes

Russian Federation, Moscow

References

Supplementary files