Refined model of thermo-visco-elastic-plastic dynamic deformation of reinforced flexible shallow shells

Full Text

Введение

Изделия из композиционных материалов (КМ) часто подвергаются интенсивному термосиловому нагружению [1–8], при котором материалы фаз композиции могут деформироваться пластически [1–3, 5–7, 9, 10]. Следовательно, актуальна проблема математического моделирования неупругого деформирования анизотропных и композитных тонкостенных элементов конструкций, которая на данный момент времени находится в стадии становления [5–7, 9–24]. Так, в работе [22] была разработана структурная модель неизотермического вязкоупругопластического деформирования многонаправленно армированного материала и проведены соответствующие расчеты неупругой динамики гибких КМ-пластин, по толщине которых температурное поле аппроксимировалось традиционно полиномом второго порядка. Однако в [23, 24] на примерах изгибного динамического поведения армированных пластин и пологих оболочек было показано, что для более точного определения теплового отклика в таких тонкостенных гибких КМ-конструкциях температуру в поперечном направлении нужно аппроксимировать полиномами 6–7-го порядков.

Для моделирования волновых процессов в динамически деформируемых тонкостенных элементах КМ-конструкций и учета их ослабленного сопротивления поперечным сдвигам традиционно используют простейшие неклассические теории Тимошенко–Рейсснера [5, 7, 16, 25–27], Амбарцумяна [22, 24, 28] и Редди–Немировского [29, 30], а также теории более высоких порядков [5, 16, 23, 31, 32], как правило, базирующиеся на гипотезе ломаной линии [5, 16, 31, 32]. В работе [23] было показано, что термоупругопластическую динамику армированных пологих оболочек следует описывать более точными соотношениями изгиба, простейшим вариантом которых является теория Амбарцумяна. Однако в [23] не учитывались вязкие свойства компонентов композиции при их упругом деформировании, поэтому разработанная в [23] теория не позволяет определять в рамках уточненной теории изгиба остаточные прогибы и остаточные деформированные состояния фаз композиции после прекращения осцилляций таких КМ-конструкций.

Нелинейные задачи динамики тонкостенных элементов, как правило, интегрируют численно с использованием явных [5, 22–24] или неявных [6, 33] схем.

На основании вышеизложенного в данной работе моделируется неизотермическая вязкоупругопластическая динамика гибких, многонаправленно армированных пологих оболочек при использовании уточненной теории их изгиба [23] и аппроксимации температуры в поперечном направлении полиномами высоких порядков [23, 24]. Численное решение этой связанной нелинейной термомеханической задачи строится с применением явной пошаговой схемы [5, 22–24].

1. Постановка начально-краевой задачи и метод ее решения

Рассматривается гибкая пологая оболочка толщиной $2h$. Ортогональная система координат $x_{i} $ введена так, что срединная поверхность конструкции — отсчетная ($|x_3|\leqslant h$); координатные линии $x_1 $, $x_2 $ — линии главной кривизны поверхности $x_3 =0$ (рис. 1, на котором искривленность оболочки в силу ее малости не изображена). Панель усилена $N$ семействами непрерывных волокон с плотностями армирования $\omega _k $ ($1\leqslant k\leqslant N$). По толщине конструкции структура перекрестного армирования однородна и в общем случае может быть пространственной [4].

Рис. 1. Элемент пологой оболочки с 2D-структурой армирования
[Figure 1. Element of a shallow shell with a 2D reinforcement structure]

Рис. 2. Взаимная ориентация глобальной и локальной (связанной с волокном $k$-го семейства) систем координат
[Figure 2. The mutual orientation of the global and local (associated with the $k$-th family of fibers) coordinate systems]

С каждым $k$-м семейством арматуры связана локальная ортогональная система координат $x_i^{(k)}$, в которой ось $x_1^{(k)} $ направлена вдоль траектории волокна, а ее ориентация в пространстве $Ox_1 x_2 x_3 $ задается углами сферической системы координат $\theta _k $, $\varphi _k $ (рис. 2). Направляющие косинусы $l_{ij}^{(k)} $ между осями $x_{i}^{(k)}$ и $x_{j}$ ($i, j=\overline{1,3}$, $1\leqslant k\leqslant N$) вычисляются по формулам [23]:
\[ \begin{equation}
\begin{array}{lll}
l_{11}^{(k)} =\sin \theta _k \cos \varphi _k ,& l_{12}^{(k)} =\sin \theta _k \sin \varphi _k , &
l_{13}^{(k)} =\cos \theta _k ,
\\
l_{21}^{(k)} =-\sin \varphi _k , & l_{22}^{(k)} =\cos \varphi _k , & l_{23}^{(k)} =0,
\\
l_{31}^{(k)} =-\cos \theta _k \cos \varphi _k , & l_{32}^{(k)} =-\cos \theta _k \sin \varphi _k , & l_{33}^{(k)} =\sin \theta _k ,\quad 1\leqslant k\leqslant N.
\end{array}
\end{equation} \tag{1} \]

Распределенными касательными напряжениями на лицевых поверхностях КМ-конструкции традиционно пренебрегаем [5–7, 16, 22–32]. В случае пространственной структуры армирования предполагаем: если арматура некоторого $k$-го семейства наклонна ($0<\theta _k <\pi /2$; см. рис. 2 и равенства (1)), то обязательно имеется другое $m$-е семейство наклонных волокон, изготовленных из того же материала, параметры армирования которого таковы: $\theta _{m} =\pi -\theta _k $, $\varphi _{m} =\varphi _k $ и $\omega _{m} =\omega _k $ ($1\leqslant k, m\leqslant N$, $m\ne k$). Структуры с таким пространственным армированием часто реализуются на практике [4]. Плоско-перекрестные структуры армирования (см. рис. 1) автоматически удовлетворяют этому требованию (в этом случае следует формально принять $\omega _{m} =\omega _k \equiv 0$).

При выполнении приведенных выше условий перемещения точек гибкой искривленной панели $U_{i} $ и осредненные деформации ее композиции $\varepsilon _{ij} $ в рамках уточненной теории изгиба аппроксимируются так (геометрическая нелинейность моделируется в приближении Кармана) [23]:
\[ \begin{equation}
\begin{array}{l}
\displaystyle U_{i}(t, \boldsymbol r)=u_{i}(t, \boldsymbol x)-x_3 \partial _{i} w(t, \boldsymbol x)+ {} \\
\displaystyle
\hspace{3cm} {}+ 2\sum _{m=0}^{M}\frac{x_3 ^{m+1} }{h^{2} }\Bigl(\frac{h^{2} }{m+1} -\frac{x_3 ^{2} }{m+3} \Bigr)\bar{\varepsilon }_{i3}^{(m)}(t, \boldsymbol x), \quad i=1, 2,
\\ 
U_3 (t, \boldsymbol r)=w(t, \boldsymbol x),\qquad \boldsymbol x=(x_1 ,x_2 ),\quad \boldsymbol r=(x_1 ,x_2 ,x_3 ); 
\end{array}
\end{equation} \tag{2} \]
\[ \begin{equation}
\begin{array}{l}
\displaystyle
\varepsilon _{ij}(t,\boldsymbol r)=\frac{1}{2}(\partial _{i} u_{j} +\partial _{j} u_{i})-x_3 \partial _{i} \partial _{j} w
+{}
\\
\hspace{3cm}
\displaystyle
{}
+\sum _{m=0}^{M}\frac{x_3 ^{m+1} }{h^{2} } \Big(\frac{h^{2} }{m+1} -\frac{x_3 ^{2} }{m+3} \Big)(\partial _{i} \bar{\varepsilon }_{j3}^{(m)} +\partial _{j} \bar{\varepsilon }_{i3}^{(m)}) +{}
\\
\displaystyle
\hspace{6cm}
{} +\frac{\delta _{ij} w}{R_{i} } +\frac{1}{2} \partial _{i} w\partial _{j} w,\quad i, j= 1, 2,
\\
\displaystyle
\varepsilon _{i3}(t, \boldsymbol r)=\frac{h^{2} -x_3 ^{2} }{h^{2} } \sum _{m=0}^{M}x_3 ^{m} \bar{\varepsilon }_{i3}^{(m)}(t,\boldsymbol x), \quad i= 1, 2, \; \boldsymbol x\in \Omega , \; |x_3 |\leqslant h,\; t\geqslant t_{0} ,
\end{array}
\end{equation} \tag{3} \]
где $u_{i}$ — тангенциальные перемещения точек срединной поверхности ($x_3 =0$) в направлениях $x_{i}$; $w$ — прогиб; $R_{i}$ — главные радиусы кривизны отсчетной поверхности; $t_{0}$ — начальный момент времени $t$; $\delta _{ij}$ — символ Кронекера; $\partial _{i}$ — частная производная по переменной $x_{i}$; $\Omega$ — область, занимаемая конструкцией в плане; $M$ — целое число, задающее количество слагаемых в степенных разложениях по переменной $x_3$. В случае $M=0$ из равенств (2), (3) получаем кинематические соотношения теорий изгиба Амбарцумяна [28] и Редди–Немировского [29, 30]. В выражениях (2) и (3) определению подлежат нестационарные двумерные функции $w$, $u_{i} $ и $\bar{\varepsilon }_{i3}^{(m)} $ ($i=1, 2$, $0\leqslant m\leqslant M$).

В настоящей работе изучается неупругая динамика искривленной КМ-панели как гибкой тонкостенной механической системы, поэтому нормальное напряжение $\sigma_{33}(t, \boldsymbol r)$ с приемлемой для практических приложений точностью по толщине оболочки можно аппроксимировать формулой [23, 27]:
\[ \begin{multline}
\sigma _{33}(t, \boldsymbol r)=\frac{\sigma _{33}^{(+)}(t, \boldsymbol x)-\sigma _{33}^{(-)}(t,\boldsymbol x)}{2h} x_3 +{} \\
{}+\frac{\sigma _{33}^{(+)}(t, \boldsymbol x)+\sigma _{33}^{(-)}(t, \boldsymbol x)}{2} ,
\quad
\boldsymbol x\in \Omega ,\; |x_3 |\leqslant h,\; t\geqslant t_{0} ,
\end{multline} \tag{4} \]
где $\sigma _{33}^{(\pm )}(t, \boldsymbol x)\equiv \sigma _{33}(t, \boldsymbol x,\pm h)$ — известные из силовых граничных условий напряжения на нижней ($-$) и верхней ($+$) лицевых поверхностях оболочки.

Как уже отмечалось во введении, решение рассматриваемой задачи будем строить, используя явную численную схему [22–24]. Согласно этому искомые функции будем вычислять в дискретные моменты времени $t_{n+1} =t_{n} +\Delta $ ($n=0, 1, 2,\dots$), где $\Delta ={\rm const}>0$ — шаг по времени. Как и в [22], предполагаем, что при $t=t_{n-1}, t_{n} $ уже определены следующие функции:
\[ \begin{equation}
\begin{array}{ll}
\mathop{\sigma _{ij}^{(k)} }\limits^{m}(\boldsymbol r)\equiv \sigma _{ij}^{(k)}(t_{m}, \boldsymbol r), &
\mathop{\dot{\sigma }_{ij}^{(k)} }\limits^{n-1}(\boldsymbol r)\equiv \dot{\sigma }_{ij}^{(k)} (t_{n-1} , \boldsymbol r),
\\
\mathop{\Theta }\limits^{m}(\boldsymbol x)\equiv \Theta(t_{m} , \boldsymbol x), &
\mathop{\dot{\Theta }}\limits^{n-1}(\boldsymbol x)\equiv \dot{\Theta }(t_{n-1},\boldsymbol x), \\
i,j=\overline{1,3},\; m=n-1,n, & 0\leqslant k\leqslant N,\; \boldsymbol x\in \Omega ,\; |x_3 |\leqslant h,
\end{array}
\end{equation} \tag{5} \]
где $\sigma _{ij}^{(k)} $ — напряжения в $k$-м компоненте композиции ($k=0$ — связующее, $k\geqslant 1$ — арматура $k$-го семейства); $\Theta $ — температура КМ-конструкции; точка сверху — частная производная по времени $t$. Тогда в текущий момент времени $t_{n} $ определяющие соотношения для $k$-го материала композиции, связывающие скорости его деформаций $\dot{\varepsilon }_{ij}^{(k)} $, скорости напряжений $\dot{\sigma }_{ij}^{(k)} $ и температуру $\Theta $, при учете обозначений, аналогичных (5), можно записать в следующей матричной форме [22]:
\[ \begin{equation}
\mathop{\dot{{\boldsymbol \sigma }}_k }\limits^{\!\!\!n} =\mathop{{\boldsymbol B} _k }\limits^{\!\!\!n} \mathop{\dot{{\boldsymbol \varepsilon }}_k }\limits^{\!\!\!n} +\mathop{\boldsymbol p_k }\limits^{\!\!\!n} ,\quad 0\leqslant k\leqslant N,
\end{equation} \tag{6} \]
где
\[ \begin{equation}
\begin{array}{l}
\boldsymbol \sigma_k =(\sigma _1 ^{(k)} \, \sigma _2 ^{(k)} \, \sigma _3 ^{(k)} \, \sigma _4 ^{(k)} \, \sigma _5 ^{(k)} \, \sigma _6 ^{(k)})^{\top}
\equiv (\sigma _{11}^{(k)} \, \sigma _{22}^{(k)} \, \sigma _{33}^{(k)} \, \sigma _{23}^{(k)} \, \sigma _{31}^{(k)} \, \sigma _{12}^{(k)} )^{\top} ,
\\
\boldsymbol \varepsilon _k =(\varepsilon _1 ^{(k)} \, \varepsilon _2 ^{(k)} \, \varepsilon _3 ^{(k)} \, \varepsilon _4 ^{(k)} \, \varepsilon _5 ^{(k)} \, \varepsilon _6 ^{(k)})^{\top}
\equiv (\varepsilon _{11}^{(k)} \, \varepsilon _{22}^{(k)} \, \varepsilon _{33}^{(k)} \, 2\varepsilon _{23}^{(k)} \, 2\varepsilon _{31}^{(k)} \, 2\varepsilon _{12}^{(k)} )^{\top} ;
\end{array}
\end{equation} \tag{7} \]
\[ \begin{equation}
\begin{array}{l}
\mathop{\boldsymbol B _k }\limits^{\!\!\!n} \equiv \mathop{\bar{\boldsymbol V}_k ^{-1} }\limits^{\!\!\!n}
\bigl(
\mathop{\bar{\boldsymbol Z}_k }\limits^{n} -
\mathop{\bar{\bar{\boldsymbol Z}}_k }\limits^{n} \bigr), \\
\mathop{\boldsymbol p_k }\limits^{n} \equiv \mathop{\bar{\boldsymbol V}_k ^{-1} }\limits^{n}
\Bigl[
\mathop{\boldsymbol V_k }\limits^{n} \mathop{\boldsymbol \sigma _k }\limits^{n-1/2} +
\dfrac{2}{\Delta } \bigl(\mathop{\Theta }\limits^{n} - \mathop{\Theta }\limits^{n-1/2} \bigr)
\mathop{\boldsymbol \beta _k }\limits^{n} \Bigr], \\
\mathop{\bar{\boldsymbol V}_k }\limits^{n} \equiv \boldsymbol I-\frac{\Delta }{2} \mathop{\boldsymbol V_k }\limits^{n} ;
\end{array}
\end{equation} \tag{8} \]
\[ \begin{equation}
\mathop{{\boldsymbol \sigma }_k }\limits^{n-1/2} =\mathop{{\boldsymbol \sigma }_k }\limits^{n-1} +\frac{\Delta }{2} \mathop{\dot{{\boldsymbol \sigma }}_k }\limits^{n-1} ,\quad \mathop{\Theta }\limits^{n-1/2} \equiv \mathop{\Theta }\limits^{n-1} +\frac{\Delta }{2} \mathop{\dot{\Theta }}\limits^{n-1} ,
\quad
0\leqslant k\leqslant N;
\end{equation} \tag{9} \]
$\boldsymbol I$ — единичная $6{\times}6$-матрица; $\bar{\boldsymbol V}_k ^{-1} $ — матрица, обратная $\bar{\boldsymbol V}_k $; $\bar{\boldsymbol Z}_k =\bigl(\bar{z}_{ij}^{(k)} \bigr)$, $\bar{\bar{\boldsymbol Z}}_k  =\bigl(\bar{\bar{z}}_{ij}^{(k)} \bigr)$ и $\boldsymbol V_k =\bigl(v_{ij}^{(k)} \bigr)$ — симметричные матрицы; ${\boldsymbol \beta }_k =\bigl(\beta _{i}^{(k)} \bigr)$ ($i,j=\overline{1,6}$) — вектор-столбец, ненулевые элементы которых определяются так:
\[ \begin{equation*}
\bar{z}_{ij}^{(k)} =2\delta _{ij} G_k +\bar{\lambda }_k ,
\quad
\bar{z}_{ll}^{(k)} =G_k ,
\end{equation*} \]
\[ \begin{equation*}
v_{ij}^{(k)} =\frac{1}{3} \Bigl[\frac{G_k }{\eta _k } \Bigl(1-\frac{\gamma _k }{1+g_k } \Bigr)-\frac{K_k }{\mu _k } \Bigr] -\delta _{ij} \frac{G_k }{\eta _k } \Bigl(1-\frac{\gamma _k }{1+g_k } \Bigr),
\end{equation*} \]
\[ \begin{equation*}
v_{ll}^{(k)} =-\frac{G_k }{\eta _k } \Big(1-\frac{\gamma _k }{1+g_k } \Big),
\end{equation*} \]
\[ \begin{equation*}
\beta _{i}^{(k)} =\frac{\gamma _k \tau _{\Theta }^{(k)} s_{i}^{(k)} }{\tau _{s}^{(k)} (1+g_k )} -3\alpha _k K_k ,
\quad
\beta _{l}^{(k)} =\frac{\gamma ^{(k)} \tau _{\Theta }^{(k)} s_{l}^{(k)} }{\tau _{s}^{(k)}(1+g_k )}, \quad i,j=\overline{1,3},\; l=\overline{4,6};
\end{equation*} \]
\[ \begin{equation*}
\bar{\bar{z}}_{ij}^{(k)} =\frac{\gamma _k G_k s_{i}^{(k)} s_{j}^{(k)} }{\tau _{s}^{(k)2} (1+g_k )}, \quad i,j=\overline{1,6};
\end{equation*} \]
\[ \begin{equation}
G_k =\frac{E_k }{2(1+\nu _k )} ,\quad K_k =\frac{E_k }{3(1-2\nu _k )} ,
\end{equation} \tag{10} \]
\[ \begin{equation*}
\bar{\lambda }_k =\frac{\nu _k E_k }{(1+\nu _k )(1-2\nu _k )} ,\quad g_k =\frac{\bar{G}_k }{G_k } ,
\end{equation*} \]
\[ \begin{equation*}
\gamma _k =
\begin{cases}
0 & \text{при } T_k <\tau _{s}^{(k)}(\chi _k ,\Theta) \text{ или } T_k =\tau _{s}^{(k)}(\chi _k ,\Theta),
W_k \leqslant 0,
\\
1 & \text{при } T_k =\tau _{s}^{(k)}(\chi _k ,\Theta), W_k >0,
\end{cases}
\end{equation*} \]
\[ \begin{equation*}
W_k \equiv \boldsymbol s_k ^{\top} \dot{{\boldsymbol \varepsilon }}_k -
\frac{\tau _{s}^{(k)2} (\chi _k ,\Theta )}{\eta _k } -
\frac{\tau _{s}^{(k)} (\chi _k ,\Theta )\tau _{\Theta }^{(k)} (\chi _k ,\Theta )}{G_k } \dot{\Theta },\quad \tau _{\Theta }^{(k)} \equiv \frac{\partial \tau _{s}^{(k)} }{\partial \Theta } ,
\end{equation*} \]
\[ \begin{equation*}
\bar{G}_k \equiv \frac{\partial \tau _{s}^{(k)} }{\partial \chi _k } ,\quad T_k ^{2} =\frac{1}{2} \sum _{i=1}^{3}(s_{i}^{(k)})^{2} +\sum _{i=4}^{6}(s_{i}^{(k)})^{2} ,\quad \chi _k =\int _{t_{0} }^{t}\sqrt{2\sum _{i,j=1}^{3}\dot{p}_{ij}^{(k)} \dot{p}_{ij}^{(k)} } dt ,
\end{equation*} \]
\[ \begin{equation*}
\dot{p}_{ij}^{(k)} =\frac{s_{ij}^{(k)} W_k }{\tau _{s}^{(k)2}(1+g_k )} ,\quad s_{ij}^{(k)} =\sigma _{ij}^{(k)} -\frac{\delta _{ij} }{3} \sum _{m=1}^{3}\sigma _{mm}^{(k)}, \quad i, j=\overline{1, 3},
\end{equation*} \]
\[ \begin{equation*}
\boldsymbol s_k =\bigl(s_1 ^{(k)} \, s_2 ^{(k)} \, s_3 ^{(k)} \, s_4 ^{(k)} \, s_5 ^{(k)} \, s_6 ^{(k)} \bigr)^{\top} \equiv
\bigl(s_{11}^{(k)} \, s_{22}^{(k)} \, s_{33}^{(k)} \, s_{23}^{(k)} \, s_{31}^{(k)} \, s_{12}^{(k)} \bigr)^{\top} ;
\end{equation*} \]
$E_k =E_k (\Theta)$, $\nu _k =\nu _k (\Theta)$ — модуль Юнга и коэффициент Пуассона материала $k$-й фазы композиции; $\eta _k =\eta _k (\Theta)$, $\mu _k =\mu _k (\Theta)$ — коэффициенты линейной сдвиговой и объемной вязкости (предполагается: вязкоупругое поведение компонентов композиции описывается соотношениями модели Максвелла–Больцмана; см. (2) в [22]); $\alpha _k =\alpha _k (\Theta)$ — коэффициент линейного температурного расширения; $\tau _{s}^{(k)}(\chi _k ,\Theta)$ — предел текучести при чистом сдвиге, который при активном пластическом деформировании равен интенсивности касательных напряжений $T_k $ и зависит от параметра упрочнения $\chi _k$ и температуры $\Theta $; $\dot{p}_{ij}^{(k)}$ — скорости пластических деформаций; $\gamma _k$ — параметр переключения, задающий при $\gamma _k =0$ термовязкоупругое деформирование, нейтральное нагружение или разгрузку, а при $\gamma _k =1$ — активное нагружение $k$-го материала композиции при его пластическом деформировании; индекс $\top$ — операция транспонирования. В соотношениях (10) по повторяющемуся индексу $l$ суммирования нет; в выражении для $W_k$ скорость температуры $\dot{\Theta }$ при $t=t_{n}$ следует заменить ее конечно-разностным аналогом, полученным на основе формулы трапеций [22]:
\[ \begin{equation}
\mathop{\dot{\Theta }}\limits^{n} =\frac{2}{\Delta } \bigl(\mathop{\Theta }\limits^{n} -\mathop{\Theta }\limits^{n-1/2} \bigr),\quad n=1, 2, 3,\dots
\end{equation} \tag{11} \]

Согласно формулам (9) и предположениям (5), в момент времени $t_{n} $ в соотношениях (8) и (11) уже известны значения функций $\mathop{{\boldsymbol \sigma }_k }\limits^{n-1/2} $ и $\mathop{\Theta }\limits^{n-1/2}$.

Как видно из выражений (10), элементы матриц $\bar{\bar{\boldsymbol Z}}_k$, $\boldsymbol V_k $, $\bar{\boldsymbol V}_k $ и вектора ${\boldsymbol \beta }_k $ (а в случае термочувствительности материалов композиции и матрицы $\bar{\boldsymbol Z}_k $) зависят от решения задачи, поэтому определяющие соотношения (6) при учете (8) и (10) являются нелинейными. Линеаризуем равенства (6), применяя итерационный процесс, который аналогичен методу переменных параметров упругости [34]. Используя линеаризованные соотношения (6) при учете (8)–(11) и повторяя рассуждения из [22], в рассматриваемый момент времени $t_{n}$ на текущей итерации получим линейное определяющее уравнение, записанное в матричной форме и характеризующее неизотермическое вязкоупругопластическое поведение КМ:
\[ \begin{equation}
\mathop{\dot{{\boldsymbol \sigma }}}\limits^{ n} =\mathop{{\boldsymbol B} }\limits^{ n} \mathop{\dot{{\boldsymbol \varepsilon }}}\limits^{ n} +\mathop{\boldsymbol p}\limits^{n} ,\quad n=0, 1, 2, \dots,
\end{equation} \tag{12} \]

где
\[ \begin{equation*}
\boldsymbol B\equiv \biggl(\omega _{0} \boldsymbol B_{0}+\sum _{k=1}^{N}\omega _k \boldsymbol B_k \boldsymbol E_k \biggr) \boldsymbol H^{-1} ,\;
\boldsymbol p\equiv \boldsymbol f- \boldsymbol B \boldsymbol g,\;
\boldsymbol f\equiv \omega _{0} \boldsymbol p_{0} +\sum _{k=1}^{N}\omega _k (\boldsymbol p_k + \boldsymbol B_k \boldsymbol r_k ),
\end{equation*} \]
\[ \begin{equation}
\boldsymbol H\equiv \omega _{0} \boldsymbol I + \sum _{k=1}^{N}\omega _k \boldsymbol E_k ,\quad
\boldsymbol g\equiv \sum _{k=1}^{N}\omega _k \boldsymbol r_k ,\quad
\omega _{0} \equiv 1-\sum _{k=1}^{N}\omega _k ,\quad
\boldsymbol r_k \equiv \boldsymbol D_k ^{-1} \boldsymbol \varsigma _k ,\qquad
\end{equation} \tag{13} \]
\[ \begin{equation*}
\boldsymbol E_k \equiv \boldsymbol D_k ^{-1} \boldsymbol C_k , \quad 1\leqslant k\leqslant N;
\end{equation*} \]
$\dot{{\boldsymbol \sigma }}$, $\dot{{\boldsymbol \varepsilon }}$ — шестикомпонентные векторы-столбцы скоростей усредненных напряжений $\dot{\sigma }_{ij} $ и деформаций $\dot{\varepsilon }_{ij} $ в композиции, по структуре аналогичные (7); $\boldsymbol B$, $\boldsymbol E_k $ и $\boldsymbol C_k $ — $6{\times 6}$-матрицы; $\boldsymbol D_k ^{-1} $ и $\boldsymbol H^{-1} $ — матрицы, обратные $6{\times} 6$-матрицам $\boldsymbol D_k $ и $\boldsymbol H$; $\boldsymbol p$, $\boldsymbol f$, $\boldsymbol g$, $\boldsymbol r_k $ и $\boldsymbol \varsigma _k $ — шестикомпонентные векторы-столбцы; $\omega _{0} $ — относительное объемное содержание связующего в репрезентативной ячейке КМ. Элементы вектора $\boldsymbol \varsigma _k =(\varsigma _{i}^{(k)})$ и матриц $\boldsymbol C_k =(c_{ij}^{(k)})$, $\boldsymbol D_k =(d_{ij}^{(k)} )$ вычисляются по формулам
\[ \begin{equation}
\begin{array}{l}
\displaystyle
c_{1j}^{(k)} =d_{1j}^{(k)} =q_{1j}^{(k)} ,\quad c_{ij}^{(k)} =\sum _{l=1}^{6}g_{il}^{(k)} b_{lj}^{(0)} ,\quad d_{ij}^{(k)} =\sum _{l=1}^{6}g_{il}^{(k)} b_{lj}^{(k)} ,
\\
\displaystyle
\varsigma _1 ^{(k)} =0,\quad \varsigma _{i}^{(k)} =\sum _{l=1}^{6}g_{il}^{(k)} (p_{l}^{(0)} -p_{l}^{(k)}) ,\quad i=\overline{2,6},\; j=\overline{1,6},\; 1\leqslant k\leqslant N,
\end{array}
\end{equation} \tag{14} \]
где
\[ \begin{equation}
\begin{array}{l}
g_{11}^{(k)} =q_{11}^{(k)} =l_{11}^{(k)} l_{11}^{(k)} ,\quad g_{12}^{(k)} =q_{12}^{(k)} =l_{12}^{(k)} l_{12}^{(k)} ,\quad \dots,
\\
g_{16}^{(k)} =2q_{16}^{(k)} =2l_{12}^{(k)} l_{11}^{(k)} ,\quad \dots,
\\
2g_{61}^{(k)} =q_{61}^{(k)} =2l_{21}^{(k)} l_{11}^{(k)} ,\quad \dots,
\\
g_{66}^{(k)} =q_{66}^{(k)} =l_{11}^{(k)} l_{22}^{(k)} +l_{12}^{(k)} l_{21}^{(k)} ,\quad 1\leqslant k\leqslant N;
\end{array}
\end{equation} \tag{15} \]
$b_{lj}^{(k)} $ и $p_{l}^{(k)} $ ($l,j=\overline{1,6}$) — линеаризованные элементы матриц $\boldsymbol B_k $ и векторов $ \boldsymbol p_k $ ($0\leqslant k\leqslant N$) в равенствах (6); направляющие косинусы $l_{ij}^{(k)} $ ($i,j=\overline{1,3}$) определены в (1). Не выписанные в явном виде в (15) элементы $6{\times} 6$-матриц $\boldsymbol G_k =(g_{ij}^{(k)})$, $\boldsymbol Q_k =(q_{ij}^{(k)})$ приведены в табл. (21.40), (21.44) в [26]. Эти матрицы соответственно определяют преобразования векторов $\dot{{\boldsymbol \sigma }}_k $ и $\dot{{\boldsymbol \varepsilon }}_k $ (см. \eqref{yan:eq6}, (7)) при переходе от глобальной системы $x_{j} $ к локальной $x_{i}^{(k)} $ (см. рис. 2). Попутно при выводе соотношений (12) и (13) получаются матричные зависимости
\[ \begin{equation}
\dot{{\boldsymbol \varepsilon }}_{0} = \boldsymbol H^{-1} \dot{{\boldsymbol \varepsilon }}- \boldsymbol H^{-1} g,
\quad
\dot{{\boldsymbol \varepsilon }}_k = \boldsymbol E_k \dot{{\boldsymbol \varepsilon }}_{0} + \boldsymbol r_k , \quad 1\leqslant k\leqslant N,
\end{equation} \tag{16} \]
которые позволяют последовательно определить скорости деформаций связующего материала $\dot{{\boldsymbol \varepsilon }}_{0} $ и армирующих волокон $k$-го семейства $\dot{{\boldsymbol \varepsilon }}_k $ через скорости деформаций КМ $\dot{{\boldsymbol \varepsilon }}$. Используя же (6) и (16), можно через $\dot{{\boldsymbol \varepsilon }}$ выразить и скорости напряжений в компонентах композиции $\dot{{\boldsymbol \sigma }}_k $ ($0\leqslant k\leqslant N$).

Согласно структуре вектор-столбцов $\dot{{\boldsymbol \sigma }}$ и $\dot{{\boldsymbol \varepsilon }}$, из третьего уравнения линеаризованной системы (12) можно определить скорость линейной поперечной деформации:
\[ \begin{equation}
\dot{\varepsilon }_{33} =b_{33}^{-1}(\dot{\sigma }_{33} -p_3 -b_{31} \dot{\varepsilon }_{11} -b_{32} \dot{\varepsilon }_{22} -2b_{34} \dot{\varepsilon }_{23} -2b_{35} \dot{\varepsilon }_{31} -2b_{36} \dot{\varepsilon }_{12}).
\end{equation} \tag{17} \]
Здесь $b_{3i} $ ($i=\overline{1,6}$) и $p_3 $ — элементы матрицы $\boldsymbol B$ и вектора $\boldsymbol p$ в определяющем соотношении (12). Значение скорости $\dot{\sigma }_{33} $ в (17) определяется за счет дифференцирования по времени аппроксимации (4), а скорости деформаций $\dot{\varepsilon }_{ij} $ в правой части (17) — путем дифференцирования по $t$ кинематических соотношений (3), т. е. в конечном итоге выражаются через искомые функции $w$, $\dot{w}$, $\dot{u}_{i} $ и $\dot{\bar{\varepsilon }}_{i3}^{(m)} $, $i=1, 2$, $0\leqslant m\leqslant M$. (В соотношениях (10), (13), (14), (16) и (17) опущен верхний индекс $n$, означающий текущий дискретный момент времени $t_{n} $.)

Как и в [22–24], температуру в поперечном направлении оболочки аппроксимируем полиномом порядка $L$:
\[ \begin{equation}
\Theta(t,\boldsymbol r)-\Theta ^{0} =\sum _{l=0}^{L}\Theta _{l}(t, \boldsymbol x)x_3 ^{l},\quad \boldsymbol x\in \Omega ,\; |x_3|\leqslant h,\; t\geqslant t_{0} ,
\end{equation} \tag{18} \]
где $\Theta _{l} $ ($0\leqslant l\leqslant L$) — неизвестные двумерные нестационарные функции; $\Theta ^{0} ={\rm const}$ — температура естественного состояния КМ-конструкции.

К указанным выше соотношениям необходимо присоединить приведенные уравнения динамического равновесия пологой оболочки, соответствующие уточненной теории изгиба, и уравнения теплового баланса, соответствующие разложению (18). (Их конечно-разностные аналоги по времени выписаны ниже.) А также нужно использовать начальные и краевые механические и тепловые граничные условия, соответствующие этим двумерным нестационарным уравнениям [23].

Как уже ранее отмечалось (см. (5)), численное решение сформулированной нелинейной начально-краевой задачи будем строить, используя явную схему по времени. Согласно этому предполагаем, что при $t=t_{m} $ кроме значений функций (5) уже определены или заданы значения следующих функций [22–24]:
\[ \begin{equation}
\begin{array}{l}
\mathop{w}\limits^{m}(\boldsymbol x)\equiv w(t_{m} , \boldsymbol x),
\quad
\mathop{u_{l}^{(p)} }\limits^{m} (\boldsymbol x)\equiv u_{l}^{(p)} (t_{m} ,\boldsymbol x),
\quad
\mathop{\sigma _{ij} }\limits^{m} (\boldsymbol r)\equiv \sigma _{ij} (t_{m} , \boldsymbol r),
\\
\mathop{\sigma _{33}^{(\pm )} }\limits^{m} (\boldsymbol x)\equiv \sigma _{33}^{(\pm )} (t_{m} , \boldsymbol x),
\quad
\mathop{U^{(\boldsymbol r)} }\limits^{n} (x)\equiv U^{(\boldsymbol r)} (t_{n} , \boldsymbol x),
\quad
\mathop{q_{i} }\limits^{n} (\boldsymbol r)\equiv q_{i} (t_{n} , \boldsymbol r),
\\
\mathop{\Theta _{s} }\limits^{m} (\boldsymbol x)\equiv \Theta _{s} (t_{m} ,\boldsymbol x),
\quad
\mathop{\dot{\Theta }_{s} }\limits^{n-1} (\boldsymbol x)\equiv \dot{\Theta }_{s}(t_{n-1} ,\boldsymbol x),
\quad
\mathop{q_{\infty }^{(\pm )} }\limits^{n}(\boldsymbol x)\equiv q_{\infty }^{(\pm )} (t_{n} , \boldsymbol x),
\\
\mathop{\varepsilon _{ij}^{(k)} }\limits^{m} (\boldsymbol r)\equiv \varepsilon _{ij}^{(k)} (t_{m} ,\boldsymbol r),
\quad
\mathop{e_{ij}^{(k)} }\limits^{m} (\boldsymbol r)\equiv e_{ij}^{(k)} (t_{m} , \boldsymbol r),
\quad
\mathop{\chi _k }\limits^{m} (\boldsymbol r)\equiv \chi _k (t_{m} , \boldsymbol r),
\\
l=1, 2,\; i,j=\overline{1,3},\; m=n-1,n,\; 0\leqslant p\leqslant M+1,\; 0\leqslant r\leqslant L-2,
\\
0\leqslant s\leqslant L,\; 0\leqslant k\leqslant N,\; \boldsymbol x\in \Omega ,\; |x_3|\leqslant h,
\end{array}
\end{equation} \tag{19} \]
где
\[ \begin{equation}
\begin{array}{l}
\displaystyle
u_{l}^{(p)} (t,\boldsymbol x)\equiv \int _{-h}^{h}U_{l}(t, \boldsymbol r)x_3 ^{p} dx_3 ,
\quad U^{(r)}(t,\boldsymbol x)\equiv \int _{-h}^{h}U(t,\boldsymbol r) x_3 ^{r} dx_3 ,
\\
\displaystyle
\sigma _{ij}(t,\boldsymbol r)=\sum _{k=0}^{N}\omega _k (\boldsymbol x)\sigma _{ij}^{(k)}(t, \boldsymbol r) ,\quad
\varepsilon _{ij}(t,\boldsymbol r)=\sum _{k=0}^{N}\omega _k (\boldsymbol x)\varepsilon _{ij}^{(k)} (t,\boldsymbol r) ,
\\
i,j=\overline{1,3},\; l=1, 2,\; 0\leqslant p\leqslant M+1,\; 0\leqslant r\leqslant L-2;
\end{array}
\end{equation} \tag{20} \]
$U$ — удельная внутренняя энергия КМ; $e_{ij}^{(k)} $ — упругие деформации $k$-го материала композиции; $q_{i} $ — компоненты усредненного вектора теплового потока в КМ, связанные с $\operatorname{grad}\Theta $ законом Фурье для армированной среды [22, 24]; $q_{\infty }^{(\pm )} $ — известные тепловые потоки через верхнюю ($+$) и нижнюю ($-$) лицевые поверхности конструкции. В выражениях (2), (3) неизвестные функции $u_{i} $, $\bar{\varepsilon }_{i3}^{(m)} $ ($i=1, 2$, $0\leqslant m\leqslant M$) связаны с $\operatorname{grad} w$ и новыми кинематическими переменными $u_{l}^{(p)} $ (см. (20)) известными матричными соотношениями [23]. Из предположений (19) при учете (18) следуют два последних предположения (5).

Для интегрирования двумерных уравнений движения пологой КМ-оболочки используем численную схему типа «крест», тогда их конечно-разностные аналоги в рамках уточненной теории изгиба будут иметь вид [23]:
\[ \begin{equation}
\begin{array}{l}
\displaystyle
\frac{2h\rho }{\Delta ^{2} } \bigl(\mathop{w}\limits^{n+1} -2\mathop{w}\limits^{n} +\mathop{w}\limits^{n-1} \bigr)= {}
\\
\displaystyle
\hspace{1.3cm}
{} =
\sum _{j=1}^{2}\partial _{j} \biggl(\mathop{M_{j3}^{(0)} }\limits^{\!\!\!\!n} +\sum _{i=1}^{2}\mathop{M_{ji}^{(0)} }\limits^{\!\!\!\!n} \partial _{i} \mathop{w}\limits^{n} \biggr)-
\sum _{i=1}^{2}R_{i}^{-1} \mathop{M_{ii}^{(0)} }\limits^{\!\!\!\!\!\!n} +\mathop{\sigma _{33}^{(+)} }\limits^{\!\!\!\!\!\!\!\!\!\!n} -\mathop{\sigma _{33}^{(-)} }\limits^{\!\!\!\!\!\!\!\!\!\!n},
\\
\displaystyle
\dfrac{\rho }{\Delta ^{2} } \bigl(\mathop{u_{i}^{(l)} }\limits^{\!\!\!n+1} -2\mathop{u_{i}^{(l)} }\limits^{\!\!\!\!\!\!n} +\mathop{u_{i}^{(l)} }\limits^{\!\!\!n-1} \bigr)=
\sum _{j=1}^{2}\partial _{j} \bigl(\mathop{M_{ij}^{(l)} }\limits^{\!\!\!n} -\mathop{M_{j3}^{(l)} }\limits^{\!\!\!n} \partial _{i} \mathop{w}\limits^{n} \bigr)- {}
\\
\displaystyle
\hspace{1.3cm}
{}
-h^{l} \bigl[
\mathop{\sigma _{33}^{(+)} }\limits^{\!\!\!\!\!\!n} -(-1)^{l} \mathop{\sigma _{33}^{(-)} }\limits^{\!\!\!\!\!\!n}
\bigr]\partial _{i} \mathop{w}\limits^{n} -
l\mathop{M_{i3}^{(l-1)} }\limits^{\!\!\!\!\!\!\!\!\!\!\!\!n} +l\mathop{M_{33}^{(l-1)} }\limits^{\!\!\!\!\!\!\!\!\!\!\!\!n} \partial _{i} \mathop{w}\limits^{n} +R_{i}^{-1} \mathop{M_{i3}^{(l)} }\limits^{\!\!\!\!\!\!n} ,
\\
\boldsymbol x\in \Omega ,\; i=1,2,\; 0\leqslant l\leqslant M+1,\; n=1, 2, 3,\dots,
\end{array}
\end{equation} \tag{21} \]
где
\[ \begin{equation}
\begin{array}{l}
\displaystyle
\rho(\boldsymbol x)\equiv \sum _{k=0}^{N}\rho _k \omega _k (\boldsymbol x),
\quad
M_{ij}^{(l)}(t,\boldsymbol x)\equiv \int _{-h}^{h}\sigma _{ij}(t, \boldsymbol r)x_3 ^{l} dx_3 , \quad i,j=\overline{1,3},
\\
\displaystyle
lM_{33}^{(l-1)}(t, \boldsymbol x)\equiv l\int _{-h}^{h} \sigma _{33} (t, \boldsymbol r)x_3 ^{l-1} dx_3 = {}
\\
\hspace{1.3cm}
{}
\displaystyle
=\frac{h^{l} }{2}\Bigl[\bigl(\sigma _{33}^{(+)} +\sigma _{33}^{(-)} \bigr)(1-(-1)^l) +
{}
\\
\hspace{2.6cm}
{}
\displaystyle
+\frac{l}{l+1} \bigl(\sigma _{33}^{(+)} -\sigma _{33}^{(-)} \bigr)(1+(-1)^{l} )\Bigr],\;\; 0\leqslant l\leqslant M+1;
\end{array}
\end{equation} \tag{22} \]
$\rho _k $ — объемная плотность материала $k$-го компонента композиции. По формулам (21) при учете выражений (22) и предположений (19) можно определить значения неизвестных функции $\mathop{w}\limits^{n+1} $ и $\mathop{u_{i}^{(l)} }\limits^{\!\!\!\!\!n+1} $ в следующий момент времени $t_{n+1} $.

Нестационарные двумерные уравнения теплопроводности пологой КМ-оболочки также интегрируются по явной схеме с шагом вперед. Дискретные по времени аналоги этих уравнений имеют вид [23, 24]:
\[ \begin{equation}
\begin{array}{l}
\displaystyle
\frac{\rho }{\Delta } \bigl(\mathop{U^{(m)} }\limits^{\!\!\!n+1} -\mathop{U^{(m)} }\limits^{\!\!\!\!\!\!n} \bigr)=
-\partial _1 \mathop{Q_1 ^{(m)} }\limits^{\!\!\!\!\!\!\!\!\!n}
-\partial _2 \mathop{Q_2 ^{(m)} }\limits^{\!\!\!\!\!\!\!\!\!n}
-\mathop{\bar{Q}_3 ^{(m)} }\limits^{\!\!\!\!\!\!\!\!\!n}
+\mathop{W^{(m)} }\limits^{\!\!\!\!\!\!\!\!\!n} ,
\\
\boldsymbol x\in \Omega ,\; 0\leqslant m\leqslant L-2,\; n=0, 1, 2,\dots;
\end{array}
\end{equation} \tag{23} \]
\[ \begin{equation}
\begin{array}{r}
\displaystyle
-\sum _{l=0}^{L}(-1)^{l} h^{l-1}(l\lambda _{33}^{(-)} +h\alpha ^{(-)})\Theta _{l}(t, \boldsymbol x)=
\alpha ^{(-)}(\Theta _{\infty }^{(-)} -\Theta ^{0})+q_{\infty }^{(-)} (t, \boldsymbol x),
\\
\displaystyle
\sum _{l=0}^{L}h^{l-1}(l\lambda _{33}^{(+)} +h\alpha ^{(+)})\Theta _{l}(t, \boldsymbol x)=
\alpha ^{(+)}(\Theta _{\infty }^{(+)} -\Theta ^{0})-q_{\infty }^{(+)} (t, \boldsymbol x),
\\
\boldsymbol x\in \Omega ,\; t\geqslant t_{0} ;
\end{array}
\end{equation} \tag{24} \]
\[ \begin{equation}
\begin{array}{r}
\displaystyle
C_{0} \sum _{i=0}^{L} H(i+m)\Theta _{i} +\frac{C_1 }{2}
\sum _{i=0}^{L} \sum _{j=0}^{L} H(i+j+m) \Theta _{i} \Theta _{j} + {} \hspace{2.7cm} {}
\\
\displaystyle
{} + \frac{C_2 }{3} \sum _{i=0}^{L} \sum _{j=0}^{L} \sum _{l=0}^{L} H(i+j+l+m) \Theta _{i} \Theta _{j} \Theta _{l} =U^{(m)}(t, \boldsymbol x),
\\
\boldsymbol x\in \Omega ,\; t\geqslant t_{0} ,\; 0\leqslant m\leqslant L-2.
\end{array}
\end{equation} \tag{25} \]
Здесь
\[ \begin{equation}
\!\!\!
\begin{array}{l}
\displaystyle
H(s)\equiv \frac{h^{s+1} }{s+1} [1-(-1)^{s+1}],
\quad
Q_{i}^{(m)}(t, \boldsymbol x)\equiv \int _{-h}^{h}q_{i}(t, \boldsymbol r)x_3 ^{m}dx_3 , i=\overline{1, 3},
\\
\displaystyle
\bar{Q}_3 ^{(m)}(t, \boldsymbol x)\equiv \int _{-h}^{h}\partial _3 q_3 (t,\boldsymbol r)x_3 ^{m} dx_3 =
h^{m}[q_3 ^{(+)}-(-1)^{m} q_3 ^{(-)}]-mQ_3 ^{(m-1)}(t, \boldsymbol x),
\\
\displaystyle
W^{(m)}(t, \boldsymbol x)\equiv\int _{-h}^{h}\sigma _{ij} \dot{\varepsilon }_{ij} x_3 ^{m} dx_3 ,
\quad
C_{l}(\boldsymbol x)\equiv \frac{1}{\rho } \sum _{k=0}^{K}c_{l}^{(k)} \rho _k \omega _k (\boldsymbol x),
\quad l=0, 1, 2,
\\
\lambda _{33}^{(\pm )} \equiv\lambda _{33} \big|_{\Theta =\Theta(t, \boldsymbol x,\pm h)} ,
\quad
q_3 ^{(\pm )} \equiv q_3 (t, \boldsymbol x,\pm h)=q_{\infty }^{(\pm )}(t, \boldsymbol x);
\end{array}
\!\!\!\!\!\!
\end{equation} \tag{26} \]
$\Theta _{\infty }^{(\pm )} $ — температура окружающей среды со стороны верхней ($+$) и нижней ($-$) лицевой поверхности конструкции; $\alpha ^{(\pm )} $ — коэффициент теплоотдачи на той же поверхности; $\lambda _{33} $ — коэффициент теплопроводности КМ в направлении $x_3 $, рассчитанный по формулам (51) из [22]. Предполагается, что в случае учета термочувствительности $k$-го компонента композиции его удельную теплоемкость $c_k $ можно аппроксимировать следующей зависимостью от температуры [22–24]:
\[ \begin{equation}
c_k (\Theta -\Theta ^{0} )=c_{0}^{(k)} +c_1 ^{(k)} \cdot (\Theta -\Theta ^{0} )+
c_2 ^{(k)} \cdot (\Theta -\Theta ^{0} )^{2} ,\quad 0\leqslant k\leqslant N,
\end{equation} \tag{27} \]
где $c_{l}^{(k)} $ — известные характеристики материала. Коэффициенты разложения (27) использованы для вычисления величин $C_{l} $ ($l=0,1,2$) в равенстве (25)
(см. выражения (26)).

Напомним, что равенства (23) — двумерные уравнения теплового баланса, соотношения (24) — граничные условия общего вида на лицевых поверхностях КМ-панели, равенства (25) выражают двумерные функции $U^{(m)} $ (см. (20)) через коэффициенты аппроксимации температуры (18) при учете разложения (27).

По формулам (23) при учете (26) и предположений (19) по явной схеме можно рассчитать значения функций $\mathop{U^{(m)} }\limits^{\!\!\!\!\!\!\!n+1} $ в момент времени $t_{n+1} $. После этого при $t=t_{n+1} $ из замкнутой системы нелинейных уравнений (24) и (25) определяются коэффициенты $\mathop{\Theta _{l} }\limits^{n+1} $ ($0\le l\le L$) в полиномиальном представлении температуры (18). На каждом шаге по времени сначала интегрируется теплофизическая составляющая рассматриваемой задачи (см. (23)–(25) при учете (26) и (27)), а затем при уже известной температуре $\mathop{\Theta }\limits^{n+1} $ — механическая составляющая задачи (см. (21) при учете (3), (12)–(17) и (22)). В остальном эта численная схема реализуется так же, как описано в [22–24].

2. Анализ результатов расчетов

Рассмотрим неизотермическое вязкоупругопластическое изгибное динамическое поведение пологой цилиндрической КМ-оболочки толщиной $2h=2$ см, которая в плане имеет прямоугольную удлиненную форму ($\Omega$: $|x_1|\leqslant a$, $|x_2|\leqslant b$, $a=3b$; $1/R_1\equiv 0$, $R_2\equiv R={\rm const}$; $b=50$ см). Стрела подъема срединной поверхности $f$ над продольными кромками ($x_2 =\pm b$) равна 10 см. Радиус кривизны $R$ при этом выражается так [23]: $R=(b^2+f^2)/(2f)$. Конструкция жестко закреплена по всем кромкам: $u_{i}^{(m)} =0$, $w=0$, $\boldsymbol x\in \Gamma $ и $t\geqslant t_{0} $ (см. (20) и (21)), где $\Gamma$ — контур, ограничивающий область $\Omega $. До момента времени $t=t_{0} =0$ искривленная панель покоится ($u_{i}^{(m)} =0$, $w=0$, $\boldsymbol x\in \Omega $, $t\leqslant t_{0} $, $i=1, 2$, $0\leqslant m\leqslant M+1$) либо при температуре естественного состоянии $\Theta =\Theta^0=20\,^{\circ}$C ($\boldsymbol x\in \Omega $, $|x_3 |\leqslant h$ и $t\leqslant t_{0} $), либо предварительно нагрета со стационарным распределением температуры (подробнее см. ниже).

При $t=t_{0} =0$ пологая КМ-оболочка начинает испытывать механическое нагружение со стороны верхней или нижней лицевой поверхности, которое условно соответствует давлению в воздушной взрывной волне [33]:
\[ \begin{equation}
\begin{array}{l}
p(t)=\begin{cases}
p_{\max } t/t_{{\rm max}} , & 0\leqslant t\leqslant t_{\max} ,
\\
p_{\max} \exp[-\beta(t-t_{\max})], & t>t_{\max} ,
\end{cases}
\\
\sigma _{33}^{(-)}(t)=
\begin{cases} -p(t), & p_{\max}>0, \\
0, & p_{\max}<0,
\end{cases}
\quad
\sigma _{33}^{(+)}(t)=
\begin{cases} 0, & p_{\max}>0,
\\
p(t), & p_{\max} <0.
\end{cases}
\end{array}
\end{equation} \tag{28} \]
Здесь
\[ \begin{equation}
\beta =-\ln(0.01)/(t_{\min} -t_{\max})>0,\quad t_{\min}\gg t_{\max};
\end{equation} \tag{29} \]
$t_{\max}$ — время, при котором $|p(t)|$ достигает максимума $|p_{\max}|$; $t_{\min}$ — время, при превышении которого уже можно не учитывать $|p(t)|$ по сравнению с $|p_{\max}|$ (выражение (29) получено при условии $p(t_{\min})=0.01p_{\max}$). Согласно равенствам (28), при $p_{\max}>0$ КМ-конструкция испытывает нагружение со стороны нижней лицевой поверхности, которая вогнута, а при $p_{\max} <0$ — со стороны верхней поверхности. Исходя из экспериментальных данных [33] примем $t_{\max} =0.1$ мс, $t_{\min}=2$ мс, $|p_{\max}|=6$ МПа.

Конструкция изготовлена из эпоксисвязующего [2] и усилена стекловолокнами [3]. Диаграмма мгновенного упругопластического активного деформирования $k$-го компонента композиции при постоянной температуре $\Theta $ аппроксимируется билинейной зависимостью:
\[ \begin{equation*}
\sigma=\begin{cases}
E_k \varepsilon , & |\varepsilon|\le\varepsilon _{{\rm s}}^{(k)} =\sigma _{{\rm s}}^{(k)} /E_k ,
\\
\operatorname{sign}(\varepsilon)\sigma_{{\rm s}}^{(k)} +E_{{\rm s}}^{(k)} (\varepsilon -\operatorname{sign}(\varepsilon)\varepsilon _{{\rm s}}^{(k)}), & |\varepsilon|>\varepsilon _{{\rm s}}^{(k)} , \; 0\leqslant k\leqslant N,
\end{cases}
\end{equation*} \]
где $\sigma $ и $\varepsilon $ — напряжение растяжения–сжатия и соответствующая деформация; $E_{{\rm s}}^{(k)} =E_{{\rm s}}^{(k)}(\Theta)$ и $\sigma _{{\rm s}}^{(k)} =\sigma _{{\rm s}}^{(k)}(\Theta)$ — модуль линейного упрочнения и предел текучести. Физико-механические характеристики материалов композиций представлены в табл. 1, где $\lambda _k $ — коэффициент теплопроводности $k$-го компонента. (Объемная вязкость фаз композиции не учитывается: $\mu _k \to \infty $, $0\leqslant k\leqslant N$; см. (10).) В расчетах использовались линейные зависимости этих характеристик от температуры $\Theta $, аппроксимированные по данным, приведенным в табл. 1.

Таблица 1. Физико-механические характеристики компонентов композиции [2, 3]
[Physical and mechanical characteristics of the composition components [2, 3]]

Material Сharacteristics	Epoxy Binder ($k=0$)		Glass Fibers ($k=1,2$)
	$\Theta =20 \,^{\circ }$C	$\Theta =100\,^{\circ}$C	$\Theta =20\,^{\circ}$C	$\Theta =100\,^{\circ}$C
$\rho _k $, kg/m$^3$	1210.0	1208.0	2520.0	2519.6
$E_k $, GPa	2.8	2.6	86.8	86.5
$\nu _k $	0.330	0.333	0.250	0.254
$\eta _k $, $\rm MPa \cdot s$	340	300	1250	1200
$\sigma_{{\rm s}}^{(k)}$, MPa	20	15	4500	4400
$E_{{\rm s}}^{(k)}$, GPa	1.114	0.763	6.230	6.079
$\lambda _k $, $\rm W / m \cdot K$	0.243	0.236	0.89	0.86
$\alpha _k \cdot10^{6}$, ${\rm K}^{-1} $	68.1	73.2	2.5	2.6
$c_k $, $\rm kJ /kg \cdot K$	1.54	1.71	0.80	0.84

Панель армирована по ортогональным направлениям $x_1$, $x_2$ двумя (${N=2}$) семействами стеклянных волокон (см. рис. 1) с плотностями армирования $\omega_1=0.1$ и $\omega_2=0.3$ соответственно. В этом случае углы армирования (см. рис. 2) имеют следующие значения: $\theta _1 =\theta _2 =\pi /2$, $\varphi _1 =0$, $\varphi _2 =\pi /2$.

На кромках КМ-панели поддерживается температура $\Theta^0$. Через лицевые поверхности теплообмен с окружающей средой осуществляется в условиях естественной конвекции ($q_{\infty }^{(\pm )} \equiv 0$ и $\alpha ^{(\pm )} =30$ Вт${}/{}$м${}\cdot{}$К [35]) при температуре воздуха $\Theta _{\infty }^{(\pm )} =\Theta ^{0} $ (см. (24)). Кроме того, дополнительно оболочка может быть предварительно нагрета внешним стационарным тепловым потоком через нижнюю ($x_3 =-h$: $q_{\infty }^{(-)} =3$ кВт/м$^2$, $q_{\infty }^{(+)} \equiv 0$) или верхнюю ($x_3 =h$: $q_{\infty }^{(+)} =-3$ кВт/м$^2$, $q_{\infty }^{(-)} \equiv 0$) лицевую поверхность. В этих случаях предполагается, что при $t<t_{0} $ конструкция квазистатически деформируется под действием только теплового нагружения и к моменту времени $t=t_{0} =0$ достигает стационарного состояния. Предварительное термомеханическое состояние при этом можно определить, используя, например, метод установления. Так как рассматриваемая конструкция является относительно тонкой ($h/b=1/50$) и ее искривленность мала, при таком предварительном нагреве к моменту времени $t=t_{0} =0$ за пределами пограничных слоев, которыми здесь пренебрегаем, стационарная температура по толщине панели изменяется практически по линейному закону. При этом наибольшее значение температуры $\Theta _{\max }^{0} \approx 93\,^{\circ }$C достигается на той лицевой поверхности, через которую подается тепловой поток, а наименьшее значение $\Theta _{\min }^{0} \approx 47\, ^{\circ }$C — в точках противоположной лицевой поверхности. Таким образом, при наличии дополнительного теплового потока через лицевые поверхности пологой КМ-оболочки предварительно наведенное температурное поле существенно неоднородно по ее толщине.

Дискретизация задачи по координатам $x_1 $ и $x_2 $ проводилась с использованием равномерной сетки с шагами $\Delta x_1 =\Delta x_2 =b/50=1$ см, шаг же по времени $\Delta =1$ мкс. При этом необходимые условия устойчивости используемой численной схемы (условия Куранта) выполняются с запасом (см. (82) в [22]).

На основании результатов работ [23, 24] порядок полинома в аппроксимации температуры (18) принят равным $L=7$ (при $L\geqslant 8$ линеаризованная система (24) и (25) плохо обусловлена, поэтому решение задачи расходится). В настоящем исследовании изучается влияние применения уточненной теории изгиба (см. (2), (3) и (21) при $M\neq 0$ [23]) на результаты расчета остаточного термомеханического состояния гибкой пологой КМ-оболочки после окончания ее осцилляций, вызванных приложением кратковременной, но интенсивной динамической нагрузки. Поэтому на рис. 3 представлены зависимости от времени наибольших значений температуры $\Theta _{{\rm m}}(t;M)=\max\limits_{\boldsymbol r}\Theta(t, \boldsymbol r;M)$ ($|x_1|\leqslant a$, $|x_2|\leqslant b$ и $|x_3|\leqslant h$) в рассматриваемой конструкции при отсутствии (рис. 3, a и 3, b: $q_{\infty }^{(\pm )}\equiv 0$) и наличии дополнительного теплового потока через верхнюю лицевую поверхность (рис. 3, c: $q_{\infty }^{(+)} =-3$ кВт/м$^2$, $q_{\infty }^{(-)} \equiv 0$). Кривые 1 и 2 получены соответственно при $M=7$ (уточненная теория изгиба) и при $M=0$ (теория Амбарцумяна). Выбор значения $M=7$ в соотношениях (2), (3) и (21) в рамках использования уточненной теории изгиба обоснован в [23]. Зависимости $\Theta _{{\rm m}}(t; M)$ на рис. 3, a соответствуют случаю, когда искривленная панель нагружается снизу ($p_{\max }=6$ МПа; см. (28), (29)) — со стороны вогнутой лицевой поверхности, а на рис. 3, b и 3, c — случаям нагружения сверху ($p_{\max }=-6$ МПа). Из поведения кривых на рис. 3 видно, что осцилляции максимальных значений температуры $\Theta _{{\rm m}} $ к моменту времени $t=100$ мс существенно уменьшаются, а при $t\approx 200$ мс их можно считать полностью затухшими (см. рис. 3 в [22]). Кривые 1 и 2 на рис. 3 заметно различаются, особенно при нагружении КМ-панели давлением (28) сверху (см. рис. 3, b и 3, c), т. е. использование уточненной теории изгиба (кривые 1) действительно приводит к уточнению температурного отклика в пологой армированной оболочке. Кривые на рис. 3, a в основном расположены ниже аналогичных кривых на рис. 3, b. Следовательно, температурный отклик в рассматриваемой конструкции является более интенсивным при ее механическом нагружении со стороны выпуклой лицевой поверхности. В частности, при $p_{\max }=-6$ МПа дополнительный нагрев искривленной панели может превышать $12\,^{\circ }$C (см. кривые на рис. 3, b и 3, c при $t\approx 9.8$ мс), а при $p_{\max } =6$ МПа не превышает $7\,^{\circ }$C (см. кривые на рис. 3, a при $t\approx 12.4$ мс).

Рис. 3. Зависимость от времени наибольших значений температуры в пологой оболочке из стеклопластика: a — случай $q_{\infty }^{(\pm )} \equiv 0$ и $p_{\max } >0$; b — случай $q_{\infty }^{(\pm )} \equiv 0$ и $p_{\max } <0$; c — случай $q_{\infty }^{(+)} =-3$ кВт/м$^2$, $q_{\infty }^{(-)} \equiv 0$ и $p_{\max } <0$
[Figure 3. Time dependence of the highest temperature values in a shallow fiberglass shell: a — case $q_{\infty }^{(\pm )} \equiv 0$ and $p_{\max } >0$; b — case $q_{\infty }^{(\pm )} \equiv 0$ and $p_{\max } <0$; c — case $q_{\infty }^{(+)} =-3$ kW/m$^2$, $q_{\infty }^{(-)} \equiv 0$ and $p_{\max } <0$]

Несмотря на то, что температурный отклик в исследуемой стеклопластиковой оболочке при рассматриваемом ее динамическом нагружении невелик (дополнительный нагрев в отдельных точках не превосходит $7 \div 12\,^{\circ}$C), он все же может оказывать заметное влияние на остаточный прогиб такой КМ-панели. Так, на рис. 4 представлены зависимости от переменной $x_2$ прогиба $w$ в центральном сечении удлиненной конструкции ($x_1 =0$), рассчитанные при отсутствии (рис. 4, a) и наличии дополнительного теплового потока через одну из лицевых поверхностей: верхнюю (рис. 4, b) или нижнюю (рис. 4, c). Кривые на рис. 4 получены при $t=500$ мс, когда осцилляции пологой КМ-оболочки практически прекратились.

Рис. 4. Зависимости остаточных прогибов от координаты $x_2$ в центральном сечении $x_1 =0$ стеклопластиковой панели: a — случай $q_{\infty }^{(\pm )} \equiv 0$; b — случай $q_{\infty }^{(+)} =-3$ кВт/м$^2$ и $q_{\infty }^{(-)} \equiv 0$; c — случай $q_{\infty }^{(-)} =3$ кВт/м$^2$, $q_{\infty }^{(+)} \equiv 0$
[Figure 4. Dependences of residual deflections on coordinate $x_2 $ in the central section $x_1 =0$ of a fiberglass panel: a — case $q_{\infty }^{(\pm )} \equiv 0$; b — $q_{\infty }^{(+)} =-3$ kW/m$^2$ and $q_{\infty }^{(-)} \equiv 0$; c — case $q_{\infty }^{(-)} =3$ kW/m$^2$ and $q_{\infty }^{(+)} \equiv 0$]

Кривые с номерами 1 (со штрихами и без них) на рис. 4 соответствуют случаям механического нагружения панели (28) со стороны вогнутой лицевой поверхности, а кривые с номерами 2 — со стороны выпуклой поверхности. Кривые с номерами без штрихов рассчитаны по уточненной теории изгиба, а кривые, номера которых помечены одним штрихом, — по теории Амбарцумяна. Кривые 1$''$ и 1$''$ на рис. 4, a рассчитаны при тех же условиях, что и кривые 1 и 2, но без учета температурного отклика в армированной панели (вязкоупругопластическое деформирование).

Кривые 1 и 2 на рис. 4 существенно отличаются от кривых 1$'$ и 2$'$ соответственно. Следовательно, использование простейшей традиционной неклассической теории изгиба Амбарцумяна (кривые 1$'$ и 2$'$) приводит к существенному искажению расчетного остаточного прогиба рассматриваемой КМ-оболочки по сравнению с расчетом по уточненной теории изгиба (кривые 1 и 2). Сравнение кривых 2 и 2$'$ на рис. 4 свидетельствует о том, что это искажение носит не только количественный, но и качественный характер. Этот вывод относится и к кривым 1, 1$'$ на рис. 4, a. Кривые на рис. 4, b и 4, c существенно отличаются от кривых с такими же номерами на рис. 4, a. А значит, дополнительное тепловое воздействие немеханического происхождения значительно влияет на величину остаточного прогиба рассматриваемой искривленной панели. Кроме того, сравнение кривых 2 на рис. 4, b и 4, c показывает, что при механическом нагружении пологой оболочки сверху на ее остаточный прогиб заметное влияние оказывает и выбор лицевой поверхности, через которую к КМ-конструкции подводится дополнительный тепловой поток. Так, максимальные по модулю ординаты точек на кривой 2 рис. 4, b на 5.9 % больше аналогичных величин на кривой 2 рис. 4, c. Сравнение же максимальных по модулю ординат точек на кривых 1 на рис. 4, b и 4, c показывает, что они различаются всего на 0.1 %. Следовательно, при нагружении искривленной панели избыточным давлением снизу выбор лицевой поверхности, через которую подводится дополнительный тепловой поток, в рассматриваемой задаче практически не оказывает влияния на форму и величину ее остаточного прогиба.

Сопоставление кривых 1 и 1$''$ на рис. 4, a демонстрирует, что неучет теплового отклика в КМ-конструкции (кривая 1$''$) может привести к качественно неверному представлению о форме остаточного прогиба пологой стеклопластиковой оболочки. Так, кривая 1 на рис. 4, a в центральной точке ($x_2 =0$) выпукла вверх, а кривая 1$''$ — вниз, хотя максимальные по модулю прогибы при этом различаются незначительно (всего на 2.7 %). Сравнение же кривых 2 и 2$''$ на рис. 4, a показывает, что при механическом нагружении искривленной панели сверху неучет температурного отклика в ней (кривая 2$''$) приводит к качественно верному представлению о форме ее остаточного прогиба, однако дает значительную ошибку при оценке максимального прогиба. Так, наибольшие по модулю значения ординат точек на кривых 2 и 2$''$ рис. 4, a различаются уже на 7.2 %.

Использование теории изгиба Амбарцумяна может приводить не только к существенному искажению расчетного остаточного прогиба пологой КМ-оболочки, но и к значительному искажению представления о деформированном состоянии компонентов ее композиции. С целью демонстрации этого факта были рассмотрены зависимости от времени максимальных значений интенсивности деформаций $\varepsilon _{*}^{(k)} $ в $k$-х фазах композиции стеклопластиковой панели ($\varepsilon _{{\rm m}}^{(k)}(t;M)=\max\limits_{\boldsymbol r} \varepsilon _{*}^{(k)}(t, \boldsymbol r; M)$, $|x_1|\leqslant a$, $|x_2|\leqslant b$, $|x_3|\leqslant h$, $k=0, 1, 2$). В табл. 2 приведены некоторые характерные значения величин $\varepsilon _{{\rm m}}^{(k)} $ для эпоксисвязующего ($k=0$) и волокон второго семейства (${k=2}$), которые в процессе осцилляций КМ-конструкции испытывают наиболее интенсивное деформирование. В табл. 2 приняты следующие обозначения: $\varepsilon _{\max}^{(k)}(M) =\max\limits_{t\geqslant 0} \varepsilon _{{\rm m}}^{(k)}(t;M)$ и $\varepsilon _{0.5}^{(k)}(M)=\varepsilon _{{\rm m}}^{(k)}(0.5; M)$, т. е. величины $\varepsilon _{0.5}^{(k)} $ можно трактовать как максимальные значения интенсивности остаточных деформаций, так как они соответствуют значениям $\varepsilon _{{\rm m}}^{(k)} $ ($k=0, 2$), полученным в момент времени $t=500$ мс, когда осцилляции искривленной панели практически полностью прекратились. Данные, представленные в табл. 2 при $L=0$, соответствуют случаям, когда температурный отклик в конструкции не учитывается (вязкоупругопластический расчет).

Таблица 2. Характерные наибольшие значения интенсивности деформаций некоторых компонентов композиции
[Characteristic maximum values of the intensity of deformations of some components of the composition]

nos.	$q_{\infty }^{(-)}$, kW/m$^2$	$q_{\infty }^{(+)}$, kW/m$^2$	$p_{\max}$, MPa	$M$	$L$	$\varepsilon _{\max }^{(0)}$, %	$\varepsilon _{0.5}^{(0)}$, %	$\varepsilon _{\max }^{(2)}$, %	$\varepsilon _{0.5}^{(2)} $, %
1	0	0	6	0	7	5.16	1.69	2.55	0.493
2	0	0	6	7	7	5.16	1.79	2.83	0.731
3	0	0	6	7	0	5.16	1.83	2.83	0.736
4	0	0	$-6$	0	7	6.48	2.18	3.53	0.644
5	0	0	$-6$	7	7	6.19	2.25	4.00	0.673
6	0	0	$-6$	7	0	6.17	2.34	3.96	0.654
7	3	—	6	0	7	5.31	1.87	2.58	0.596
8	3	—	6	7	7	5.08	1.77	2.84	0.731
9	3	—	$-6$	0	7	6.54	2.23	3.39	0.779
10	3	—	$-6$	7	7	6.24	2.27	3.89	0.711
11	—	$-3$	6	0	7	5.32	1.87	2.58	0.549
12	—	$-3$	6	7	7	5.22	1.86	2.80	0.699
13	—	$-3$	$-6$	0	7	6.90	2.14	3.43	0.698
14	—	$-3$	$-6$	7	7	6.38	2.33	3.84	0.674

Сравнение значений $\varepsilon _{\max }^{(2)} $ в расчетах с номерами 4 и 5 (см. табл. 2) показывает, что использование теории Амбарцумяна в этом случае приводит к занижению максимальной величины интенсивности деформаций второго семейства волокон на 11.3 %, в расчетах 13 и 14 — на 12 %, а в расчетах 9 и 10 — на 14.7 %. В еще большей степени теория Амбарцумяна может приводить к искажению расчетного остаточного деформированного состояния фаз композиции. Так, значение $\varepsilon _{0.5}^{(2)}$ в расчете 1 на 48.3 % меньше, чем в расчете 2.

Сравнение величин $\varepsilon _{\max }^{(k)} $ и $\varepsilon _{0.5}^{(k)} $ ($k=0, 2$) в расчетах 2, 3 и 5, 6 демонстрирует, что неучет теплового отклика в пологой КМ-оболочке (см. расчеты 3 и 6) не оказывает существенного влияния на изменение этих величин, в отличие от его влияния на величину и форму остаточного прогиба. Так, значения $\varepsilon _{0.5}^{(0)} $ в расчетах 5 и 6 различаются всего на 4 %, а остальные величины $\varepsilon _{\max }^{(k)} $ и $\varepsilon _{0.5}^{(k)} $ в указанных расчетах различаются менее чем на 3 %.

Заключение

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

Полученные результаты показали, что, как и в [23], для расчета неизотермического вязкоупругопластического динамического поведения искривленных стеклопластиковых панелей необходимо использовать уточненную теорию их изгиба. Использование же простейшего варианта этой теории — традиционной неклассической теории Амбарцумяна [22, 24, 28] — может привести к существенному как количественному, так и качественному искажению расчетного остаточного прогиба и, в еще большей степени, расчетного остаточного деформированного состояния компонентов композиции по сравнению с расчетами, выполненными по уточненной теории изгиба [23].

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

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

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

Конкурирующие интересы. У меня нет конфликта интересов в авторстве и публикации этой статьи.
Авторская ответственность. Я несу полную ответственность за предоставление окончательной версии рукописи в печать. Окончательная версия рукописи мною одобрена.
Финансирование. Работа выполнена в рамках государственного задания (№ госрегистрации 121030900260-6).

About the authors

Andrei P. Yankovskii

Khristianovich Institute of Theoretical and Applied Mechanics Siberian Branch of the Russian Academy of Sciences

Author for correspondence.
Email: yankovsky_ap@itam.nsc.ru
ORCID iD: 0000-0002-2602-8357
SPIN-code: 9972-3050
Scopus Author ID: 7003288442
ResearcherId: J-9106-2018
https://www.mathnet.ru/person28373

Dr. Phys. & Math. Sci.; Leading Research Scientist; Lab. of Fast Processes Physics

Russian Federation, 630090, Novosibirsk, Institutskaya st., 4/1

References

Supplementary files