Bending of a circular disk: from cylinder to ultrathin membrane

Full Text

Введение

Расчет и оптимизация физико-механических параметров круглых пластин и мембран часто используются при разработке разнообразных микросистем. Как правило, при проектировании микроэлектронных устройств основное внимание уделяется электромагнитным полям и связанным с ними параметрам. Однако при разработке микроэлектромеханических систем (МЭМС) необходимо также учитывать механические деформации и напряжения. В зависимости от назначения упругие элементы МЭМС могут представлять собой как пластины значительной относительной толщины (например, сенсоры акселерометров), так и ультратонкие мембраны (например, газоанализаторы) [1; 2]. При таком разбросе параметров выбор уравнений и их решений для описания напряженно-деформированного состояния (НДС) в каждом конкретном случае является неочевидным. Ситуация становится более сложной при переходе к ультратонким изгибаемым элементам, изгиб которых сильно зависит от натяжений в срединной плоскости, причем эта зависимость существенно нелинейная [3]. Кроме того, в силу особенностей технологических процессов изготовления ультратонких элементов (напыления, травления, роста), в них возникают значительные остаточные напряжения [4-6]. В этой связи выбор адекватной модели и решения соответствующей начально-краевой задачи представляет актуальный вопрос в современной теории расчета тонкостенных систем. Ему посвящена настоящая статья.

Основным результатом является “шкала относительных толщин”, на которой отмечены интервалы применимости четырех математических моделей:

• линейной трехмерной теории упругости;
• теории изгиба Кирхгофа;
• нелинейной теории изгиба Феппля Кармана;
• нелинейной теории гиперупругих мембран Адкинса.

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

1. Деформируемый диск

В статье рассматривается осесимметричное деформирование кругового диска толщиной $h$ и радиуса $R$, выполненного из упругого материала. Боковая поверхность диска жестко закреплена, а его основания свободны от напряжений (рис. 1.1). Эту задачу целесообразно рассматривать в цилиндрических координатах $\left(r\mathrm{,}\phi \mathrm{,}z\right)$, которые связаны с декартовыми координатами $\left(x\mathrm{,}y\mathrm{,}z\right)$ следующими соотношениями^[1]:

$r=\sqrt{x^2+y^2}\mathrm{,}\phi =arctg\left(y\mathrm{,}x\right)\mathrm{,}z=z\mathrm{.}$

Рис 1.1. Расчетная схема пластины

Fig. 1.1. The plate design scheme

Начало координат $O$ поместим в центре нижнего основания диска, а ось $OZ$ направим вдоль его оси так, чтобы координаты точек диска принадлежали множеству:

$\mathcal{D}=\{\left(r\mathrm{,}\phi \mathrm{,}z\right)\in {ℝ}^3:r\in ]0$$R[\mathrm{,}\phi \in \left]\mathrm{0,}2\pi \right[\mathrm{,}z\in \left]\mathrm{0,}h\right[$

Векторные поля будем представлять в виде разложений по локальным базисам $\left({\text{e}}_r\mathrm{,}{\text{\hspace{0.05em}e}}_{\phi }\mathrm{,}{\text{\hspace{0.05em}e}}_z\right)$, связанным с ортонормированным базисом $\left(\text{i}\mathrm{,}\text{\hspace{0.05em}j}\mathrm{,}\text{\hspace{0.05em}k}\right)$ преобразованиями:

${\text{e}}_r\text{=i}\cos \phi +\text{j}\sin \phi \mathrm{,}{\text{\hspace{1em}e}}_{\phi }\text{=}-\text{i}\sin \phi +\text{j}\cos \phi \mathrm{,}{\text{\hspace{1em}e}}_z\text{=k}\mathrm{.}$

В частности, векторное поле перемещений точек диска $\text{u}$ имеет вид:

$\text{u=}u{\text{e}}_r+w{\text{e}}_z\mathrm{,}$

где $u\mathrm{,}w$ функции, представляющие размерные радиальные и осевые перемещения:

$u:\left]\mathrm{0,}R\right[\times \left]\mathrm{0,}h\right[\ni \left(r\mathrm{,}z\right)\mapsto u\left(r\mathrm{,}z\right)\in ℝ\mathrm{,}$$\hspace{1em}w:\left]\mathrm{0,}R\right[\times \left]\mathrm{0,}h\right[\ni \left(r\mathrm{,}z\right)\mapsto w\left(r\mathrm{,}z\right)\in ℝ\mathrm{.}$

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

$\mathcal{P}=\mathcal{P}\left(r\right){\text{e}}_z\text{=}\left[p-\frac{p{R}^2}{R^2-a^2}\theta \left(r-a\right)\right]{\text{e}}_z\mathrm{,}$$\underset{0}{\overset{h}{\int }}\underset{0}{\overset{2\pi }{\int }}\underset{0}{\overset{R}{\int }}\mathcal{P}\left(r\right)rdrd\phi dz=0$ (1.1)

Рис. 1.2. Распределение внешних объемных сил

Fig. 1.2. Distribution of external volume forces

Здесь $p$ объемная плотность заданного (внешнего) силового поля, $a$ кольцевой области, $\theta$ Хевисайда. Нагрузка на основной части диска может быть ассоциирована, например, с внешним электростатическим воздействием, а нагрузка в кольцевой области моделирует опорную реакцию. Такая специфическая форма нагружения позволяет в дальнейшем получить более простые разложения и более эффективные для вычисления аналитические формулы, а погрешность, связанная с неучетом особенностей конкретной реализации краевых условий, незначительна в силу принципа Сен-Венана [7].

В рамках настоящей статьи анализируются четыре математические модели, одна из которых рассматривает материал как гиперупругий (с потенциалом Муни), а остальные линейно-упругий. Для того чтобы эти модели были согласованы, характеристики упругого и гиперупругого материалов должны быть связаны между собой. Это можно сделать, если записать гиперупругий потенциал Муни [8] через параметры Ламе $\mu \mathrm{,}\lambda$:

$\mathcal{W}=\frac{\mu }{4}\left[\left(1-\beta \right)\left(\frac{I_1}{I_3^{{\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}}}-3\right)+\left(1+\beta \right)\left(\frac{I_2}{I_3^{2/3}}-3\right)\right]+\frac{2\mu +3\lambda }{6}{\left(\sqrt{I_3}-1\right)}^2\mathrm{,}$ (1.2)

где $\beta$ материальная константа второго порядка, а $I_1\mathrm{,}{I}_2\mathrm{,}{I}_3$ инварианты левого тензора деформаций Коши-Грина²:

$\text{B=I}+\nabla {\text{u}}^{\text{T}}+\nabla \text{u}+\nabla \text{u}\cdot \text{}\nabla {\text{u}}^{\text{T}}\mathrm{,}$

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

$I_1\text{=}3+2q_1+q_2\mathrm{,}$

$I_2\text{=}3+4q_1+\left(2q_1^2+q_2-q_3\right)+2\left(q_1{q}_2-q_4\right)+\frac{1}{2}\left(q_2^2-q_6\right)\mathrm{,}$

$I_3\text{=}1+2q_1+\left(2q_1^2-q_3\right)+\frac{2}{3}\left(2q_1^3-3q_1{q}_3+q_5\right)+$

$+\frac{1}{2}\left(4q_1^2{q}_2-8q_1{q}_4-q_2^2-2q_2{q}_3+q_6+4q_7+2q_8\right)+\left(q_1{q}_2^2-q_1{q}_6-2q_2{q}_4+2q_9\right)+$

$+\frac{1}{6}\left(q_2^3-3q_2{q}_6+2q_{10}\right)\mathrm{,}$

где

$\text{}{q}_1\text{=tr\hspace{0.05em}}\nabla \text{u}\mathrm{,}{q}_2=\nabla \text{u:}\nabla \text{u}\mathrm{,}{q}_3=\nabla {\text{u}}^{\text{T}}\text{:}\nabla \text{u}\mathrm{,}$

$\text{}{q}_4\text{=}\left(\nabla \text{u}\cdot \text{}\nabla \text{u}\right)\text{:}\nabla \text{u}\mathrm{,}{q}_5=\left(\nabla \text{u}\cdot \text{}\nabla \text{u}\right):\text{}\nabla {\text{u}}^{\text{T}}\mathrm{,}{q}_6=\left(\nabla {\text{u}}^{\text{T}}\text{}\cdot \text{}\nabla \text{u}\right):\text{}\left(\nabla {\text{u}}^{\text{T}}\text{}\cdot \text{}\nabla \text{u}\right)\mathrm{,}$

$\text{}{q}_7\text{=}\left(\nabla \text{u}\cdot \text{}\nabla \text{u}\right)\text{:}\left(\nabla {\text{u}}^{\text{T}}\text{}\cdot \text{}\nabla \text{u}\right)\mathrm{,}{q}_8=\left(\nabla {\text{u}}^{\text{T}}\text{}\cdot \text{}\nabla \text{u}\right):\text{}\left(\nabla \text{u}\cdot \text{}\nabla {\text{u}}^{\text{T}}\right)\mathrm{,}{q}_9=\left(\nabla \text{u}\cdot \text{}\nabla {\text{u}}^{\text{T}}\text{}\cdot \text{}\nabla \text{u}\right):\text{}\left(\nabla {\text{u}}^{\text{T}}\text{}\cdot \text{}\nabla \text{u}\right)\mathrm{,}$

$\text{\hspace{0.05em}\hspace{0.05em}}{q}_{10}\text{=}\left(\nabla \text{u}\cdot \text{}\nabla {\text{u}}^{\text{T}}\text{}\cdot \text{}\nabla \text{u}\right)\text{:}\left(\nabla \text{u}\cdot \text{}\nabla {\text{u}}^{\text{T}}\text{}\cdot \text{}\nabla \text{u}\right)\mathrm{.}$

Для линеаризации записанных выше соотношений представим градиент перемещений в виде произведения приведенного градиента $\tilde{\nabla }\text{u}$ и малого параметра, в качестве которого возьмем норму тензора $\nabla \text{u}$:

$\nabla \text{u=}ϵ\tilde{\nabla }\text{u}\mathrm{,}ϵ=||\nabla \text{u}||,\tilde{\nabla }\text{u=}\frac{\nabla \text{u}}{||\nabla \text{u}||}\mathrm{,}$

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

$\frac{I_1}{I_3^{1/3}}\text{=}3+\nabla {\text{u}}^{\text{T}}\text{:}\nabla \text{u}+\nabla \text{u:}\nabla \text{u}-\frac{2}{3}{\left(\text{tr\hspace{0.05em}}\nabla \text{u}\right)}^2+o\left({ϵ}^2\right)\mathrm{,}$

$\frac{I_2}{I_3^{2/3}}\text{=}3+\nabla {\text{u}}^{\text{T}}\text{:}\nabla \text{u}+\nabla \text{u:}\nabla \text{u}-\frac{2}{3}{\left(\text{tr\hspace{0.05em}}\nabla \text{u}\right)}^2+o\left({ϵ}^2\right)\mathrm{,}$ (1.3)

${\left(\sqrt{I_3}-1\right)}^2={\left(\text{tr\hspace{0.05em}}\nabla \text{u}\right)}^2+o\left({ϵ}^2\right)\mathrm{.}$

Подставляя (1.3) в (1.2), можно получить линейно-упругий потенциал с точностью до $o\left({ϵ}^2\right)$:

$\mathcal{W}=\frac{\mu }{2}\left(\nabla {\text{u}}^{\text{T}}\text{:}\nabla \text{u}+\nabla \text{u:}\nabla \text{u}\right)+\frac{\lambda }{2}{\left(\text{tr\hspace{0.05em}}\nabla \text{u}\right)}^2+o\left({||\nabla \text{u}||}^2\right)=\mu \left(\epsilon :\text{}\epsilon \right)+\frac{\lambda }{2}{\left(\text{tr\hspace{0.05em}}\epsilon \right)}^2+o\left({ϵ}^2\right)\mathrm{,}$

где $\epsilon$ малых деформаций, представляющий линейную часть тензора деформаций Грина-Венана $\text{E}$, который, в свою очередь, выражается через правый тензор Коши-Грина $\text{C}$:

$\text{E=}\frac{1}{2}\left(\text{C}-\text{I}\right)\text{=}\epsilon +\text{o}\left(\nabla \text{u}\right)\mathrm{,}\epsilon =\frac{1}{2}\left(\nabla {\text{u}}^{\text{T}}+\nabla \text{u}\right)\mathrm{.}$

Компоненты тензора $\epsilon$связаны с компонентами вектора $\text{u }$перемещений соотношениями Коши, которые в случае осевой симметрии записываются в цилиндрических координатах следующим образом³:

$\left[\epsilon \right]=\left(\begin{array}{ccc}\frac{\partial u}{\partial r}& 0& \frac{1}{2}\left(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial r}\right)\\ 0& \frac{u}{r}& 0\\ \frac{1}{2}\left(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial r}\right)& 0& \frac{\partial w}{\partial z}\mathrm{.}\end{array}\right)\mathrm{.}$ (1.4)

Для построения решения и его последующего анализа целесообразно ввести безразмерные переменные, а именно

• пространственные переменные: $\tilde{r}=r/R\mathrm{,}\widetilde{z=}z/R$
• перемещения: $\widetilde{u=}u/R\mathrm{,}\widetilde{w=}w/R$
• масштабные параметры: $\tilde{R}=1,\widetilde{h=}h/R$
• модули Ламе: $\widetilde{\mu =1,}\tilde{\lambda }\text{=}\lambda /\mu$
• параметры нагружения: $\tilde{p}=pR/\mu \mathrm{,}\widetilde{a=}a/R\mathrm{.}$

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

При решении прикладных задач используются различные совокупности материальных констант, в частности, в трехмерных задачах теории упругости традиционно используются модули Ламе, а в технической механике Юнга $E$и коэффициент Пуассона $\nu$. Для единообразия в настоящей статье все материальные константы выразим через один безразмерный параметр $k$:

$\tilde{\lambda }=k-\mathrm{1,}\tilde{E}\text{=}\frac{E}{\mu }\text{=}\frac{3k-1}{k}\mathrm{,}\nu =\frac{k-1}{2k}\mathrm{.}$

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

${\tilde{\sigma }}_{eqv}=\frac{1}{\sqrt{2}}\sqrt{{\left({\tilde{\sigma }}_{rr}-{\tilde{\sigma }}_{\phi \phi }\right)}^2+{\left({\tilde{\sigma }}_{\phi \phi }-{\tilde{\sigma }}_{zz}\right)}^2+{\left({\tilde{\sigma }}_{zz}-{\tilde{\sigma }}_{rr}\right)}^2+6{\tilde{\sigma }}_{rz}^2}\mathrm{,}$ (1.5)

где ${\tilde{\sigma }}_{rr}\mathrm{,}{\tilde{\sigma }}_{\phi \phi }\mathrm{,}{\tilde{\sigma }}_{zz}\mathrm{,}{\tilde{\sigma }}_{rz}$ тензора напряжений $\sigma$.

В дальнейшем все соотношения и результаты будут представлены в безразмерном виде, поэтому, чтобы не усложнять систему обозначений, знак “тильда” будем опускать.

2. Решение в трехмерной постановке

Рассмотрим постановку задачи о НДС круговой плиты и ее решение в рамках линейной теории упругости. Задача об упругом цилиндре является классическим примером, допускающим аналитическое решение, различные варианты которого развивались с конца XIX века (решения Похгаммера [10], Кри [11] и их последователей⁴). Основная сложность в построении таких решений заключается в удовлетворении краевым условиям, форма которых диктуется не изначальной постановкой задачи, а возможностью разделить переменные и свести задачу для уравнений с частными производными к независимым двухточечным задачам для обыкновенных дифференциальных уравнений. В решениях Похгаммера краевые условия специального вида, соответствующие гладко-жесткому контакту (roller contact), задаются на основаниях цилиндра, что не позволяет применить их непосредственно к рассматриваемой в статье проблеме. Однако их модификация, предложенная Файлоном [17] и развитая Сайто [18; 19], позволяет несколько иначе разделить краевые задачи на две двухточечные задачи, перенеся обременительные гладко-жесткие условия на боковую поверхность и предоставив свободу выбора условий на основаниях. При этом условия на боковой поверхности цилиндра как нельзя лучше соответствуют жесткому закреплению в терминах теории пластин, а возможная неточность его реализации (например, можно пожелать абсолютно жесткого закрепления точек цилиндрической поверхности, при котором все компоненты перемещений на нем обращаются в ноль) играет исчезающе малую роль при уменьшении относительной толщины в силу принципа Сен-Венана.

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

${\nabla }^{2}u-\frac{u}{r^2}+k\left(\frac{{\partial }^{2}u}{\partial r^2}+\frac{1}{r}\frac{\partial u}{\partial r}-\frac{u}{r^2}+\frac{{\partial }^{2}w}{\partial r\partial z}\right)=0,$${\nabla }^{2}w+k\left(\frac{{\partial }^{2}u}{\partial r\partial z}+\frac{1}{r}\frac{\partial u}{\partial z}+\frac{{\partial }^{2}w}{\partial z^2}\right)=-\mathcal{P}\left(r\right)\mathrm{,}$ (2.1)

${\left.{\sigma }_{zz}\right|}_{z=0,h}{\left.=\left(k+1\right)\frac{\partial w}{\partial z}+\left(k-1\right)\left(\frac{\partial u}{\partial r}+\frac{u}{r}\right)\right|}_{z=0,h}\text{=0,}$ (2.2)

${\left.{\sigma }_{rz}\right|}_{z\text{\hspace{0.05em}=0,}h}\text{}{\left.\text{=}\frac{\partial u}{\partial z}+\frac{\partial w}{\partial r}\right|}_{z=0,h}\text{=0, \hspace{1em}\hspace{1em}}{\left.{\sigma }_{rz}\right|}_{r=1}{\left.\text{=}\frac{\partial u}{\partial z}+\frac{\partial w}{\partial r}\right|}_{r=1}\text{=0,\hspace{1em}\hspace{1em}}{\left.u\right|}_{r=1}\text{=0,}$

где ${\nabla }^2$Лапласа в цилиндрических координатах:

${\nabla }^2=\frac{{\partial }^2}{\partial r^2}+\frac{1}{r}\frac{\partial }{\partial r}+\frac{{\partial }^2}{\partial z^2}\mathrm{.}$

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

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

$u=f_1\left(r\right)g_1\left(z\right)\mathrm{,}w=f_2\left(r\right)g_2\left(z\right)$

при соблюдении условий⁵

$f_2^{\text{'}}=-\alpha f_1\mathrm{,}{f}_1^{\text{'}}+\frac{f_1}{r}=\alpha f_2\mathrm{,}$ (2.3)

где $\alpha$константа, то система уравнений (2.1) может быть приведена к виду:

$-\left(f_1^{\text{'}\text{'}}+\frac{f_1^{\text{'}}}{r}-\frac{f_1}{r^2}\right)\frac{1}{f_1}=\frac{g_1^{\text{'}\text{'}}}{g_1}\frac{1}{k+1}-\frac{\alpha k}{k+1}\frac{g_2^{\text{'}}}{g_1}\text{=}{\lambda }_1\mathrm{,}$(2.4)

$-\left(f_2^{\text{'}\text{'}}+\frac{f_2^{\text{'}}}{r}\right)\frac{1}{f_2}=\frac{g_2^{\text{'}\text{'}}}{g_2}\left(k+1\right)+\alpha k\frac{g_1^{\text{'}}}{g_2}=\text{}{\lambda }_2\mathrm{,}$

в котором ${\lambda }_1\mathrm{,}{\lambda }_2$ обозначают константы разделения. В действительности эти константы оказываются равными, поскольку из шести уравнений (2.3), (2.4) только четыре являются независимыми, и если из системы (2.3) поочередно исключать $f_1\mathrm{,}{f}_2$, тогда получим два уравнения, подобных (2.4):

$f_1^{\text{'}\text{'}}+\frac{f_1^{\text{'}}}{r}-\frac{f_1}{r^2}=-{\alpha }^2{f}_1\mathrm{,}{f}_2^{\text{'}\text{'}}+\frac{f_2^{\text{'}}}{r}=-{\alpha }^2{f}_2\mathrm{.}$

Следовательно, ${\lambda }_1={\lambda }_2={\alpha }^2$. Таким образом, имеем четыре независимых дифференциальных уравнения:

$r^2{f}_1^{\text{'}\text{'}}+r{f}_1^{\text{'}}+\left({\alpha }^2{r}^2-1\right)f_1\text{=}\mathrm{0,}{r}^2{f}_2^{\text{'}\text{'}}+r{f}_2^{\text{'}}+{\alpha }^2{r}^2{f}_2=0,$ (2.5)

$g_1^{\text{'}\text{'}}-{\alpha }^2\left(k+1\right)g_1-\alpha k{g}_2^{\text{'}}=\text{}\mathrm{0,}\left(k+1\right)g_2^{\text{'}\text{'}}-{\alpha }^2{g}_2+\alpha k{g}_1^{\text{'}}=0,$

которым отвечают разделяющиеся краевые условия:

${\left.\left(k+1\right)\frac{g_2^{\text{'}}}{g_1}+\frac{\left(k-1\right)}{f_2}\left(f_1^{\text{'}}+\frac{f_1}{r}\right)\right|}_{z=0,h}=0{\left.\frac{g_1^{\text{'}}}{g_2}+\frac{f_2^{\text{'}}}{f_1}\right|}_{z=0,h}\text{=0,\hspace{1em}\hspace{1em}}{\left.f_1\right|}_{r=1}\text{=0,\hspace{1em}\hspace{1em}}{\left.f_2^{\text{'}}\right|}_{r=1}\text{=0,}$

или с учетом условий (2.3):

$\left.\left(k+1\right)g_2^{\text{'}}+\alpha \left(k-1\right)g_1\right|=0,{\left.g_1^{\text{'}}-\alpha g_2\right|}_{z=0,h}\text{=0,\hspace{1em}\hspace{1em}}{\left.f_1\right|}_{r=1}\text{=0,\hspace{1em}\hspace{1em}}{\left.f_2^{\text{'}}\right|}_{r=1}\text{=0.}$

Легко видеть, что первые два уравнения (2.5) являются уравнениями Бесселя, соответственно, первого и нулевого порядка; их общие решения могут быть представлены линейными комбинациями функций Бесселя первого и второго рода:

$f_1=C_1{J}_1\left(\alpha r\right)+C_2{Y}_1\left(\alpha r\right)\mathrm{,}{f}_2=C_3{J}_0\left(\alpha r\right)+C_4{Y}_0\left(\alpha r\right)\mathrm{,}$ (2.6)

где $C_1\mathrm{,}{C}_2\mathrm{,}{C}_3\mathrm{,}{C}_4$ - неизвестные константы. Из условия регулярности решения в полюсе константы $C_2\mathrm{,}{C}_4$ следует положить равными нулю. Константы $C_1\mathrm{,}{C}_3$ равны друг другу в силу (2.3) и для краткости могут быть приняты равными единице, так как при разделении переменных отделенные функции определяются с точностью до множителя.

Поскольку параметр $\alpha$ может принимать любое значение из $ℝ\text{}\backslash \text{}\left\{0\right\}$ (случай нулевого $\alpha$ будет рассмотрен отдельно), соотношения (2.6) задают континуальное множество решений, элементы которого в общем случае не удовлетворяют краевым условиям, заданным на боковой поверхности диска. Вместе с тем, если ограничить произвол в выборе параметра $\alpha$только корнями трансцендентного уравнения:

$J_1\left(\alpha \right)=0$ (2.7)

то удается получить счетное семейство решений, удовлетворяющих краевым условиям⁶:

$f_{1n}=J_1\left({\alpha }_{n}r\right)\mathrm{,}{f}_{2n}=J_0\left({\alpha }_{n}r\right)\mathrm{,}n\in \mathbb{N}\mathrm{.}$ (2.8)

Корни трансцендентного уравнения (2.7) функции Бесселя первого порядка, достаточно хорошее асимптотическое приближение которых дается формулой:

${\alpha }_n\approx \frac{{\left(\pi +4n\pi \right)}^2{J}_2\left(\left(1/4+n\right)\pi \right)}{16J_1\left(\left(1/4+n\right)\pi \right)-4\left(\pi +4n\pi \right)J_0\left(\left(1/4+n\right)\pi \right)}\mathrm{.}$

Отдельно рассмотрим случай $\alpha =0$ (решение для этого случая будем помечать нулевым индексом):

$r^2{f}_{10}^{\text{'}\text{'}}+r{f}_{10}^{\text{'}}-f_{10}=0,r^2{f}_{20}^{\text{'}\text{'}}+r{f}_{20}^{\text{'}}=0.$

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

$f_{10}=\frac{C_{5}r}{2}+\frac{C_6}{r}\mathrm{,}{f}_{20}=C_7+C_8\ln r\mathrm{.}$

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

$f_{10}=0,f_{20}=1.$ (2.9)

Учитывая (2.8) и (2.9), вектор-функцию $\left\{f_1\mathrm{,}{f}_2\right\}$ можно представить в виде разложения по системе функций $\left\{\left(\mathrm{0,}1\right)\mathrm{,}\left(J_1\left({\alpha }_{n}r\right)\mathrm{,}{J}_0\left({\alpha }_{n}r\right)\right)\right\}$, и поскольку оператор, порождаемый уравнениями (2.1) и краевыми условиями (2.2), самосопряженный, то эта система функций составляет базис в гильбертовом пространстве двухкомпонентных вектор-функций, заданных над интервалом $\left(\mathrm{0,}1\right)$. Соответственно, будем искать решение исходной неоднородной системы (2.1) в виде разложений:

$u=\sum _n^{\infty }{g}_{1n}{J}_1\left({\alpha }_{n}r\right)\mathrm{,}w=g_{20}+\sum _n^{\infty }{g}_{2n}{J}_0\left({\alpha }_{n}r\right)\mathrm{.}$ (2.10)

Подстановка этих разложений в уравнения (2.1):

$\sum _n^{\infty }\left[g_{1n}^{\text{'}\text{'}}\text{}-{\alpha }^2\left(k+1\right)g_{1n}-\alpha k{g}_{2n}^{\text{'}}\right]J_1\left({\alpha }_{n}r\right)\text{=}\mathrm{0,}$

$\left(k+1\right)g_{20}^{\text{'}\text{'}}+\sum _n^{\infty }\left[\left(k+1\right)g_{2n}^{\text{'}\text{'}}-{\alpha }^2{g}_{2n}+\alpha k{g}_{1n}^{\text{'}}\right]J_0\left({\alpha }_{n}r\right)\text{=}-\mathcal{P}\left(r\right)\mathrm{,}$

и краевые условия (2.2):

${\left.\sum _n^{\infty }\left(g_{1n}^{\text{'}}-\alpha g_{2n}\right)J_1\left({\alpha }_{n}r\right)\right|}_{z=0,h}=\text{}\mathrm{0,}$

${\left.\left(k+1\right)g_{20}^{\text{'}}+\sum _n^{\infty }\left[\left(k+1\right)g_{2n}^{\text{'}}+\alpha \left(k-1\right)g_{1n}\right]J_0\left({\alpha }_{n}r\right)\right|}_{z=0,h}=\text{}0$

приводит к уравнениям относительно коэффициентов разложений, которые сами являются функциями переменной $z$. Для нахождения этих функций подействуем на каждое уравнение системы операторами проектирования:

$P_0\left(\text{f}\right)\underset{0}{\overset{1}{=\int }}\text{f}\mathrm{.}\left(\mathrm{0,}1\right)rdr\mathrm{,}{P}_n\left(\text{f}\right)\underset{0}{\overset{1}{=\int }}\text{f}\mathrm{.}\left(J_1\left({\alpha }_{n}r\right)\mathrm{,}{J}_0\left({\alpha }_{n}r\right)\right)rdr\mathrm{.}$

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

$g_{1n}^{\text{'}\text{'}}-{\alpha }^2\left(k+1\right)g_{1n}-\alpha k{g}_{2n}^{\text{'}}=0,\left(k+1\right)g_{2n}^{\text{'}\text{'}}-{\alpha }^2{g}_{2n}+\alpha k{g}_{1n}^{\text{'}}=-b_n\mathrm{,}$ (2.11)

и одно уравнение для $g_{20}$:

$\left(k+1\right)g_{20}^{\text{'}\text{'}}=-b_0\mathrm{,}$ (2.12)

где $b_0\mathrm{,}{b}_n$ правой части на собственные подпространства, вычисляемые по формулам:

$b_0=2\underset{0}{\overset{1}{\int }}\mathcal{P}\left(r\right)rdr=0,b_n=\frac{2}{J_0{\left({\alpha }_n\right)}^2}\underset{0}{\overset{1}{\int }}\mathcal{P}\left(r\right)r{J}_0\left({\alpha }_{n}r\right)dr=\frac{2pa{J}_1\left({\alpha }_{n}a\right)}{{\alpha }_n\left(1-a^2\right)J_0{\left({\alpha }_n\right)}^2}\mathrm{.}$

Аналогично спроектируем на собственные подпространства краевые условия:

${\left.g_{1n}^{\text{'}}-\alpha g_{2n}\right|}_{z=0,h}\text{=0,\hspace{1em}\hspace{1em}}{\left.\left(k+1\right)g_{2n}^{\text{'}}+\alpha \left(k-1\right)g_{1n}\right|}_{z=0,h}\text{=0,}$ (2.13)

${\left.\left(k+1\right)g_{20}^{\text{'}}\right|}_{z=0,h}=0.\text{}$ (2.14)

Заметим, что сходимость (в среднеквадратичном) разложений к правым частям уравнений (2.1) обеспечивается самосопряженностью задачи Штурма, порождаемой первыми двумя уравнениями (2.5), однородными краевыми условиями при $r=1$ и тем фактом, что каждое собственное значение этой задачи двукратное.

Рассмотрим решение полученных уравнений относительно функций, зависящих от переменной $z$. Уравнение для случая $\alpha =0$ (2.12) может быть сразу проинтегрировано:

$g_{20}=C_9+C_{10}z-\frac{b_0{z}^2}{2\left(k+1\right)}\mathrm{.}$

Из этого соотношения следует, что краевые условия (2.14) для $g_{20}$ не могут быть удовлетворены, если $b_0$отлично от нуля или, что то же самое, если объемные силы не самоуравновешены. Если же это условие выполнено, то $g_{20}$равняется произвольной константе, неопределяемой из краевых условий. Эта константа характеризует смещение цилиндра как жесткого целого вдоль его оси. После построения полного решения исходной неоднородной системы мы сможем выбрать эту константу так, чтобы прогиб в какой-нибудь наперед заданной точке, например, на контуре, был равен нулю.

Осталось только построить решения для систем (2.11). Общее решение для каждой системы может быть построено как сумма общего решения однородной системы:

$g_{1n}^0=\text{}{C}_{1n}{e}^{\alpha z}+C_{2n}{e}^{-\alpha z}-C_{3n}{e}^{\alpha z}\left(k+k\alpha z\right)-C_{4n}{e}^{-\alpha z}\left(k-k\alpha z\right)\mathrm{,}$ (2.15)

$g_{2n}^0=\text{}-C_{1n}{e}^{\alpha z}+C_{2n}{e}^{-\alpha z}-C_{3n}{e}^{\alpha z}\left(2-k\alpha z\right)+C_{4n}{e}^{-\alpha z}\left(2+k\alpha z\right)$

и частного решения неоднородной системы, определяемого методом вариации произвольных постоянных (Лагранжа):

$g_{1n}^{*}=\text{}\frac{-b_n}{2\left(k+1\right)}\underset{0}{\overset{z}{\int }}\left[k\left(z-\zeta \right)sh{\alpha }_n\left(z-\zeta \right)\right]d\zeta =-b_n\frac{k\left({\alpha }_{n}zch{\alpha }_{n}z-sh{\alpha }_{n}z\right)}{2\left(k+1\right){\alpha }_n^2}\mathrm{,}$

$g_{2n}^{*}=\text{}\frac{-b_n}{2{\alpha }_n\left(k+1\right)}\underset{0}{\overset{z}{\int \text{[}}}\left(k+2\right)sh{\alpha }_n\left(z-\zeta \right)-k{\alpha }_n\left(z-\zeta \right)ch{\alpha }_n\left(z-\zeta \right)\text{\hspace{0.05em}]}d\zeta =$(2.16)

$\hspace{1em}=\hspace{1em}{b}_n\frac{k{\alpha }_{n}zsh{\alpha }_{n}z+2\left(1+k\right)\left(1-ch{\alpha }_{n}z\right)}{2\left(k+1\right){\alpha }_n^2}\mathrm{.}$

Константы интегрирования $C_{1n}\mathrm{,}{C}_{2n}\mathrm{,}{C}_{3n}\mathrm{,}{C}_{4n}$ определяются из краевых условий (2.13):

$C_{1n}\text{=}-\frac{K_n}{k+1}\left[e^{2{\alpha }_{n}h}{k}^2-e^{{\alpha }_{n}h}\left({\alpha }_{n}h{k}^2-{\alpha }_{n}hk+k+1\right)-k^2+k+1\right]\mathrm{,}$

$C_{2n}\text{=}\frac{K_n}{k+1}\left[e^{2{\alpha }_{n}h}\left(k^2-k-1\right)-e^{{\alpha }_{n}h}\left({\alpha }_{n}h{k}^2-{\alpha }_{n}hk-k-1\right)-k^2\right]\mathrm{,}$

$C_{3n}\text{=}-\frac{K_n}{k+1}\left[e^{2{\alpha }_{n}h}k-2e^{{\alpha }_{n}h}\left({\alpha }_{n}hk-k-1\right)-3k-2\right]/2$,

$C_{4n}\text{=}\frac{K_n}{k+1}\left[e^{2{\alpha }_{n}h}\left(3k+2\right)-2e^{{\alpha }_{n}h}\left({\alpha }_{n}hk+k+1\right)-k\right]\text{/2,}$

где символом $K_n$обозначено выражение:

$K_n=\frac{-b_n}{2{\alpha }_n^{2}k\left(2e^{{\alpha }_{n}h}{\alpha }_{n}h-e^{2{\alpha }_{n}h}+1\right)}\mathrm{.}$

Подставляя константы интегрирования в общее решение (2.15) и прибавляя частное решение (2.16), получим общее решение для систем неоднородных уравнений (2.11). После подстановки этого решения в соотношения (2.10) получим решение исходной краевой задачи (2.1), (2.2):

$u\text{}\underset{i}{\overset{\infty }{\text{=}\sum }}{K}_n{e}^{-{\alpha }_{n}z}\left\{e^{2{\alpha }_{n}h}\right[k\left({\alpha }_{n}z-1\right)-1]+e^{{\alpha }_{n}h}\left[k\left({\alpha }_{n}h-{\alpha }_{n}z+1\right)+1\right]-$

$-e^{{\alpha }_n\left(h+2z\right)}\left[k\left({\alpha }_{n}h-{\alpha }_{n}z-1\right)-1\right]-e^{2{\alpha }_{n}z}\left[k\left({\alpha }_{n}z+1\right)+1\right]J_1\left({\alpha }_{n}r\right)\mathrm{,}$

$w=\text{}{W}_0+\sum _i^{\infty }{K}_n{e}^{-{\alpha }_{n}z}\{e^{2{\alpha }_{n}h}\left({\alpha }_{n}zk+1\right)+e^{{\alpha }_{n}h}\left[k\left({\alpha }_{n}h-{\alpha }_{n}z\right)-1\right]+e^{{\alpha }_n\left(h+2z\right)}\left[k\left({\alpha }_{n}h-{\alpha }_{n}z\right)+1\right]-$

$-4e^{{\alpha }_n\left(h+z\right)}{\alpha }_{n}hk+2e^{{\alpha }_n\left(2h+z\right)}k-2e^{{\alpha }_{n}z}k+e^{2{\alpha }_{n}z}\left({\alpha }_{n}zk-1\right)J_0\left({\alpha }_{n}r\right)\mathrm{,}$

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

$W_0=-\sum _{i=1}^{\infty }{K}_n\left[\left(e^{2{\alpha }_{n}h}-1\right)\left(2k+1\right)-2e^{{\alpha }_{n}h}{\alpha }_{n}hk\right]J_0\left({\alpha }_n\right)\mathrm{,}$

Напряжения выражаются через перемещения из закона Гука с помощью соотношений Коши:

${\sigma }_{rr}\text{}\text{=}\sum _{i=1}^{\infty }2K_n{e}^{-{\alpha }_{n}z}\left\{\left({\alpha }_n{J}_0\left({\alpha }_{n}r\right)-\frac{J_1\left({\alpha }_{n}r\right)}{r}\right)\right[e^{2{\alpha }_{n}h}\left[k\left({\alpha }_{n}z-1\right)-1\right]+e^{{\alpha }_{n}h}\left[k\left({\alpha }_{n}h-{\alpha }_{n}z+1\right)+1\right]-$

$-e^{{\alpha }_n\left(h+2z\right)}\left[k\left({\alpha }_{n}h-{\alpha }_{n}z-1\right)-1\right]-e^{2{\alpha }_{n}z}\left[k\left({\alpha }_{n}z+1\right)+1\right]-$

$-\left(k-1\right){\alpha }_n{J}_0\left({\alpha }_{n}r\right)\left(e^{2{\alpha }_{n}h}-e^{{\alpha }_{n}h}-e^{{\alpha }_n\left(h+2z\right)}+e^{2{\alpha }_{n}z}\right)\}$

${\sigma }_{\phi \phi }=\text{}\sum _i^{\infty }2K_n{e}^{-{\alpha }_{n}z}\{\left[e^{2{\alpha }_{n}h}k\left({\alpha }_{n}z-1\right)-1\right]+e^{{\alpha }_{n}h}\left[k\left({\alpha }_{n}h-{\alpha }_{n}z+1\right)+1\right]-$

$-e^{{\alpha }_n\left(h+2z\right)}\left[k\left({\alpha }_{n}h-{\alpha }_{n}z-1\right)-1\right]-e^{2{\alpha }_{n}z}\left[k\left({\alpha }_{n}z+1\right)+1\right]]\frac{J_1\left({\alpha }_{n}r\right)}{r}-$

$-\left(k-1\right){\alpha }_n{J}_0\left({\alpha }_{n}r\right)\left(e^{2{\alpha }_{n}h}-e^{{\alpha }_{n}h}-e^{{\alpha }_n\left(h+2z\right)}+e^{2{\alpha }_{n}z}\right)\},$

${\sigma }_{zz}=\sum _{i=1}^{\infty }2K_n{e}^{-{\alpha }_{n}z}[e^{2{\alpha }_{n}h}{\alpha }_{n}z+e^{{\alpha }_{n}h}\left({\alpha }_{n}h-{\alpha }_{n}z\right)-e^{{\alpha }_n\left(h+2z\right)}\left({\alpha }_{n}h-{\alpha }_{n}z\right)-e^{2{\alpha }_{n}z}{\alpha }_{n}z]\left(-k{\alpha }_n{J}_0\left({\alpha }_{n}r\right)\right)\mathrm{,}$

${\sigma }_{rz}=\text{}\sum _{i=1}^{\infty }2K_n{e}^{-{\alpha }_{n}z}[e^{2{\alpha }_{n}h}\left({\alpha }_{n}z-1\right)+e^{{\alpha }_{n}h}\left({\alpha }_{n}h-{\alpha }_{n}z+1\right)+e^{{\alpha }_n\left(h+2z\right)}\left({\alpha }_{n}h-{\alpha }_{n}z-1\right)-$

$-2e^{{\alpha }_n\left(h+z\right)}{\alpha }_{n}h+e^{{\alpha }_n\left(2h+z\right)}-e^{{\alpha }_{n}z}+e^{2{\alpha }_{n}z}\left({\alpha }_{n}z+1\right)]\left(-k{\alpha }_n{J}_1\left({\alpha }_{n}r\right)\right)\mathrm{.}$

Поскольку решение представлено в рядах, то желательно исследовать их сходимость и оценить погрешность частичных сумм в зависимости от числа учитываемых слагаемых. Для этой цели положим $p=1$ и подставим частичные суммы с различным числом слагаемых в уравнения (2.1). Получаемые в результате разложения (которые из общих теоретических рассуждений сходятся в среднеквадратичном) сравним с оригинальной кусочно-постоянной правой частью. Результаты приведены на рис. 2.1. В верхнем ряду показано сопоставление оригинальной кусочно-линейной правой части с результатами действия левой части уравнения (2.1) на частичные суммы 10, 60 и 110 слагаемых, а во втором ряду их отличие от частичной суммы 200 слагаемых. Из рисунка видно, что ряд быстро сходится, но в точке разрыва функции внешних сил имеется неустранимая погрешность (эффект Гиббса). Этот эффект может быть сглажен, если суммирование производить методом средних арифметических (Фейера) [20], однако в этом случае частичные суммы сходятся медленнее, как видно из графиков, приведенных на рис. 2.2.

Рис. 2.1. Оценка сходимости частичных сумм

Fig. 2.1. Estimation of partial sum convergence

Рис. 2.2. Оценка сходимости частичных сумм (суммирование по Фейеру)

Fig. 2.2. Estimation of partial sum convergence (Fejer summation)

3. Модель Кирхгофа-Лява

Перейдем к постановке задачи в рамках теории пластин Кирхгофа-Лява. В ее основе лежат гипотезы, ограничивающие возможные перемещения и распределения напряжений по толщине, а именно [21-23]:

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

Замечание 1. Часто к кинематической гипотезе добавляют дополнительное требование о сохранении длины материального отрезка [24; 25]. Для того чтобы избежать путаницы, кинематическую гипотезу с дополнительным условием, следуя Н.А. Кильчевскому [26], будем называть гипотезой прямых неизменяемых нормальных элементов.

Ограничения на возможные перемещения сужают класс вектор-функций, представляющих пространственное поле перемещений $\text{u}\left(r\mathrm{,}z\right)$, так что оно может быть выражено через одну скалярную функцию прогиба , которая в случае осесимметричной деформации зависит от одной переменной $\omega \left(r\right)$:

$\text{u=}-z{\omega }^{\text{'}}{\text{e}}_r+\omega {\text{e}}_z\mathrm{.}$ (3.1)

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

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

$\left[\sigma \right]=\frac{E}{1+\nu }\left(\begin{array}{ccc}{\epsilon }_{rr}& 0& 0\\ 0& {\epsilon }_{\phi \phi }& 0\\ 0& 0& 0\end{array}\right)+\widehat{E}\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& {\epsilon }_{zz}\end{array}\right)+\frac{E\nu }{1-{\nu }^2}\left({\epsilon }_{rr}+{\epsilon }_{\phi \phi }\right)\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 0\end{array}\right)+q\left(\begin{array}{ccc}0& 0& 1\\ 0& 0& 0\\ 1& 0& 0\end{array}\right)\mathrm{.}$ (3.2)

Здесь ${\epsilon }_{rr}\mathrm{,}{\epsilon }_{\phi \phi }\mathrm{,}{\epsilon }_{zz}$ тензора малых деформаций (1.4), - реакция идеальных связей, сохраняющих прямой угол, определяемая для каждого элементарного объема, $\widehat{E}$- модуль Юнга для растяжения волокон нормальных к срединной поверхности. Таким образом, если материал пластины полагается трансверсально-изотропным, то условие сохранения длины нормали является следствием остальных гипотез.

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

${\epsilon }_{rr}=-z{\omega }^{\text{'}\text{'}}\mathrm{,}{\epsilon }_{\phi \phi }=-\frac{z{\omega }^{\text{'}}}{r}$

и подставим результат в закон состояния (3.2):

${\sigma }_{rr}=-\frac{E}{1-{\nu }^2}z\left({\omega }^{\text{'}\text{'}}+\nu \frac{{\omega }^{\text{'}}}{r}\right)\mathrm{,}{\sigma }_{\phi \phi }=-\frac{E}{1-{\nu }^2}z\left(\frac{{\omega }^{\text{'}}}{r}+\nu {\omega }^{\text{'}\text{'}}\right)\mathrm{,}{\sigma }_{rz}={\sigma }_{zr}=q\mathrm{.}$ (3.3)

Поскольку ${\sigma }_{rz}={\sigma }_{zr}=q$, то из уравнений равновесия:

$\frac{\partial {\sigma }_{rr}}{\partial r}+\frac{{\sigma }_{rr}-{\sigma }_{\phi \phi }}{r}+\frac{\partial {\sigma }_{rz}}{\partial z}=0,\frac{\partial {\sigma }_{zr}}{\partial r}+\frac{{\sigma }_{zr}}{r}+\mathcal{P}\left(r\right)=0$ (3.4)

найдем:

${\sigma }_{rz}={\sigma }_{zr}=q=-\int \left(\frac{\partial {\sigma }_{rr}}{\partial r}+\frac{{\sigma }_{rr}-{\sigma }_{\phi \phi }}{r}\right)dz=\frac{E}{2\left(1-{\nu }^2\right)}{z}^2\frac{\partial }{\partial r}{\nabla }^2\omega +C\left(r\right)\mathrm{,}$ (3.5)

Неизвестная функция $C\left(r\right)$ определяется из условия отсутствия касательных напряжений на лицевых поверхностях:

$q=\frac{E}{2\left(1-{\nu }^2\right)}\left(z^2-\frac{h^2}{4}\right)\frac{\partial }{\partial r}{\nabla }^2\omega \mathrm{.}$

Замечание 2. В настоящем разделе координата $z$ отсчитывается от срединной поверхности пластины (в отличие от раздела 2, в котором осевая координата отсчитывалась от нижнего основания).

Сократив число искомых функций до одной, $\omega$, получаем из полной системы уравнений равновесия (3.4) и краевых условий переопределенную задачу относительно , которая в общем случае решений не имеет. Однако, допуская выполнение краевых условий на лицевой поверхности приближенно (в интегральном смысле) и вводя результирующие усилия (результанты) по следующим формулам:

$T_{rr}\text{=}{\int }_{-h/2}^{h/2}{\sigma }_{rr}dz=0,M_{rr}\text{=}{\int }_{-h/2}^{h/2}z{\sigma }_{rr}dz=-D\left({\omega }^{\text{'}\text{'}}+\nu \frac{{\omega }^{\text{'}}}{r}\right)\mathrm{,}$

$T_{\phi \phi }=\text{}\underset{-h/2}{\overset{h/2}{\int }}{\sigma }_{\phi \phi }dz=0,M_{\phi \phi }\underset{-h/2}{\overset{h/2}{=\int }}z{\sigma }_{\phi \phi }dz=-D\left(\frac{{\omega }^{\text{'}}}{r}+\nu {\omega }^{\text{'}\text{'}}\right)\mathrm{,}$ (3.6)

$Q_{rz}\text{=}{Q}_{zr}=\underset{-h/2}{\overset{h/2}{\int }}{\sigma }_{rz}dz=-D\frac{\partial }{\partial r}{\nabla }^2\omega \mathrm{,}$

где символ $D$обозначает цилиндрическую жесткость:

$D=\frac{E{h}^3}{12\left(1-{\nu }^2\right)}=\frac{k{h}^3}{3\left(k+1\right)}\mathrm{,}$

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

$T_{rr}^{\text{'}}+\frac{T_{rr}-T_{\phi \phi }}{r}=0,Q_{zr}^{\text{'}}+\frac{Q_{zr}}{r}+\mathcal{P}\left(r\right)h=0,M_{rr}^{\text{'}}+\frac{M_{rr}-M_{\phi \phi }}{r}-Q_{rz}=0,$ (3.7)

и краевым условиям, соответствующим жесткому закреплению на круговом контуре:

${\left.\omega \right|}_{r=1}=0,{\left.{\omega }^{\text{'}}\right|}_{r=1}\text{=0.}$ (3.8)

В модели Кирхгофа уравнения равновесия усилий и моментов в радиальном направлении тождественно удовлетворены; осталось только уравнение равновесия перерезывающих сил:

$-Q_{zr}^{\text{'}}-\frac{Q_{zr}}{r}=D{\nabla }^2{\nabla }^2\omega =\mathcal{P}\left(r\right)h\mathrm{,}$ (3.9)

где ${\nabla }^2{\nabla }^2$оператор, который в полярных координатах и в случае осевой симметрии приводится к дифференциальному выражению:

${\nabla }^2{\nabla }^{2}f=f^{\left(4\right)}+\frac{2f^{\text{'}\text{'}\text{'}}}{r}-\frac{f^{\text{'}\text{'}}}{r^2}+\frac{f^{\text{'}}}{r^3}\mathrm{.}$

Располагая решением краевой задачи (3.8), (3.9), по формулам (3.3), (3.5) могут быть определены все компоненты тензора напряжений, а из соотношений (1.5) эффективное напряжение (Мизеса).

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

$\left(\mathrm{1,}{r}^2\mathrm{,}\ln r\mathrm{,}{r}^2\ln r\right)\mathrm{,}$

и используя метод Лагранжа, получим:

$\omega =\text{}\frac{h}{D}\underset{0}{\overset{r}{\int }}\frac{\mathcal{P}\left(\zeta \right)}{8}\text{[2}\left({\zeta }^2+r^2\right)\ln r+\left({\zeta }^2+1\right)\left(1-r^2\right)\text{\hspace{0.05em}]}\zeta d\zeta +$

$+\frac{h}{D}\underset{r}{\overset{1}{\int }}\frac{\mathcal{P}\left(\zeta \right)}{8}\text{[2}\left(r^2+{\zeta }^2\right)\ln \zeta +\left(r^2+1\right)\left(1-{\zeta }^2\right)\text{\hspace{0.05em}]}\zeta d\zeta \mathrm{.}$

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

$\omega =\left\{\begin{array}{cc}\frac{ph}{64D}\left\{r^4+\frac{a^2}{1-a^2}\left[4\left(a^2+2r^2\right)\ln a+\left(3+2r^2\right)\left(1-a^2\right)\right]\right\}& r\le a\\ -\frac{ph{a}^2}{64D\left(1-a^2\right)}\left\{{\left(r^2-1\right)}^2-4\left(a^2+2r^2\right)\ln r-2\left(a^2+2\right)\left(1-r^2\right)\right\}& ra\end{array}\right.$

Заметим, что если устремить параметр $a$к единице, то, используя правило Лопиталя-Бернулли, в пределе можно получить решение для случая однородного нагружения:

$\omega =\underset{a\to 1}{lim}\frac{ph}{64D}\left\{r^4+\frac{a^2}{1-a^2}\left[4\left(a^2+2r^2\right)\ln a+\left(3+2r^2\right)\left(1-a^2\right)\right]\right\}=\frac{ph}{64D}{\left(r^2-1\right)}^2\mathrm{.}$ (3.10)

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

${\sigma }_{rr}\text{=}-\frac{4k}{k+1}z\left({\omega }^{\text{'}\text{'}}+\frac{k-1}{2k}\frac{{\omega }^{\text{'}}}{r}\right)\mathrm{,}{\sigma }_{\phi \phi }=-\frac{4k}{k+1}z\left(\frac{{\omega }^{\text{'}}}{r}+\frac{k-1}{2k}{\omega }^{\text{'}\text{'}}\right)\mathrm{,}$

${\sigma }_{rz}={\sigma }_{zr}=\frac{2k}{k+1}\left(z^2-\frac{h^2}{4}\right)\frac{d}{dr}{\nabla }^2\omega \mathrm{.}$

4. Модель Феппля-фон Кармана

Линейное уравнение изгиба Кирхгофа-Лява не учитывает влияние напряжений в плоскости пластины на ее изгибную жесткость, что оправдано для пластин средней толщины, но ведет к существенным погрешностям (а для ультратонких пластин к ошибкам на несколько порядков) при значительном уменьшении относительной толщины и приближении изгибного НДС к мембранному. Конечно, классическая линейная теория мембран, развитая еще в XVIII веке, приводит к очень простому уравнению Пуассона. Но в него в качестве параметра входит усилие натяжения на контуре, определение которого как раз и представляет основную проблему. С достаточной степенью точности задача определения натяжения решается уравнениями Феппля-фон Кармана, предложенными в начале XX века в работах [29; 30]. Эти уравнения в некотором смысле представляют собой модификацию уравнения Кирхгофа-Лява и плоской теории упругости, к которым добавляются нелинейные термы, содержащие взаимное влияние изгиба и напряжений в плоскости. При этом перемещения представляются разложениями, подобными (3.1):

$\text{u=}\left(\upsilon -z{\omega }^{\text{'}}\right){\text{e}}_r+\omega {\text{e}}_z\mathrm{,}$

где $\upsilon$функция двух координат $r\mathrm{,}\phi$ (в случае осевой симметрии - одной координаты, $r$).

В модели Феппля-фон Кармана допускаются большие повороты и следует различать отсчетную и актуальную формы пластины, поэтому при описании ее деформаций следует использовать один из тензоров конечных деформаций. Для вывода системы уравнений Феппля-фон Кармана используется редуцированный тензор Грина-Сен-Венана ${\text{E}}^{*}$ (деформации фон Кармана):

$\left[{\text{E}}^{*}\right]=\left(\begin{array}{ccc}{\upsilon }^{\text{'}}+\frac{1}{2}{\left({\omega }^{\text{'}}\right)}^2-z{\omega }^{\text{'}\text{'}}& 0& 0\\ 0& \frac{\upsilon }{r}-z\frac{{\omega }^{\text{'}}}{r}& 0\\ 0& 0& 0\end{array}\right)\mathrm{.}$

Компоненты этого тензора связаны с компонентами второго тензора напряжений Пиола-Кирхгофа $\sigma$ уравнениями состояния специального вида [31], которые могут быть получены формальной подстановкой редуцированного тензора Сен-Венана в уравнение состояния (3.2):

${\sigma }_{rr}=\frac{E}{1-{\nu }^2}\left[{\upsilon }^{\text{'}}+\frac{1}{2}{\left({\omega }^{\text{'}}\right)}^2+\nu \frac{\upsilon }{r}-z\left({\omega }^{\text{'}\text{'}}+\nu \frac{{\omega }^{\text{'}}}{r}\right)\right]\mathrm{,}{\sigma }_{\phi \phi }=\frac{E}{1-{\nu }^2}\left[\frac{\upsilon }{r}+\nu \left({\upsilon }^{\text{'}}+\frac{1}{2}{\left({\omega }^{\text{'}}\right)}^2\right)-z\left(\frac{{\omega }^{\text{'}}}{r}+\nu {\omega }^{\text{'}\text{'}}\right)\right]\mathrm{.}$

Касательные напряжения, как и в предыдущей модели, не определяются из уравнений состояния, но могут быть получены из уравнений равновесия, которые могут быть сформулированы в форме (3.4) для компонент первого тензора Пиола-Кирхгофа $\text{P}$ [9]:

$\text{P=F}\cdot \text{}\sigma \mathrm{.}$ (4.1)

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

$\left[\text{F}\right]\approx \left(\begin{array}{ccc}1& 0& -{\omega }^{\text{'}}\\ 0& 1& 0\\ {\omega }^{\text{'}}& 0& 1\end{array}\right)\mathrm{.}$ (4.2)

Подставляя эту матрицу в соотношение выше, получим выражения для компонент первого тензора Пиола-Кирхгофа:

$\left[\text{P}\right]\approx \left(\begin{array}{ccc}{\sigma }_{rr}-{\omega }^{\text{'}}{\sigma }_{zr}& 0& {\sigma }_{rz}\\ 0& {\sigma }_{\phi \phi }& 0\\ {\sigma }_{zr}+{\omega }^{\text{'}}{\sigma }_{rr}& 0& {\omega }^{\text{'}}{\sigma }_{rz}\end{array}\right)\mathrm{.}$

Теперь можно воспользоваться уравнениями равновесия в форме (3.4):

$\frac{\partial {\sigma }_{rr}}{\partial r}+\frac{{\sigma }_{rr}-{\sigma }_{\phi \phi }}{r}+\frac{\partial {\sigma }_{rz}}{\partial z}-\frac{\partial \left({\omega }^{\text{'}}{\sigma }_{zr}\right)}{\partial r}-\frac{\left({\omega }^{\text{'}}{\sigma }_{zr}\right)}{r}\text{=}\mathrm{0,}$

$\frac{\partial {\sigma }_{zr}}{\partial r}+\frac{{\sigma }_{zr}}{r}+{\omega }^{\text{'}}\frac{\partial {\sigma }_{rz}}{\partial z}+\frac{\partial \left({\omega }^{\text{'}}{\sigma }_{rr}\right)}{\partial r}+\frac{\left({\omega }^{\text{'}}{\sigma }_{rr}\right)}{r}+\mathcal{P}\left(r\right)=\text{}0.$

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

$\frac{\partial {\sigma }_{rr}}{\partial r}+\frac{{\sigma }_{rr}-{\sigma }_{\phi \phi }}{r}+\frac{\partial {\sigma }_{rz}}{\partial z}=0,\frac{\partial {\sigma }_{zr}}{\partial r}+\frac{{\sigma }_{zr}}{r}+\frac{\partial \left({\omega }^{\text{'}}{\sigma }_{rr}\right)}{\partial r}+\frac{\left({\omega }^{\text{'}}{\sigma }_{rr}\right)}{r}+\mathcal{P}\left(r\right)=0.$

При сделанных предположениях, касательные напряжения в модели Феппля-фон Кармана полностью совпадают с моделью Кирхгофа:

${\sigma }_{rz}={\sigma }_{zr}=\frac{E}{2\left(1-{\nu }^2\right)}\left(z^2-\frac{h^2}{4}\right)\frac{\partial }{\partial r}{\nabla }^2\omega \mathrm{.}$

Замечание 3. Как видно, для вывода уравнений Феппля-фон Кармана приходится несколько вольно обращаться с нелинейными слагаемыми, что ставит вопрос о корректности этой модели. Этот вопрос был подробно рассмотрен в работах Сьярле [31-34] и др., где он показал, что уравнения Феппля-фон Кармана естественно возникают как первый член асимптотического разложения уравнений нелинейной теории упругости.

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

$\text{}{T}_{rr}^{\text{'}}+\frac{T_{rr}-T_{\phi \phi }}{r}=0,M_{rr}^{\text{'}}+\frac{M_{rr}-M_{\phi \phi }}{r}-Q_{rz}=0,$ (4.3)

$\text{}{Q}_{zr}^{\text{'}}+\frac{Q_{zr}}{r}+{\left({\omega }^{\text{'}}{T}_{rr}\right)}^{\text{'}}+\frac{\left({\omega }^{\text{'}}{T}_{rr}\right)}{r}+\mathcal{P}\left(r\right)h=0.$

Здесь моменты и поперечные силы совпадают с (3.6), а мембранные усилия становятся ненулевыми:

$T_{rr}=\frac{Eh}{1-{\nu }^2}\left[{\upsilon }^{\text{'}}+\frac{1}{2}{\left({\omega }^{\text{'}}\right)}^2+\nu \frac{\upsilon }{r}\right]\mathrm{,}{T}_{\phi \phi }=\frac{Eh}{1-{\nu }^2}\left[\frac{\upsilon }{r}+\nu \left({\upsilon }^{\text{'}}+\frac{1}{2}{\left({\omega }^{\text{'}}\right)}^2\right)\right]\mathrm{.}$ (4.4)

Как и в теории Кирхгофа-Лява, соответствующий выбор функции распределения касательных напряжений позволяет тождественно удовлетворить уравнению равновесия моментов. Уравнение равновесия сил в осевом направлении дополняется нелинейными слагаемыми, которые выражают связь между изгибом пластины и ее растяжением в своей плоскости. Нелинейные слагаемые, в свою очередь, определяются из уравнения равновесия сил в радиальном направлении. Для того чтобы удовлетворить этому уравнению, удобно ввести функцию напряжений Эри $\psi$:

$\frac{{\psi }^{\text{'}}}{r}=\frac{T_{rr}}{h}\mathrm{,}{\psi }^{\text{'}\text{'}}=\frac{T_{\phi \phi }}{h}\mathrm{.}$ (4.5)

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

${\nabla }^2{\nabla }^2\psi +\frac{E}{2}L\left[\omega \mathrm{,}\omega \right]=0,$ (4.6)

где $L$дифференциальный оператор:

$L\left[f\mathrm{,}g\right]=\frac{1}{r}\left(f^{\text{'}}{g}^{\text{'}\text{'}}+g^{\text{'}}{f}^{\text{'}\text{'}}\right)\mathrm{.}$

После подстановки перерезывающей силы и функции Эри в оставшееся уравнение равновесия получим уравнение, замыкающее систему уравнений Феппля-фон Кармана:

$D{\nabla }^2{\nabla }^2\omega -hL\left[\psi \mathrm{,}\omega \right]=\mathcal{P}\left(r\right)h\mathrm{.}$ (4.7)

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

$\upsilon =\frac{1}{E}\left({\psi }^{\text{'}\text{'}}r-\nu {\psi }^{\text{'}}\right)\mathrm{.}$

Более удобной может оказаться формулировка системы Феппля-фон Кармана в перемещениях. Такая формулировка может быть получена путем подстановки соотношений (3.6) и (4.4) в уравнения равновесия (4.3):

$\text{}D{\nabla }^2{\nabla }^2\omega -\frac{Eh}{1-{\nu }^2}\left[\left({\upsilon }^{\text{'}}+\frac{1}{2}{\left({\omega }^{\text{'}}\right)}^2\right)\left({\omega }^{\text{'}\text{'}}+\frac{{\omega }^{\text{'}}}{r}\right)+\left({\upsilon }^{\text{'}\text{'}}+\nu \frac{{\upsilon }^{\text{'}}}{r}\right){\omega }^{\text{'}}+\left({\left({\omega }^{\text{'}}\right)}^2+\nu \frac{\upsilon }{r}\right){\omega }^{\text{'}\text{'}}\right]=\mathcal{P}\left(r\right)h\mathrm{,}$ (4.8)

$\text{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}}{\upsilon }^{\text{'}\text{'}}+\frac{{\upsilon }^{\text{'}}}{r}-\frac{\upsilon }{r^2}+{\omega }^{\text{'}}{\omega }^{\text{'}\text{'}}+\frac{{\left({\omega }^{\text{'}}\right)}^2}{2r}\left(1-\nu \right)=0.$

При этом первое уравнение (4.8) может быть приведено к более простому виду, если выразить ${\upsilon }^{\text{'}\text{'}}$ из второго уравнения и подставить результат в его левую часть:

$D{\nabla }^2{\nabla }^2\omega -\frac{Eh}{1-{\nu }^2}\left[\left({\upsilon }^{\text{'}}+\frac{1}{2}{\left({\omega }^{\text{'}}\right)}^2\right)\left({\omega }^{\text{'}\text{'}}+\nu \frac{{\omega }^{\text{'}}}{r}\right)+\frac{\upsilon }{r}\left(\frac{{\omega }^{\text{'}}}{r}+\nu {\omega }^{\text{'}\text{'}}\right)\right]=\mathcal{P}\left(r\right)h\mathrm{.}$ (4.9)