Splitting of a strip consisting of two identical orthotropic half-strips with isotropy axes symmetrically inclined to the interface

Толық мәтін

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

Задача о трещине, параллельной границам полосы (или полуплоскости, в предельном случае) имеет давнюю историю, начиная с работы Обреимова [1], в которой была посчитана скорость высвобождения энергии (СВЭ) для отслаиваемого от полуплоскости слоя слюды, моделируемого балкой в рамках элементарной теории Бернулли–Эйлера, и сравнением с экспериментальными данными вычислена адгезионная прочность слюды. Задача о полубесконечной трещине параллельной границе упругой изотропной полуплоскости была решена в [2], где были получены точные аналитические выражения для отрывной и сдвиговой мод КИН для трех независимых мод нагружения – продольной и поперечной – силами и изгибающим моментом, прикладываемыми к отделяемому слою вдали от вершины трещины. В работе [3] данное решение было обобщено на случай различных упругих свойств отделяемой полосы и полуплоскости, хотя и при наличии ограничения на сочетание упругих параметров, выражающееся в равенстве нулю второго параметра упругого несоответствия Дундурса [4].

Аналитические решения задачи о центральной полубесконечной трещине в упругой полосе были получены многими авторами. В частности, следует отметить следующие работы. В работе [5] получено решение для случая нормального нагружения берегов трещины парой сил, приложенных на произвольном расстоянии от вершины трещины; однако данное решение нуждается в коррекции. Решение для случая нагружения берегов трещины специально распределенными нормальными усилиями получено в [6]. Наиболее строгое решение для случая нагружения берегов трещины парой нормальных сил получено в [7]. Нагружение парой не только нормальных, но и тангенциальных сил рассмотрено в [8]. В [9] получено решение для случая различных упругих модулей слоев, разделяемых трещиной (хотя и при сохранении равенства нулю второго параметра Дундурса).

Задача о полосе, составленной из двух изотропных слоев произвольной толщины и упругих свойств с полубесконечной интерфейсной трещиной, была рассмотрена в [10] (см. также [11, 12]), где было получено полуаналитическое решение для нагружения осевыми силами и изгибающими моментами, приложенными вдали от вершины трещины. Показано, что произвольная нагрузка данного вида может быть представлена в виде суперпозиции двух базисных мод нагружения, в качестве которых были выбраны нагружение симметрично приложенными моментами и нагружение парой осевых сил и компенсирующим изгибающим моментом, приложенным к одной из отслаиваемых полос.

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

Аналогичные решения были получены для ортотропной полосы с главными осями тензора упругости, совпадающими с геометрическими осями как численными методами [13, 14, 16–19], так и аналитическими [20, 21]. Так в [20] получено решение для случая нагружения парой сил, приложенных к берегам трещины, и таким образом позволяющее получить решение для нагружения произвольной системой нормальных напряжений, симметрично приложенных к берегам трещины. Аналитические выражения для КИН для всех четырех мод нагружения получены в [21]. В работе [22] с помощью решений [3, 9, 21] и процедуры масштабирования были получены решения для составных ортотропных слоев для некоторых комбинаций упругих констант и отношения толщин, заключающихся в удовлетворении дополнительных условий, связывающих отношения упругих констант и толщин составляющих слоев.

Другое возможное обобщение состоит в рассмотрении полосы, состоящей из двух полуполос равной толщины и одного ортотропного материала, но с главными осями тензора упругости симметрично наклоненными к границе раздела. В работе [11] для мод нагружения, соответствующих изгибу моментами и действию продольных сил КИН были получены на основе элементарного балочного решения и соотношения типа ирвиновского для связи КИН и СВЭ [23, 24]. Однако, аналогично изотропному варианту, КИН и СВЭ в данном случае не могут быть получены из элементарных решений для мод нагружения, соответствующих поперечным силам.

Среди работ, в которых рассматриваются межфазные трещины и трещины около границы раздела, следует также отметить такие работы, как [25, 26].

Кроме точных аналитических и численных решений в ряде работ были получены решения с использованием некоторых упрощающих предположений. В ряде случаев такие решения могут быть асимптотически точными для некоторого соотношения параметров. Так существенное упрощение может быть получено при рассмотрении отслоения жесткого тонкого слоя от упругого основания в приближении одномерной теории (изгиба и растяжения брусьев) [27–30], или a priori пренебрегая некоторыми членами в полном наборе уравнений [31–33]. В частности, в [34] было показано, что значения СВЭ полученные на основе теории балок (пластин) асимптотически точны для больших отношений модулей упругости тонкого отслаиваемого слоя и основания.

Методы, используемые для решения задач, включают в себя численные методы, в первую очередь метод конечных элементов (МКЭ) [12–14, 18, 19], численное решение интегральных уравнений [11], использование элементарных решений теории пластин (балок), дающих в ряде случаев точные решения [10, 11, 16], сведение системы к интегральному уравнению (или системе интегральных уравнений) и последующим применением метода Винера–Хопфа [2, 6–9, 21], что позволяет получить выражение для КИН в замкнутой форме. Данный подход будет использован в настоящей работе.

Настоящее исследование посвящено решению плоской задачи теории упругости о полосе, составленной из двух слоев равной толщины из одинакового ортотропного материала с главными осями тензора упругости, наклоненными симметрично относительно границы раздела. Для четырех независимых мод нагружения путем применения двустороннего преобразования Лапласа задача сведена к двум скалярным задачам Винера Хопфа. Получены аналитические выражения для двух мод КИН в терминах элементарных функций и однократных интегралов от комбинаций элементарных функций, зависящие от трех безразмерных параметров: двух комбинаций упругих параметров и угла наклона. Решение может рассматриваться как обобщение решения [21] на случай симметрично наклоненных главных осей анизотропии.

2. Формулировка задачи. Геометрическая конфигурация и система прикладываемых нагрузок. Рассмотрим полосу -h < y < h составленную из двух упругих ортотропных полос с главными осями тензора упругости направленными симметрично относительно границы раздела под углом ö, и центральной полу-бесконечной трещиной y = 0, x < 0 (рис. 1). Здесь введена декартова система координат xy с центром в вершине трещины и осью х, параллельной границам полосы. Механическое поведение системы определяется системой двумерных уравнений теории упругости, включающих

• уравнения равновесия:

$\frac{\partial {\sigma }_{xx}}{\partial x}+\frac{\partial {\sigma }_{xy}}{\partial y}=\mathrm{0,}\frac{\partial {\sigma }_{xy}}{\partial x}+\frac{\partial {\sigma }_{yy}}{\partial y}=\mathrm{0,}$ (2.1)

где ${\sigma }_{xx}(x, y), {\sigma }_{yy}(x, y),{\sigma }_{xy}(x, y)$ – компоненты тензора напряжений;

• соотношения Коши, связывающие компоненты тензора деформации ${\epsilon }_{xx}(x, y), {\epsilon }_{yy}(x, y),{\epsilon }_{xy}(x, y)$ и смещения u(x, y), v(x, y):

${\epsilon }_{xx}=\frac{\partial u}{\partial x},{\epsilon }_{yy}=\frac{\partial v}{\partial y},{\epsilon }_{xy}=\frac{1}{2}\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right);$ (2.2)

• закон Гука, записываемый для анизотропного тела в виде:

$\begin{array}{l}{\epsilon }_{xx}={\beta }_{11}{\sigma }_{xx}+{\beta }_{12}{\sigma }_{yy}+{\beta }_{16}{\sigma }_{xy},\\ {\epsilon }_{yy}={\beta }_{12}{\sigma }_{xx}+{\beta }_{22}{\sigma }_{yy}+{\beta }_{26}{\sigma }_{xy},\\ 2{\epsilon }_{xy}={\beta }_{16}{\sigma }_{xx}+{\beta }_{26}{\sigma }_{yy}+{\beta }_{66}{\sigma }_{xy}.\end{array}$ (2.3)

Здесь ${\beta }_{jk}$ – модифицированные коэффициенты податливости в глобальной системе координат xy, совпадающие для условий плоского напряженного состояния с обычными коэффициентами податливости s_jk, а для условий плоской деформации определяемые как

${\beta }_{jk}^{}=s_{jk}^{}-\frac{s_{j3}^{}{s}_{k3}^{}}{s_{33}^{}}.$ (2.4)

Уравнение совместности деформации получается из (2.2) исключением компонент смещения

$\frac{{\partial }^2}{\partial x^2}{\epsilon }_{yy}+\frac{{\partial }^2}{\partial y^2}{\epsilon }_{xx}=2\frac{{\partial }^2}{\partial x\partial y}{\epsilon }_{xy}.$ (2.5)

Границы y = ±h и y = 0, x < 0 свободны от напряжений:

$\begin{array}{l}{\sigma }_{yy}\left(x,\pm h\right)={\sigma }_{xy}^{}\left(x,\pm h\right)=\mathrm{0,}\left|x\right|<\infty ,\\ {\sigma }_{yy}\left(x,\pm 0\right)={\sigma }_{xy}\left(x,\pm 0\right)=\mathrm{0,}x<0.\end{array}$ (2.6)

 

Рис. 1. Конфигурация и система приложенных нагрузок.

 

Нагрузка в виде трех изгибающих моментов M₁, M₂, M₃, трех продольных P₁, P₂, P₃ и трех поперечных сил V₁, V₂, V₃ (рис. 1) приложена на достаточно большом расстоянии от вершины трещины $(I_1\to \infty ,I_2\to \infty )$так, чтобы считать его бесконечным. Для компенсации моментов, создаваемых поперечными силами, приложены дополнительные моменты $V_1{I}_1, V_2{I}_1, V_3{I}_2$.

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

$P_3=P_1-P_2,M_3=M_1-M_2+\frac{P_1+P_2}{2}h,V_3=V_1-V_2.$ (2.7)

Более того, напряженное состояние вблизи вершины трещины полностью определяется лишь четырьмя независимыми силовыми параметрами [13, 14, 35], так как существуют две комбинации силовых параметров не вызывающие раскрытие трещины. Следуя [9, 21, 22] в качестве четырех силовых параметров выбраны интегральные значения напряжений, действующих на линии продолжения трещины

$\begin{array}{c}M=-\underset{0}{\overset{\infty }{\int }}x{\sigma }_{yy}\left(x\mathrm{,0}\right)dx,V=\underset{0}{\overset{\infty }{\int }}{\sigma }_{yy}\left(x\mathrm{,0}\right)dx,\\ T=\underset{0}{\overset{\infty }{\int }}{\sigma }_{xy}\left(x\mathrm{,0}\right)dx,{\tau }_{\infty }=\underset{x\to \infty }{\lim }{\sigma }_{xy}\left(x\mathrm{,0}\right).\end{array}$ (2.8)

Данные величины отличаются от обычно используемых [10–13], однако связь между различными наборами параметров находится элементарно из условий баланса [21, 35]. Силовые параметры, входящие в (2.8), могут быть выражены через параметры нагружения следующим образом

$M=\frac{M_1+M_2}{2},T=\frac{\left(P_1+P_2\right)}{8}+\frac{3\left(M_2-M_1\right)}{4h},V=\frac{V_1+V_2}{2},{\tau }_{\infty }=\frac{3\left(V_2-V_1\right)}{4h}.$ (2.9)

Здесь первые три соотношения следуют из условий баланса сил и моментов, четвертое соотношение получается из рассмотрения распределения поля напряжений в составной полосе.

Требуется найти поле напряжений вблизи вершины трещины, при условии их интегрируемости

$\left\{{\sigma }_{xy}\right\}=O\left(x^{-n}\right),n<\mathrm{1,}x\to +0.$

3. Коэффициенты податливости. Связь со значениями в собственной системе координат. В системе координат, связанной с главными осями тензора упругости и наклоненной к глобальной системе координат xy на угол $\varphi$(рис. 1), в уравнениях закона Гука остаются четыре константы ${\beta }_{11}^0,{\beta }_{22}^0,{\beta }_{12}^0,{\beta }_{66}^0$ причем решение задачи в напряжениях определяется двумя безразмерными параметрами [17]

$\lambda ={\beta }_{11}^0/{\beta }_{22}^0,\rho =\left({\beta }_{66}^0+2{\beta }_{12}^0\right)/\left(2\sqrt{{\beta }_{11}^0{\beta }_{22}^0}\right),$ (3.1)

на которые наложены термодинамические ограничения $0<\lambda <\infty ,-1<p\le \infty .$

Используя стандартные формулы преобразования [36], коэффициенты податливости в глобальной системе координат xy с учетом записываются в виде

$\begin{array}{c}{\beta }_{11}={\lambda }^{-1}{\beta }_{11}^0\left(\lambda {\cos }^4\varphi +2\sqrt{\lambda }\rho {\cos }^2\varphi {\text{\hspace{0.17em}sin}}^2\varphi +{\text{sin}}^4\varphi \right),\\ {\beta }_{22}={\lambda }^{-1}{\beta }_{11}^0\left({\cos }^4\varphi +2\sqrt{\lambda }\rho {\cos }^2\varphi {\text{\hspace{0.17em}sin}}^2\varphi +\lambda {\sin }^4\varphi \right),\\ {\beta }_{12}={\lambda }^{-1}{8}^{-1}{\beta }_{11}^0\left(1+\lambda +6\sqrt{\lambda }\rho -\left(1+\lambda -2\sqrt{\lambda }\rho \right)\cos 4\varphi \right)-{\beta }_{66}^0/\mathrm{2,}\\ {\beta }_{66}={\lambda }^{-1}{\beta }_{11}^0\left(1+\lambda -2\sqrt{\lambda }\rho \right){\sin }^{2}2\varphi +{\beta }_{66}^0,\\ {\beta }_{16}=-{\lambda }^{-1}{2}^{-1}{\beta }_{11}^0\left(-1+\lambda +\left(1+\lambda -2\sqrt{\lambda }\rho \right)\cos 2\varphi \right)\sin 2\varphi ,\\ {\beta }_{26}={\lambda }^{-1}{2}^{-1}{\beta }_{11}^0\left(1-\lambda +\left(1+\lambda -2\sqrt{\lambda }\rho \right)\cos 2\varphi \right)\sin 2\varphi .\end{array}$ (3.2)

4. Вспомогательная задача об ортотропной полосе. Рассмотрим на плоскости xy верхнюю полуполосу 0 < y < h, | x | < ∞ при следующих граничных условиях:

${\sigma }_{yy}\left(x,h\right)={\sigma }_{xy}\left(x,h\right)=\mathrm{0,}{\sigma }_{xy}\left(x\mathrm{,0}\right)=q_x\left(x\right),{\sigma }_{yy}\left(x\mathrm{,0}\right)=q_y\left(x\right),$ (4.1)

где q_x(x), q_y(y) – заданные функции, и найдем решение упругой задачи при заданных граничных условиях.

Введем функцию напряжений Эри F(x, y), определяемую как

${\sigma }_{xx}=\frac{{\partial }^{2}F}{\partial y^2},{\sigma }_{yy}=\frac{{\partial }^{2}F}{\partial x^2},{\sigma }_{xy}=-\frac{{\partial }^{2}F}{\partial x\partial y}.$ (4.2)

Уравнения при этом тождественно удовлетворяются, а подстановка в (4.2) затем в (2.5) сводит систему уравнений упругости к одному уравнению относительно одного неизвестного – функции напряжений [36]

${\beta }_{22}\frac{{\partial }^{4}F}{\partial x^4}-2{\beta }_{26}\frac{{\partial }^{4}F}{\partial x^3\partial y}+\left(2{\beta }_{12}+{\beta }_{66}\right)\frac{{\partial }^{4}F}{\partial x^2\partial y^2}-2{\beta }_{16}\frac{{\partial }^{4}F}{\partial x\partial y^3}+{\beta }_{11}\frac{{\partial }^{4}F}{\partial y^4}=0.$ (4.3)

Общее решение уравнения (4.3) может быть получено путем применения двустороннего преобразования Лапласа, определяемого как

$\widehat{f}\left(p,y\right)=\underset{-\infty }{\overset{\infty }{\int }}f\left(x,y\right)e^{-px}dx,p\in C.$ (4.4)

Обратное преобразование определяется следующим образом

$f\left(x,y\right)=-\frac{1}{2\pi i}\underset{L}{\int }\widehat{f}\left(p,y\right)e^{px}dp,$ (4.5)

где контур интегрирования L соответствует мнимой оси, и направление интегрирования сверху вниз.

Применение преобразования (4.4) к (4.3) приводит к уравнению для образа

$\begin{array}{c}{\beta }_{22}{p}^4\stackrel{⌢}{F}\left(p,y\right)-2{\beta }_{26}{p}^3\frac{d\stackrel{⌢}{F}\left(p,y\right)}{dy}+\left(2{\beta }_{12}+{\beta }_{66}\right)p^2\frac{d^2\stackrel{⌢}{F}\left(p,y\right)}{d{y}^2}-\\ -2{\beta }_{16}p\frac{d^3\stackrel{⌢}{F}\left(p,y\right)}{d{y}^3}+{\beta }_{11}\frac{d^4\stackrel{⌢}{F}\left(p,y\right)}{d{y}^4}=\mathrm{0,}\end{array}$ (4.6)

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

$\widehat{F}(p,y)=C_1{e}^{k_{1}py}+C_2{e}^{k_{2}py}+C_3{e}^{k_{3}py}+C_4{e}^{k_{4}py},$ (4.7)

где C_i (i = 1÷4) – константы, подлежащие определению из граничных условий, а k_i (i = 1÷4) – корни характеристического уравнения (наиболее просто получить данное решение из решения для случая $\varphi$ = 0, а затем преобразовать корни для произвольного $\varphi$ [36])

$k_{\mathrm{1,2}}=\frac{i{\lambda }^{-\frac{1}{4}}\sqrt{\rho \pm \sqrt{{\rho }^2-1}}\cos \varphi -\sin \varphi }{\cos \varphi +i{\lambda }^{-\frac{1}{4}}\sqrt{\rho \pm \sqrt{{\rho }^2-1}}\sin \varphi },k_{\mathrm{3,4}}=\frac{-i{\lambda }^{-\frac{1}{4}}\sqrt{\rho \pm \sqrt{{\rho }^2-1}}\cos \varphi -\sin \varphi }{\cos \varphi -i{\lambda }^{-\frac{1}{4}}\sqrt{\rho \pm \sqrt{{\rho }^2-1}}\sin \varphi }.$ (4.8)

Отделяя в k_i действительную часть от мнимой, корни можно записать в виде:

$k_{\mathrm{1,3}}=A_1\pm i{B}_1,k_{\mathrm{2,4}}=A_2\pm i{B}_2.$ (4.9)

Для $p$ > 1 коэффициенты A₁, B₁, A₂, B₂ принимают вид

$\begin{array}{l}{A}_1=\frac{{\lambda }^{-\frac{1}{2}}\left(\rho +\sqrt{{\rho }^2\text{}-1}\right)\sin \varphi \cos \varphi -\sin \varphi \cos \varphi }{{\cos }^2\varphi +{\lambda }^{-\frac{1}{2}}\left(\rho +\sqrt{{\rho }^2-1}\right){\sin }^2\varphi },B_1=\frac{{\lambda }^{-\frac{1}{4}}\sqrt{\rho +\sqrt{{\rho }^2-1}}}{{\cos }^2\varphi +{\lambda }^{-\frac{1}{2}}\left(\rho +\sqrt{{\rho }^2-1}\right){\sin }^2\text{}\varphi },\\ A_2=\frac{{\lambda }^{-\frac{1}{2}}\left(\rho -\sqrt{{\rho }^2-1}\right)\sin \varphi \cos \varphi -\sin \varphi \cos \varphi }{{\cos }^2\text{}\varphi +{\lambda }^{-\frac{1}{2}}\left(\rho -\sqrt{{\rho }^2-1}\right){\sin }^2\varphi },B_2=\frac{{\lambda }^{-\frac{1}{4}}\sqrt{\rho -\sqrt{{\rho }^2-1}}}{{\cos }^2\text{}\varphi +{\lambda }^{-\frac{1}{2}}\left(\rho -\sqrt{{\rho }^2\text{}-1}\right){\sin }^2\text{}\varphi }.\end{array}$ (4.10)

Для -1 < $p$ < 1 вид данных коэффициентов достаточно громоздок, однако структура (4.9) сохраняется.

Константы C_i находятся из системы уравнений, определяющих граничные условия (4.1), с учетом (4.2) принимают вид

$\begin{array}{c}{p}^2\stackrel{⌢}{F}\left(p\mathrm{,0}\right)=p^2\left(C_1+C_2+C_3+C_4\right)={\stackrel{⌢}{q}}_y\left(p\right),\\ -p\frac{\partial \stackrel{⌢}{F}\left(p,y\right)}{\partial y}\left|{}_{y=0}\right.=-p^2\left(k_1{C}_1+k_2{C}_2+k_3{C}_3+k_4{C}_4\right)={\stackrel{⌢}{q}}_x\left(p\right),\\ p^2\stackrel{⌢}{F}\left(p,h\right)=p^2\left(C_1{e}^{k_{1}ph}+C_2{e}^{k_{2}ph}+C_3{e}^{k_{3}ph}+C_4{e}^{k_{4}ph}\right)=\mathrm{0,}\\ -p\frac{\partial \stackrel{⌢}{F}\left(p,y\right)}{\partial y}\left|{}_{y=h}\right.=-p^{2}h\left(C_1{k}_1{e}^{k_{1}ph}+C_2{k}_2{e}^{k_{2}ph}+C_3{k}_3{e}^{k_{3}ph}+C_4{k}_4{e}^{k_{4}ph}\right)=0.\end{array}$ (4.11)

Подстановка найденных значений C_i в выражение для образа функции напряжений (4.7) дает решение задачи.

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

$\frac{\partial u}{\partial x},\frac{\partial v}{\partial x}$

выражаются через функцию F(x, y) с использованием (4.2), (2.2), (2.3) как

$\begin{array}{c}{u}^{\text{'}}=\frac{\partial u}{\partial x}={\beta }_{11}\frac{{\partial }^{2}F}{\partial y^2}+{\beta }_{12}\frac{{\partial }^{2}F}{\partial x^2}-{\beta }_{16}\frac{{\partial }^{2}F}{\partial x\partial y},\\ v^{\text{'}}=\frac{\partial v}{\partial x}={\beta }_{26}\frac{{\partial }^{2}F}{\partial x^2}-\left({\beta }_{66}+{\beta }_{12}\right)\frac{{\partial }^{2}F}{\partial x\partial y}+2{\beta }_{16}\frac{{\partial }^{2}F}{\partial y^2}-{\beta }_{11}\int \frac{{\partial }^{3}F}{\partial y^3}dx.\end{array}$ (4.12)

Применение к (4.12) преобразования Лапласа (4.4) после подстановки (3.2), (4.7), (4.9) и коэффициентов C_i, полученных из решения системы (4.11) дает искомую связь

$\begin{array}{l}{\widehat{u}}^{\text{'}}\left(p\mathrm{,0}\right)=C_{11}{\widehat{q}}_x^{}\left(p\right)+C_{12}{\widehat{q}}_y^{}\left(p\right),\\ {\widehat{v}}^{\text{'}}\left(p\mathrm{,0}\right)=C_{21}{\widehat{q}}_x^{}\left(p\right)+C_{22}{\widehat{q}}_y^{}\left(p\right).\end{array}$ (4.13)

Выражения для C₁₁, C₁₂, C₂₁, C₂₂ здесь не выписаны, поскольку они слишком громоздки, при дальнейших вычислениях будут использованы некоторые их комбинации.

5. Сведения задачи к системе уравнений Винера–Хопфа. Рассмотрим преобразование Лапласа (4.4) от следующих величин

$F_{+}\left(p\right)=\left\{\begin{array}{l}{F}_{1+}\left(p\right)\\ F_{2+}\left(p\right)\end{array}\right\}=\underset{0}{\overset{\infty }{\int }}\left\{\begin{array}{l}{\sigma }_{yy}^{}\left(x\mathrm{,0}\right)\\ {\sigma }_{xy}^{}\left(x\mathrm{,0}\right)\end{array}\right\}{e}^{-px}dx,$ (5.1)

$F_{-}\left(p\right)=\left\{\begin{array}{l}{F}_{1-}\left(p\right)\\ F_{2-}\left(p\right)\end{array}\right\}=\underset{-\infty }{\overset{0}{\int }}\frac{\partial }{\partial x}\left\{\begin{array}{l}{v}^{+}\left(x\mathrm{,0}\right)-v^{-}\left(x\mathrm{,0}\right)\\ u^{+}\left(x\mathrm{,0}\right)-u^{-}\left(x\mathrm{,0}\right)\end{array}\right\}{e}^{-px}dx.$ (5.2)

Здесь интервалы интегрирования сокращены, поскольку для x < 0 равны нулю напряжения, а для x > 0 разности смещений, а следовательно и их производные по x. Из равенства нулю подынтегральных выражений в (5.2), (5.1) для отрицательных и положительных аргументов следует голоморфность F₊(p), F_-(p) для правой (Re p > 0) и левой (Re p < 0) полуплоскостей комплексного переменного p, соответственно, включая границу, где данные функции могут иметь полюса [8] (см. также [37]). Верхние индексы “+” и “–” в (5.2) относятся к смещениям для верхнего и нижнего берега трещины. Соотношения между образами производных смещений и напряжениями для верхней полуполосы определяются выражениями (4.13), аналогичные выражения для нижней полуполосы определяются из (4.13) заменой h на -h и ö на -ö. Подстановка данной разности в (5.2) приводит к векторному уравнению

$F_{-}\left(p\right)=K\left(p\right)F_{+}\left(p\right),p\in L$ (5.3)

причем перекрестные члены, определяемые K₁₂, K₂₁ соответствующие влиянию нормальных напряжений на разность тангенциальных смещений и касательных напряжений на разность нормальных смещений сокращаются, что очевидно ввиду симметрии задачи. Здесь контур L соответствует мнимой оси, или может быть трансформирован в соответствие правилам интегрирования по комплексной плоскости. Ненулевые компоненты матрицы K(p) определяются как

$K_{11}\left(p\right)=-\frac{d_1\mathrm{sh}\left(A_1-A_2\right)p+d_2\sin B_{1}p\cos B_{2}p+d_3\cos B_{1}p\sin B_{2}p}{d_5\left(2B_1{B}_2\cos B_{1}p\cos B_{2}p-2B_1{B}_2\mathrm{ch}\left(A_1-A_2\right)p+d_4\sin B_{1}p\sin B_{2}p\right)},$ (5.4)

$K_{22}\left(p\right)=-\frac{d_6\mathrm{sh}\left(A_1-A_2\right)p+d_7\cos B_{2}p\sin B_{1}p+d_8\cos B_{1}p\sin B_{2}p}{d_9\cos B_{1}p\cos B_{2}p-d_9\mathrm{ch}\left(A_1-A_2\right)p+d_{10}\sin B_{1}p\sin B_{2}p},$ (5.5)

где

$\begin{array}{c}{d}_1=4B_1{B}_2\left(A_2(A_1^2-A_1{A}_2+B_1^2)-A_1{B}_2^2\right),\\ d_2=-2B_2\left(A_1^4-2A_1^3{A}_2-2A_1{A}_2{B}_1^2+B_1^2(-A_2^2+B_1^2-B_2^2)+A_1^2(A_2^2+2B_1^2+B_2^2)\right),\\ d_3=-2B_1\left(A_2^2({(A_1-A_2)}^2+B_1^2)-B_2^2(A_1^2+2A_1{A}_2-2A_2^2+B_1^2)+B_2^4\right),\\ d_4={\left(A_1-A_2\right)}^2+B_1^2+B_2^2,d_5=(A_1^2+B_1^2)(A_2^2+B_2^2).\end{array}$ (5.6)

$\begin{array}{c}{d}_6=4\left(A_1-A_2\right)B_1{B}_2,d_7=-2B_2\left({\left(A_1-A_2\right)}^2-B_1^2+B_2^2\right),\\ d_8=-2B_1\left({\left(A_1-A_2\right)}^2+B_1^2-B_2^2\right),d_9=2B_1{B}_2,d_{10}={\left(A_1-A_2\right)}^2+B_1^2+B_2^2.\end{array}$ (5.7)

Кроме того, должны удовлетворяться условия в ключевых точках (нуле и бесконечности)

$F_{+}\left(p\right)=\left\{\begin{array}{c}V+Mp+o\left(p\right)\\ {\tau }_{\infty }{p}^{-1}+T+o\left(1\right)\end{array}\right\},\mathrm{Re}p\to 0+,$ (5.8)

$F_{+}\left(p\right)=O\left(p^{-\nu }\right),\nu <\mathrm{1,}\mathrm{Re}p\to +\infty .$ (5.9)

Первое из условий следует из (2.8), второе – из условия интегрируемости напряжений (2.10).

Уравнения (5.3)–(5.9) определяют матричную задачу Римана (Винера–Хопфа), которая в силу условия K₁₂ = K₂₁ = 0 сводится к двум скалярным задачам: для нормальной и касательной составляющей, в которых положим временно без нарушения общности h = 1.

6. Решение для нормальной составляющей. Ключевым моментом решения является факторизация функции K₁₁(p), т.е. ее представление в виде произведения (отношения) двух функций, голоморфных в левой (Re p > 0) и правой (Re p < 0) полуплоскости комплексного переменного p, за исключением, быть может, точки (0, 0),

$K\left(p\right)={\Lambda }_{-}^{-1}\left(p\right){\Lambda }_{+}\left(p\right).$ (6.1)

По нахождении функций ${\Lambda }_{\pm }^{}$(p), решение находится с помощью обобщенной теоремы Лиувилля (напр., [37])

$F_{\pm }\left(p\right)={\Lambda }_{\pm }^{-1}\left(p\right)\Pi \left(p\right).$ (6.2)

Здесь P(p) – функция, возможно имеющая нули и полюса в нуле и на бесконечности, определяемая из условий (5.8), (5.9) в данных точках.

Используя (4.8), (5.4), уравнение (6.1) может быть представлено в виде

${\Lambda }_{1-}^{-1}\left(p\right){\Lambda }_{1+}\left(p\right)=-{\text{\hspace{0.17em}ctg}}^3\left(p\text{\hspace{0.17em}}\right)g_1\left(\lambda ,\rho ,\varphi \right)G_1\left(p\right),$ (6.3)

$G_1\left(p\right)={\mathrm{tg}}^3\left(p\text{}\right)\text{}{g}_1^{-1}\left(\lambda ,\rho ,\varphi \right)K_{11}\left(p\right),$ (6.4)

$g_1=-\frac{16\sqrt{2}{\lambda }^{1/4}\sqrt{\rho +1}\left({\cos }^2\varphi +\sqrt{\lambda }{\sin }^2\varphi \right)}{3+3\lambda +2\sqrt{\lambda }\rho -4\left(\lambda -1\right)\cos 2\varphi +\left(1+\lambda -2\sqrt{\lambda }\rho \right)\cos 4\varphi }.$ (6.5)

Здесь коэффициенты tg³(p) и g₁$(\lambda ,p,\varphi )$ в (6.3) выбраны так, чтобы функция G₁(p) была голоморфна на всей мнимой оси включая точку (0, 0) и стремилась к единице при p → ±∞.

Для -1 < p < 1 радикалы $\sqrt{p-1}$ в (6.4) и формулах получающихся в дальнейшем становятся чисто мнимыми, однако они появляются в сочетании с другими мнимыми коэффициентами, так что окончательные выражения остаются действительными. Поэтому нет необходимости отдельно рассматривать случай -1 < p < 1.

Использование интегралов Коши для факторизации G₁(p) и стандартное представление для котангенса через гамма-функцию Г(p) (напр., [37]) позволяет осуществить факторизацию в следующем виде

${\Lambda }_{1+}\left(p\right)=\frac{{\Gamma }^3\left(1+p{\pi }^{-1}\right)}{{\Gamma }^3\left(1/2+p{\pi }^{-1}\right)}{J}_{1+}\left(p\right),$ (6.6)

${\Lambda }_{1-}\left(p\right)=g_1\left(\lambda ,\rho ,\varphi \right)\frac{{\Gamma }^3\left(1/2-p{\pi }^{-1}\right)}{{\Gamma }^3\left(-p{\pi }^{-1}\right)}{J}_{1-}\left(p\right),$ (6.7)

$J_{1\pm }\left(p\right)=\exp \left(\frac{1}{2\pi }\underset{-\infty }{\overset{\infty }{\int }}\ln \frac{G_1\left(s\right)ds}{is-p}\right).$ (6.8)

При этом ${\Lambda }_{1+}^{}$(p) в нуле и на бесконечности имеет следующие асимптотические разложения

${\Lambda }_{1+}^{-1}\left(p\right)={\pi }^{3/2}{R}_1\left[1-p{Y}_1\left(\rho \right)\right]+O\left(p^2\right),$ (6.9)

${\Lambda }_{1+}^{-1}\left(p\right)={\left(\frac{\pi }{p}\right)}^{3/2}\text{}+O\left(p^{-1/2}\right),\mathrm{Re}p\to +\infty .$ (6.10)

Здесь

$R_1={\left[\frac{6\sqrt{2}\left(\lambda {\cos }^4\varphi +2\sqrt{\lambda }\rho {\cos }^2\varphi {\sin }^2\varphi +{\sin }^4\varphi \right)}{{\lambda }^{1/4}\sqrt{1+\rho }\left({\cos }^2\varphi +\sqrt{\lambda }{\sin }^2\varphi \right)}\right]}^{1/2}$ (6.11)

квадратный корень из отношения асимптотик G₁(p) в нуле и на бесконечности, остальные величины есть

$\begin{array}{c}{Y}_1\left(\lambda ,\rho ,\varphi \right)=\frac{1}{2\pi }\underset{-\infty }{\overset{\infty }{\int }}\frac{d{L}_1\left(s\right)}{ds}\frac{ds}{s}=\frac{1}{2\pi }\underset{-\infty }{\overset{\infty }{\int }}\left[L_1\left(s\right)-L_1\left(0\right)\right]\frac{ds}{s^2},\\ L_1\left(s\right)=\ln \left[s^3\text{}\frac{d_1\sin \left(A_1-A_2\right)s+d_2\mathrm{ch}{B}_{2}s\mathrm{sh}{B}_{1}s+d_3\mathrm{ch}{B}_{1}s\mathrm{sh}{B}_{2}s}{d_5\left(2\left(\cos \left(A_1-A_2\right)s-\mathrm{ch}{B}_{1}s\mathrm{ch}{B}_{2}s\right)B_1{B}_2+d_4\mathrm{sh}{B}_{1}s\mathrm{sh}{B}_{2}s\right)}\right].\end{array}$ (6.12)

Анализ (6.2), (6.9), (6.10) вместе с (5.8), (5.9) приводит к виду функции П(p)

$\Pi \left(p\right)={\pi }^{-3/2}{R}_1^{-1}\left[V+\left(M+Y_1(\lambda ,\rho ,\varphi )V\right)p\right].$ (6.13)

Появление в (6.13) членов с отрицательны ми степенями p привело бы к нарушению условия (5.9), в то время как появление членов со степенями больше единицы привело бы к нарушению условия (5.8).

Таким образом (6.2), (6.6)–(6.8), (6.13) дают решение поставленной задачи.

Рассмотрим асимптотику p → +∞. Подстановка (6.10), (6.13) в (6.2) дает

$F_{1+}\left(p\right)=R_1^{-1}\left(M+Y_1\left(\lambda ,\rho ,\varphi \right)V\right)p^{-1/2}+o\left(p^{-1/2}\right),\mathrm{Re}p\to +\infty .$ (6.14)

Использование теоремы абелевого типа [37] дает асимптотику напряжения при x → 0+:

${\sigma }_{yy}={\pi }^{-1/2}{R}_1^{-1}\left(M+Y_1(\lambda ,\rho ,\varphi )V\right)x^{-1/2}+o\left(x^{-1/2}\right).$ (6.15)

Откуда, по возвращении к произвольному h, КИН находится как

$K_I{h}^{3/2}=\sqrt{2}{R}_1^{-1}\left(M+Y_1\left(\lambda ,\rho ,\varphi \right)Vh\right).$ (6.16)

Структура выражения (6.16) совпадает со структурой аналогичных выражений [8, 17, 21, 22] для случая ненаклоненных слоев. Для члена, пропорционального моменту, выражение (6.16) совпадает с результатом [10].

7. Решение задачи для сдвиговой моды. Ан алогично предыдущему случаю, используя (4.8), (5.5) аналог выражения записывается следующим образом

${\Lambda }_{2-}^{-1}\left(p\right){\Lambda }_{2+}\left(p\right)=-\text{ctg\hspace{0.17em}}\left(p\text{}\right)\text{\hspace{0.17em}}{g}_2\left(\lambda ,\rho ,\varphi \right)G_2\left(p\right),$ (7.1)

$G_2\left(p\right)=-\text{tg}\left(p\text{}\right)\text{\hspace{0.17em}}{g}_2^{-1}\left(\lambda ,\rho ,\varphi \right)K_{22}\left(p\right),$ (7.2)

$g_2=\frac{2\sqrt{2}{\lambda }^{1/4}\sqrt{1+\rho }\left(\sqrt{\lambda }{\cos }^2\varphi +{\sin }^2\varphi \right)}{\lambda {\cos }^4\varphi +2\sqrt{\lambda }\rho {\cos }^2\varphi {\sin }^2\varphi +{\sin }^4\varphi }.$ (7.3)

Факторизация имеет вид

${\Lambda }_{2+}\left(p\right)=\frac{\Gamma \left(1+p{\pi }^{-1}\right)}{\Gamma \left(1/2+p{\pi }^{-1}\right)}{J}_{2+}\left(p\right),$ (7.4)

${\Lambda }_{2-}\left(p\right)=-\frac{\Gamma \left(1/2-p{\pi }^{-1}\right)}{\Gamma \left(-p{\pi }^{-1}\right)}{g}_1\left(\lambda ,\rho ,\varphi \right)J_{2-}\left(p\right),$ (7.5)

$\begin{array}{c}{J}_{2\pm }\left(p\right)=\\ =\exp \left\{\frac{1}{2\pi }\right.\underset{-\infty }{\overset{\infty }{\int }}\ln \left.\left[s\frac{d_6\sin \left(A_1-A_2\right)t+d_7\mathrm{ch}{B}_{2}t\mathrm{sh}{B}_{1}t+d_8\mathrm{ch}{B}_{1}t\mathrm{sh}{B}_{2}t}{d_9\left(\cos \left(A_1-A_2\right)t-\mathrm{ch}{B}_{1}t\mathrm{ch}{B}_{2}t\right)+d_{10}\mathrm{sh}{B}_{1}t\mathrm{sh}{B}_{2}t}\right]\frac{ds}{is-p}\right\}.\end{array}$ (7.6)

Со следующими асимптотическими представлениями

${\Lambda }_{2+}^{-1}\left(p\right)={\pi }^{1/2}{R}_2\left[1-p{Y}_2\left(\lambda ,\rho ,\varphi \right)\right]+O\left(p^2\right),$ (7.7)

${\Lambda }_{2+}^{-1}\left(p\right)={\left(\frac{\pi }{p}\right)}^{1/2}\text{}+O\left(p^{-1/2}\right),\mathrm{Re}p\to +\infty ,$ (7.8)

$R_2=\frac{2^{3/4}}{{\lambda }^{1/8}{\left(1+\rho \right)}^{1/4}}{\left[\frac{\lambda {\cos }^4\varphi +2\sqrt{\lambda }\rho {\cos }^2\varphi {\sin }^2\varphi +{\sin }^4\varphi }{\sqrt{\lambda }{\cos }^2\varphi +{\sin }^2\varphi }\right]}^{1/2},$ (7.9)

$Y_2\left(\rho \right)=\frac{1}{2\pi }\underset{-\infty }{\overset{\infty }{\int }}\frac{d{L}_2\left(s\right)}{ds}\frac{ds}{s}=\frac{1}{2\pi }\underset{-\infty }{\overset{\infty }{\int }}\left[L_2\left(s\right)-L_2\left(0\right)\right]\frac{ds}{s^2},$ (7.10)

$L_2\left(s\right)=\ln \left[s\frac{d_6\sin \left(A_1-A_2\right)s+d_7\mathrm{ch}{B}_{2}s\mathrm{sh}{B}_{1}s+d_8\mathrm{ch}{B}_{1}s\mathrm{sh}{B}_{2}s}{d_9\left(\cos \left(A_1-A_2\right)s-\mathrm{ch}{B}_{1}s\mathrm{ch}{B}_{2}s\right)+d_{10}\mathrm{sh}{B}_{1}s\mathrm{sh}{B}_{2}s}\right].$ (7.11)

Сравнение (6.2), (7.7), (7.8) вместе с (5.8), (5.9) приводит к

$\Pi \left(p\right)={\pi }^{-3/2}{R}_2^{-1}\left[{\tau }_{\infty }{p}^{-1}+\left(T+Y_2\left(\lambda ,\rho ,\varphi \right){\tau }_{\infty }\right)\right].$ (7.12)

Появление в (7.12) членов степени p^-k, k > 1 привело бы к нарушению условия (5.8), в то время как появление членов с положительными степенями p привело бы к нарушению условия (5.9).

Таким образом, формулы (6.2), (7.4)–(7.6), (7.12) дают решение поставленной задачи.

Рассмотрим асимптотику решения при $p\to +\infty$. Подстановка (7.8), (7.12) в (6.2) дает

$F_{1+}\left(p\right)=R_2^{-1}\left(T+Y_2\left(\lambda ,\rho ,\varphi \right){\tau }_{\infty }\right)p^{-1/2}+o\left(p^{-1/2}\right),\mathrm{Re}p\to +\infty .$ (7.13)

Используя теорему абелевого типа [37] и возвращаясь к произвольному h, получаем асимптотические формулы для напряжений при $x\to 0+$

${\sigma }_{xy}={\pi }^{-1/2}{R}_2^{-1}\left(T+h{Y}_2\left(\lambda ,\rho ,\varphi \right){\tau }_{\infty }\right)h^{-1/2}{x}^{-1/2}+o\left(x^{-1/2}\right).$ (7.14)

Откуда КИН,

$K_{II}{h}^{1/2}=\sqrt{2}{R}_2^{-1}\left(T+Y_2\left(\lambda ,\rho ,\varphi \right){\tau }_{\infty }h\right).$ (7.15)

Структура выражения (7.15) совпадает со структурой (6.16).

8. Результаты численных расчетов. Зависимости КИН от изгибающего момента и пары продольных сил с компенсирующими моментами определяются согласно (6.16), (7.15) алгебраическими функциями от трех параметров $\lambda ,p,\varphi$ и с точностью до обозначений совпадают с выражениями [11], полученными из анализа элементарных балочных решений с использованием результатов [38, 39]. Зависимости КИН от симметрично и несимметрично приложенных поперечных сил не могут быть найдены из элементарных решений и определяются произведением констант (6.11), (7.9) на функции Y₁, Y₂ от тех же трех параметров. Их зависимости от p для $\varphi =\pi /3$ при различных значениях $\lambda$представлены на рис. 2 и 4, зависимости от $\varphi$для различных p и $\lambda$ – на рис. 3 и 5. Отметим, что подобно случаю слоев с главными осями тензоров упругости, расположенных нормально [21], значения функции Y₂ для r близких к -1 становится отрицательным. Аналогично [21], произвольная нагрузка, приложенная достаточно далеко от вершины трещины, может быть представлена как суперпозиция рассмотренных случаев. Существенно, что функции Y₁, Y₂ являются четными функциями от угла наклона слоев $\varphi$, что позволяет избежать ошибок при определении направления его отсчета. Для случая отсутствия наклона слоев $\varphi$= 0 выражения (6.16), (7.15) сводятся к результату [21], однако в отличие от рассмотренного там случая в рассматриваемой задаче оказалось затруднительным выделить некоторую комбинацию параметров в функции Y₁$(\lambda , p, \varphi )$, так что подынтегральное выражение осталось зависящим от всех трех параметров. Тем не мене е, функции Y₁, Y₂ обладают некоторой симметрией: благодаря зависимостям (3.2) можно заключить, что$Y_k(\lambda , p, \varphi )=Y_k(1/\lambda , p, \pi /2-\varphi )$, что видно из рис. 3, 5.

 

Рис. 2. Зависимость функции Y₁ от p для $\varphi =\pi /3$; сплошные линии $-\lambda =1$, пунктирные линии $-\lambda =4$, точечные линии $-\lambda =1/4$.

 

Рис. 3. Зависимость функции Y₁ от $\varphi$: (a) для p = 2, (b) для p = 5, (c) для p = -0.7; сплошные линии – l = 1, пунктирные линии $-\lambda =4$, штрих-пунктирные линии $-\lambda =8$, точечные линии $-\lambda =1/4$.

 

Рис. 4. Зависимость функции Y₂ от p для $\varphi =\pi /3$; сплошные линии $-\lambda =1$, пунктирные линии $-\lambda =4$, точечные линии $-\lambda =1/4$.

 

Рис. 5. Зависимость функции Y₂ от $\varphi$: (a) для p = 2, (b) для p = 5, (c) для p = -0.7; сплошные линии $-\lambda =1$, пунктирные линии $-\lambda =4$, штрих-пунктирные линии $-\lambda =8$, точечные линии $-\lambda =1/4$.

 

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

1) парой изгибающих моментов, симметрично приложенным к отслаиваемым частям полуполос;

2) противоположно направленными продольными силами, приложенными к отслаиваемым частям полуполос;

3) поперечными силами, симметрично приложенными к отслаиваемым частям полуполос;

4) поперечными силами, приложенными в одном направлении к отслаиваемым частям полуполос, и поперечной силой, приложенной к противоположной части составной полосы в противоположном направлении.

Для случаев 2)–4) также прикладываются компенсирующие моменты, так чтобы суммарный момент относительно вершины трещины был равным нулю.

Произвольная нагрузка может быть получена как суперпозиция данных мод нагружения и двух «неактивных» мод, не вызывающих появления КИН.

Ставится задача нахождения КИН.

Путем применения двустороннего преобразования Лапласа задача сведена к матричной задаче Римана, в силу симметрии распадающейся на две скалярные задачи, для которых получены точные аналитические решения. Полученное решение справедливо для всех термодинамически допустимых значений упругих констант. Для четырех независимых активных мод нагружения найдены выражения для коэффициентов интенсивности напряжений в виде комбинаций элементарных функций либо однократных интегралов от комбинаций элементарных функций, зависящих от трех независимых параметров. Значения КИН от изгибающего момента и пары продольных сил с компенсирующими моментами определяются алгебраическими функциями трех параметров $\lambda , p, \varphi$ и с точностью до обозначений совпадают с выражениями [11], полученными из анализа элементарных решений, в то время как зависимости КИН от симметрично и несимметрично приложенных поперечных сил не могут быть найдены из элементарных решений и определяются произведением комбинации констант (6.11), (7.9) на функции Y₁, Y₂, выражающиеся через однократные интегралы от комбинаций алгебраических функций, зависящих от тех же трех параметров $\lambda , p, \varphi$. Аналогично случаю слоев с главными осями тензоров упругости, расположенных нормально [21], значения функции Y₂ для r близких к -1 могут становиться отрицательным. Функции Y₁, Y₂ являются четными функциями от угла наклона слоев$\varphi$. Для случая отсутствия наклона слоев $\varphi$ = 0 выражения (6.16), (7.15), сводятся к результату [21], однако в отличие от рассмотренного там случая, при произвольном угле наклона ö оказалось затруднительным выделить некоторую комбинацию параметров в функции $Y_1(\lambda , p, \varphi )$, так что подынтегральное выражение осталось зависящим от всех трех параметров.

Работа выполнена при финансовой поддержке государственного задания (№ госрегистрации 124012500441-6) для КБУ.

Авторлар туралы

K. Ustinov

A.Yu. Ishlinsky Institute for problem in Mechanics RAS

Хат алмасуға жауапты Автор.
Email: ustinov@ipmnet.ru
Ресей, Moscow

N. Borisova

Federal State Educational Institution of Higher Education “Prince Alexander Nevsky Military University” of the Ministry of Defense of the Russian Federation

Email: nbolo@yandex.ru
Ресей, Moscow

Әдебиет тізімі

Қосымша файлдар