Modeling of non-isothermal elastic-plastic behavior of reinforced shallow shells in the framework of a refined bending theory

Full Text

Введение

Изделия из композитных материалов (КМ) находят широкое применение в современной инженерной практике [1–15] и часто испытывают интенсивное термосиловое нагружение [2, 3, 5, 8–10, 12, 15], при котором компоненты композиции деформируются неупруго [5, 8, 12, 16–18]. Следовательно, актуальной является проблема математического моделирования термоупру-гопластического деформирования КМ-пластин и КМ-оболочек, находящаяся на сегодняшний день в стадии становления [5, 8, 12, 19–23]. Так, в работе [23] разработана модель неизотермического упругопластического деформирования материала, многонаправленно армированного непрерывными волокнами, и выполнены расчеты динамического изгиба КМ-пластин из таких материалов.

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

Для интегрирования нелинейных задач динамики тонкостенных элементов конструкций используют как явные [5, 22, 23], так и неявные [8, 26] численные схемы.

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

1. Постановка задачи и метод расчета

Рассмотрим деформирование гибкой пологой КМ-оболочки толщиной $2h$, с которой свяжем ортогональную систему координат $x_i$ так, что ее срединная поверхность является отсчетной ($\big|x_3 \big|\leqslant h$), а координатные линии $x_1$ и $x_2$ совпадают с линиями главной кривизны этой поверхности $x_3=0$ (рис. 1, где искривленность панели не изображена в силу ее малости). Оболочка многонаправленно армирована (возможно, и пространственно) $K$ семействами волокон с плотностями армирования $\omega _k$ ($1\leqslant k\leqslant K$). В поперечном направлении конструкции $Ox_3$ структура армирования однородна.

Рис. 1. Элемент волокнистой пологой оболочки
[Figure 1. Element of fibrous shallow shell]

Рис. 2. Локальная прямоугольная система координат, связанная с армирующим волокном $k$-го семейства
[Figure 2. Local cartesian coordinate system associated with $k$-th family reinforcing fiber]

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

Предполагается, что внешними касательными силами на лицевых поверхностях оболочки можно пренебречь, и в случае пространственной структуры армирования выполняются требования, изложенные в замечании в [23]. (В случаях плоско-перекрестного армирования, как на рис. 1, указанные требования заведомо выполняются.) При этом, согласно [22], усредненные деформации композиции $\varepsilon_{ij}$ и перемещения точек гибкой пологой оболочки $U_i$ в рамках уточненной теории изгиба можно аппроксимировать так (геометрическая нелинейность задачи учитывается в приближении Кармана):
\[ \begin{equation}
\varepsilon_{ij} \big( t,{\bf r} \big)=\frac{1}{2}\big(\partial_i u_j +\partial _j u_i \big) - x_3 \partial _i \partial _j w
+ \sum\limits_{m=0}^M \frac{x_3^{m+1}}{h^2} \Bigl(\frac{h^2}{m+1}-\frac{x_3^2 }{m+3}\Bigr) 
\bigl(\partial_i \varepsilon_{j3}^{(m)} + \partial_j \varepsilon _{i3}^{(m)}\bigr)+\frac{\delta _{ij} w}{R_i }+\frac{1}{2}\partial _i w\partial _j w,
\label{yank:eq1}
\end{equation} \tag{2} \]
\[ \begin{equation*}
\varepsilon_{i3} \big(t,{\bf r} \big)=\frac{h^2-x_3^2 }{h^2}\sum\limits_{m=0}^M x_3^m \varepsilon 
_{i3}^{(m)} \big(t,{\bf x} \big),\quad i,j=1,2; 
\end{equation*} \]
\[ \begin{equation}
\begin{aligned}
& U_i \big(t,{\bf r} \big)=u_i \big(t,{\bf x} \big) - x_3 \partial_i w + 2 \sum\limits_{m=0}^M \frac{x_3^{m+1}}{h^2} 
\Bigl( \frac{h^2}{m+1}-\frac{x_3^2}{m+3}\Bigr)\varepsilon _{i3}^{(m)} \big(t,{\bf x} \big), \\ 
& U_3 \big( t,{\bf r} \big) = w\big(t,{\bf x}\big),\quad i=1,2,\quad {\bf x} \in G,\\
&\big| x_3 \big|\leqslant h, \quad t \geqslant t_0,\quad {\bf x}=\big( x_1, x_2 \big),\quad {\bf r}=\big(x_1, x_2, x_3 \big), 
\end{aligned}
\end{equation} \tag{3} \]
где $w$ — прогиб; $u_i$ — перемещения точек отсчетной поверхности ($x_3 =0$) в тангенциальных направлениях $x_i$; $R_i$ — главные радиусы кривизны срединной поверхности; $t_0$ — начальный момент времени $t$; $\partial _i$ — оператор частного дифференцирования по переменной $x_i$ ($i=1, 2$); $\delta _{ij}$ — символ Кронекера; $M$ — целое число, которое задает количество слагаемых, удерживаемых в частичных суммах по степеням $x_3$; $G$ — область, занимаемая искривленной панелью в плане. При $M=0$ из соотношений (2) и (3) вытекают кинематические гипотезы теорий Амбарцумяна [1] и Редди [6, 9]. В равенствах (2) и (3) подлежат определению двумерные функции $w$, $u_i$, $\varepsilon_{i3}^{(m)}$ ($i=1, 2$; $0\leqslant m\leqslant M$).

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

К соотношениям (2)–(4) следует присоединить двумерные уравнения движения гибкой пологой оболочки (см. (10)–(12) в [22]) и определяющие уравнения, связывающие усредненные скорости деформаций $\dot {\varepsilon }_{ij}$, напряжений $\dot {\sigma }_{ij}$ и температуры $\dot \Theta$ в композиции конструкции. Последние уравнения целесообразно представить в матричном виде [23]:
\[ \begin{equation}
\dot {\boldsymbol \sigma }={\rm {\bf B}}\dot {\boldsymbol \varepsilon }+{\bf p},
\end{equation} \tag{5} \]
где
\[ \begin{equation}
\begin{aligned}
& \dot {\boldsymbol \sigma }= ( {\dot {\sigma }_1 , \dot {\sigma }_2 , \dot 
{\sigma }_3 , \dot {\sigma }_4 , \dot {\sigma }_5 , \dot {\sigma 
}_6 } )^{\top}\equiv ( {\dot {\sigma }_{11} , \dot {\sigma 
}_{22} , \dot {\sigma }_{33} , \dot {\sigma }_{23} , \dot {\sigma 
}_{31} , \dot {\sigma }_{12} } )^{\top}, \\ 
& \dot {\boldsymbol \varepsilon }= ( {\dot {\varepsilon }_1 , \dot {\varepsilon 
}_2 , \dot {\varepsilon }_3 , \dot {\varepsilon }_4 , \dot 
{\varepsilon }_5 , \dot {\varepsilon }_6 } )^{\top}\equiv 
( {\dot {\varepsilon }_{11} , \dot {\varepsilon }_{22} , \dot 
{\varepsilon }_{33} , 2\dot {\varepsilon }_{23} , 2\dot {\varepsilon 
}_{31} , 2\dot {\varepsilon }_{12} } )^{\top};
 \end{aligned}
\end{equation} \tag{6} \]
\[ \begin{equation}
\begin{aligned}
& {\bf B} \equiv \biggl( \omega _0 {\bf Z}_0 +\sum\limits_{k=1}^N {\omega _k } \rm {\bf Z}_k {\bf E}_k \biggr) {\bf H}^{-1},\quad 
{\bf p}\equiv {\bf f}-{\bf Bg}, \\
& {\bf f}\equiv \omega _0 {\bf p}_0 +\sum\limits_{k=1}^N {\omega _k } ( {\bf p}_k +{\bf Z}_k {\bf r}_k ), \quad
{\bf H}\equiv \omega _0 {\bf I}+\sum\limits_{k=1}^N {\omega _k } {\bf E}_k , \\
& {\bf g}\equiv \sum\limits_{k=1}^N {\omega _k } {\bf r}_k ,\quad 
{\bf r}_k \equiv {\bf D}_k^{-1} \varsigma _k , \quad {\bf E}_k \equiv {\bf D}_k^{-1} {\bf C}_k\quad (1\leqslant k\leqslant N), \\ 
& {\bf Z}_k = {\bf \bar {Z}}_k -G_k {\bf \bar {\bar {Z}}}_k , \quad 
{\bf p}_k \equiv \beta _k \dot {\Theta }\quad (0\leqslant k\leqslant N),\quad \omega _0 \equiv 
1-\sum\limits_{k=1}^N {\omega _k } ; 
\end{aligned}
\end{equation} \tag{7} \]
$\textbf{I}$ — единичная $6{\times} 6$-матрица; $\textbf{B}$, ${\bf E}_k$, ${\bf C}_k$, ${\bf Z}_k$, ${\bf \bar {Z}}_k$, ${\bf \bar {\bar {Z}}}_k$ — $6{\times} 6$-матрицы; ${\bf D}_k^{-1} $, ${\bf H}^{-1}$ — матрицы, обратные $6{\times} 6$-матрицам ${\bf D}_k $, $\textbf{H}$; $\textbf{p}$, $\textbf{f}$, $\textbf{g}$, ${\bf r}_k$, $\varsigma_k$, ${\bf p}_k$, $\beta _k$ — шестикомпонентные векторы-столбцы. Элементы матриц ${\bf C}_k = ( c_{ij}^{(k)} )$, ${\bf D}_k = ( d_{ij}^{(k)} )$, ${\bf \bar {Z}}_k = ( {\bar {z}_{ij}^{(k)} } )$, ${\bf \bar {\bar {Z}}}_k = ( \bar {\bar {z}}_{ij}^{(k)} )$ и вектор-столбцов $\varsigma _k = ( \varsigma _i^{(k)} )$, $\beta _k =( \beta _i^{(k)} )$ вычисляются по формулам [23]:
\[ \begin{equation}
\begin{aligned}
& \bar {z}_{ij}^{(k)} =2\delta _{ij} G_k +\lambda _k , \quad \bar {z}_{ll}^{(k)} =G_k , 
\\
&
\beta _i^{(k)} =\frac{K_\Theta ^{(k)} }{3K_k }\sum\limits_{m=1}^3 {\sigma _{mm}^{(k)} } +3K_k \alpha _k 
+\frac{s_i^{(k)} }{G_k }\bigl[ {G_\Theta ^{(k)} } -\tau _s^{(k)} ( \tau _s^{(k)} G_\Theta ^{(k)} -\tau _\Theta ^{(k)} G_k )A_k \bigr],
\\
&
 \beta _l^{(k)} =\frac{s_l^{(k)} }{G_k }\bigl[ G_\Theta ^{(k)} -\tau _s^{(k)} ( \tau 
_s^{(k)} G_\Theta ^{(k)} -\tau _\Theta ^{(k)} G_k )A_k \bigr], \quad (i,j=\overline {1,3}, \; l=\overline {4,6} ),
\\
& 
\bar {\bar {z}}_{ij}^{(k)} =A_k s_i^{(k)} s_j^{(k)}\quad (i,j=\overline {1,6} ), \\
& 
\lambda _k =\frac{\nu _k E_k }{ ( 1+\nu _k ) ( 1-2\nu _k )}, \quad
2G_k =\frac{E_k }{1+\nu _k }, \quad 
3K_k =\frac{E_k }{1-2\nu _k }, \\
&
A_k =\frac{\gamma _k G_k }{ ( G_k +\bar {G}_k )\tau _s^{(k)\, 2} }, \quad 
K_\Theta ^{(k)} =\frac{dK_k }{d\Theta }, \quad G_\Theta ^{(k)} =\frac{dG_k }{d\Theta },
\\
& \tau _\Theta ^{(k)} =\frac{\partial \tau _s^{(k)} }{\partial \Theta }, \quad 
\bar {G}_k =\frac{\partial \tau _s^{(k)} }{\partial \chi _k }, \quad 
\chi _k = \int_{t_0 }^t \sqrt {2\dot {p}_{ij}^{(k)} \dot {p}_{ij}^{(k)} } dt, \\ 
&
 \gamma _k =
 \begin{cases}
 0, & \text{если $T_k <\tau _s^{(k)}$ или $T_k =\tau _s^{(k)}$, $W_k \leqslant 0$}, 
\\ 
 1, & \text{если $T_k =\tau _s^{(k)}$, $W_k >0$}, 
\\ 
 \end{cases} \\ 
&
W_k =G_k {\bf s}_k^{\top} \dot {\boldsymbol \varepsilon }_k +\tau _s^{(k)} G_k^{-1} 
\bigl( \tau _s^{(k)} G_\Theta ^{(k)} -\tau _\Theta ^{(k)} G_k \bigr)\dot {\Theta },
\\
& T_k^2 =\frac{1}{2}\sum\limits_{i=1}^3 s_i^{(k)\, 2} +\sum\limits_{i=4}^6 s_i^{(k)\, 2} , \quad 
 0\leqslant k\leqslant N; 
 \end{aligned}
\end{equation} \tag{8} \]
\[ \begin{equation}
\begin{aligned}
& c_{1j}^{(k)} =d_{1j}^{(k)} =q_{1j}^{(k)} , \quad 
 c_{ij}^{(k)} =\sum\limits_{l=1}^6 g_{il}^{(k)} b_{lj}^{(0)} , \quad 
 d_{ij}^{(k)} =\sum\limits_{l=1}^6 g_{il}^{(k)} b_{lj}^{(k)} , 
\\
& 
\varsigma _1^{(k)} =0, \quad 
\varsigma _i^{(k)} =\sum\limits_{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{aligned}
\end{equation} \tag{9} \]
\[ \begin{equation}
\begin{aligned}
& 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{aligned}
\end{equation} \tag{10} \]
$\alpha _k$ — коэффициент линейного температурного расширения $k$-го компонента композиции ($k=0$ — связующая матрица, $k\geqslant 1$ — волокна $k$-го семейства); $E_k$, $\nu_k$ — модуль Юнга и коэффициент Пуассона того же материала; $\tau _s^{(k)} =\tau _s^{(k)} ( \chi _k ,\Theta )$ — предел текучести того же материала при чистом сдвиге, зависящий от параметра Одквиста $\chi_k$ и температуры $\Theta$ и при активном пластическом деформировании равный интенсивности касательных напряжений $T_k$; $\bar {G}_k$ — касательный модуль сдвига при пластическом деформировании $k$-й фазы композиции; $\sigma _{ij}^{(k)}$ —  напряжения в том же материале; ${\bf s}_k = ( s_l^{(k)} )$, $\dot {\boldsymbol \varepsilon }_k = ( \dot {\varepsilon }_l^{(k)} )$ — шестикомпонентные векторы-столбцы с компонентами девиатора напряжений $s_{ij}^{(k)}$ и тензора скоростей деформаций $\dot {\varepsilon }_{ij}^{(k)}$ в $k$-м материале композиции, по структуре аналогичные векторам $\dot {\sigma }$ и $\dot {\varepsilon }$ соответственно (см. (6)); $\dot {p}_{ij}^{(k)}$ — скорости пластических деформаций того же материала; $\gamma_k$ — параметр переключения, который при $\gamma_k =0$ определяет термоупругое деформирование, разгрузку или нейтральное нагружение, а при $\gamma_k =1$ — активное нагружение при пластическом деформировании $k$-го компонента композиции; $p_l^{(k)}$ — элементы вектора-столбца ${\bf p}_k$ ($l=\overline {1,6} $, $0\leqslant k\leqslant N$); точка — производная по времени $t$; индекс $^\top$ — операция транспонирования. По повторяющимся индексам $k$ и $l$ в равенствах (8) суммирования нет. Не выписанные в (10) элементы $6{\times} 6$-матриц ${\bf G}_k = ( g_{ij}^{(k)} )$ и ${\bf Q}_k = ( q_{ij}^{(k)} )$, зависящие от направляющих косинусов (1), приведены в табл. 21.40 и 21.44 в [27]. Матрицы ${\bf G}_k$ и ${\bf Q}_k$ определяют преобразования шестикомпонентных векторов-столбцов скоростей напряжений $\dot {\boldsymbol \sigma }_k = ( \dot {\sigma }_l^{(k)} )$ и деформаций $\dot {\boldsymbol \varepsilon }_k = ( \dot {\varepsilon }_l^{(k)} )$ ($l=\overline {1,6}$) в $k$-м компоненте композиции при переходе от глобальной ортогональной системы координат $x_j$ к локальной системе $x_i^{(k)}$ (см. рис. 2). Векторы $\dot {\boldsymbol \sigma }_k$ и $\dot {\boldsymbol \varepsilon }_k$, имеющие соответственно структуру, аналогичную (6), связаны определяющими соотношениями [23]:
\[ \begin{equation}
\dot {\boldsymbol \sigma }_k ={\bf Z}_k \dot {\boldsymbol \varepsilon }_k +{\bf p}_k 
,\quad 0\leqslant k\leqslant N,
\end{equation} \tag{11} \]
где матрица ${\bf Z}_k$ и вектор-столбец ${\bf p}_k$ вычисляются по формулам (7). В соотношениях (8) учтено, что упругие характеристики $k$-го материала композиции могут зависеть от температуры $\Theta$ (термочувствительность).

При выводе равенств (5) и (7) попутно получаются следующие структурные матричные соотношения:
\[ \begin{equation}
\dot {\boldsymbol \varepsilon }_0 ={\bf H}^{-1}\dot {\boldsymbol \varepsilon }- {\bf H}^{-1}{\bf g},
\quad 
\dot {\boldsymbol \varepsilon }_k ={\bf E}_k \dot {\boldsymbol \varepsilon }_0 +{\bf r}_k ,\quad 1\leqslant k\leqslant N.
\end{equation} \tag{12} \]
Первое выражение (12) определяет скорости деформаций связующей матрицы $\dot {\boldsymbol \varepsilon }_0$ через скорости осредненных деформаций композиции $\dot {\boldsymbol \varepsilon }$, а второе соотношение — скорости деформаций волокон $k$-го семейства $\dot {\boldsymbol \varepsilon }_k$ через скорости деформаций связующего материала $\dot {\boldsymbol \varepsilon }_0$. По формулам же (11) при учете (12) можно выразить и скорости напряжений $\dot {\sigma }_{ij}^{(k)}$ в $k$-м компоненте композиции через осредненные скорости деформаций $\dot {\varepsilon }_{ij}$ ($i,j=\overline {1,3}$, $0\leqslant k\leqslant N$).

В силу структуры векторов-столбцов $\dot {\boldsymbol \sigma }$ и $\dot {\boldsymbol \varepsilon }$ (см. (6)) из третьего равенства системы (5) определим скорость линейной деформации в направлении $x_3$:
\[ \begin{equation}
\dot {\varepsilon}_{33} =b_{33}^{-1} \big(\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} \big),
\end{equation} \tag{13} \]
где $p_3$, $b_{3i}$ ($i=\overline {1,6}$) — элементы вектора-столбца $\textbf{p}$ и матрицы $\textbf{B}$ в соотношении (5); скорость $\dot {\sigma }_{33}$ выражается из (4) путем дифференцирования по времени. Скорости деформаций $\dot {\varepsilon }_{ij}$ в правой части (13) вычисляются за счет дифференцирования по $t$ соотношений (2), т. е. определяются через двумерные функции $w$, $\dot {w}$, $\dot {u}_i$, $\dot {\varepsilon }_{i3}^{(m)}$ ($i=1,2$, $0\leqslant m\leqslant M$).

Используя результаты работы [23], температуру конструкции $\Theta$ аппроксимируем по ее толщине полиномом порядка $L$:
\[ \begin{equation}
\Theta (t, {\bf r} )-\Theta^0 = \sum\limits_{l=0}^L 
\Theta _l (t, {\bf x} ) x_3^l,\quad {\bf x}\in G,\quad |x_3 | \leqslant 
h,\quad t\geqslant t_0,
\end{equation} \tag{14} \]
где $\Theta^0=\rm{const}$ — температура естественного состояния искривленной панели; $\Theta_l$ — подлежащие определению двумерные функции ($0\leqslant l\leqslant L$).

Для завершения формулировки связанной задачи неупругого термомеханического деформирования пологой КМ-оболочки к соотношениям (2)–(14) нужно присоединить уравнения теплового баланса. Так как метрика в искривленной панели с приемлемой для приложений точностью может быть отождествлена с метрикой в декартовой прямоугольной системе координат [1], в качестве двумерных уравнений, описывающих теплопроводность в рассматриваемой пологой КМ-оболочке, можно применять формулы (4.11), (4.13), (4.14) и (4.18) из работы [23]. Кроме того, следует использовать необходимые механические (см. (22)–(24) и (26) в [22]) и теплофизические (см. (3.20) в [23]) граничные условия, а также начальные условия, заданные при $t=t_0$ (см. (25), (26) в [22] и (3.22) в [23]).

Для численного интегрирования поставленной нелинейной начально-краевой задачи используем явную схему шагов по времени, т. е. неизвестные функции будем определять в дискретные моменты времени ${t=t_n}$ ($n=0, 1, 2, \ldots$). При этом считаем, что в моменты времени $t_m$ уже вычислены значения следующих функций [22, 23]:
\[ \begin{equation*}
{\mathop w\limits^m} ({\bf x} )\equiv w (t_m , {\bf x} ),\quad 
{\mathop u \limits^m }{} _l^{(p)} ({\bf x})\equiv u_l^{(p)} (t_m ,{\bf x} ),\quad 
{\mathop \sigma \limits^m}{} _{ij} ( {\bf r})\equiv \sigma _{ij} (t_m ,{\bf r}), 
\end{equation*} \]
\[ \begin{equation}
\begin{aligned}
&
{\mathop \sigma \limits^m }{}_{33}^{(\pm )} ( {\bf x} ) \equiv \sigma _{33}^{(\pm )} (t_m ,{\bf x} ),
\quad 
{\mathop U \limits^n} {}^{(r)} ({\bf x} )\equiv U^{(r)} (t_n ,{\bf x} ),\quad 
{\mathop q \limits^n }_i ( {\bf r} )\equiv q_i (t_n , {\bf r} ), 
\\ 
&
{\mathop \Theta \limits^m} _s ( {\bf x} )\equiv \Theta _s (t_m , {\bf x}),\quad 
{\mathop {\dot {\Theta }} \limits^{n-1} }_s ( {\bf x} )\equiv \dot {\Theta }_s ( t_{n-1} , {\bf x}),
\quad 
{\mathop q \limits^n} {}_\infty ^{(\pm )} ( {\bf x} )\equiv q_\infty ^{(\pm )} (t_n , {\bf x}), 
\\ 
& 
{\mathop \sigma \limits^m }{}_{ij}^{(k)} ( {\bf r} )\equiv \sigma _{ij}^{(k)} (t_m , {\bf r} ),
\quad 
{\mathop \varepsilon \limits^m}{}_{ij}^{(k)} ({\bf r})\equiv \varepsilon _{ij}^{(k)} (t_m ,{\bf r} ),
\quad 
{\mathop \chi \limits^m }{}^{(k)} ( {\bf r}) \equiv \chi ^{(k)} (t_m ,{\bf r} ), \\ 
& 
l=1,2,\quad i,j=\overline {1,3},\quad m=n-1,n,\quad 0\leqslant p \leqslant M+1,\quad 0\leqslant r\leqslant L-2, \\ 
& 0\leqslant s\leqslant L,\quad 0\leqslant k\leqslant K,\quad {\bf x} \in G,\quad | x_3 |\le h, 
\end{aligned}
\end{equation} \tag{15} \]
где 
\[ \begin{equation}
\begin{aligned}
& u_l^{(p)} (t,{\bf x} ) \equiv \int _{-h}^h U_l (t, {\bf r} )x_3^p dx_3,
\quad 
U^{(r)} (t, {\bf x} )\equiv \int _{-h}^h U (t,{\bf r} ) x_3^r dx_3, \\ 
& l=1,2,\quad 0\leqslant p\leqslant M+1,\quad 0\leqslant r\leqslant L-2; 
\end{aligned}
\end{equation} \tag{16} \]
$U$ — удельная внутренняя энергия композиции оболочки; $\sigma _{ij}^{(k)}$, $\varepsilon _{ij}^{(k)}$ — напряжения и деформации в $k$-й фазе композиции ($k=0$ — связующий материал, $k\geqslant 1$ — арматура $k$-го семейства); $\chi ^{(k)}$ — параметр Одквиста в той же фазе; $q_i$ — компоненты вектора теплового потока в композиции; $q_\infty ^{(\pm )}$ — заданные значения тепловых потоков через нижнюю $(-)$ и верхнюю $(+)$ лицевые поверхности оболочки. Искомые функции $w$, $u_i$ и $\varepsilon _{i3}^{(m)} $ ($i=1,2$, $0\leqslant m\leqslant M$) в соотношениях (2) и (3) однозначно определяются через кинематические переменные $u_l^{(p)}$ (см. выражения (16)) за счет матричного равенства (30) в [22].

В механической составляющей рассматриваемой связанной неизотермической упругопластической задачи производные по времени $t$ аппроксимируем их конечно разностными аналогами на трехточечном шаблоне $\{ t_{n-1}, t_n , t_{n+1} \}$. Это позволяет построить явную численную схему [22]. Так, после замены в двумерных уточненных уравнениях динамического равновесия гибкой пологой КМ-оболочки вторых производных по $t$ от кинематических переменных $w$ и $u_i^{(p)}$ их конечными разностями при учете соотношений (3), (4), (16) и обозначений, аналогичных (15), будем иметь [22] 
\[ \begin{equation}
\begin{aligned}
&
\frac{2h\rho }{\tau ^2} \bigl(\mathop w\limits^{n+1} - 2\mathop w\limits^n + \mathop w\limits^{n-1} \bigr)=
\sum\limits_{j=1}^2 \partial _j \biggl( 
{\mathop M \limits^n} {}_{j3}^{(0)} + 
\sum\limits_{i=1}^2 {\mathop M \limits^n} {}_{ji}^{(0)} \partial _i \mathop w\limits^n \biggr) - \sum\limits_{i=1}^2 R_i^{-1} \mathop M \limits^n{}_{ii}^{(0)} 
+ \mathop \sigma \limits^n{} _{33}^{(+)} -\mathop \sigma \limits^n {} _{33}^{(-)}, 
\\ 
&
\frac{\rho }{\tau ^2}\bigl( 
\mathop u\limits^{n+1}{}_i^{(l)} - 
2\mathop u\limits^n{}_i^{(l)} +
\mathop u\limits^{n-1}{}_i^{(l)} \bigr) = 
\sum\limits_{j=1}^2 \partial_j \bigl(
\mathop M \limits^n{}_{ij}^{(l)} -
\mathop M \limits^n{}_{j3}^{(l)} \partial _i \mathop w\limits^n 
\bigr)-h^l\bigl[ 
\mathop \sigma \limits^n{} _{33}^{(+)} -(-1 )^l \mathop \sigma \limits^n {}_{33}^{(-)} \bigr] 
\partial_i \mathop w \limits^n - 
l\mathop M \limits^n {}_{i3}^{(l-1)}+ 
l\mathop M \limits^n {}_{33}^{(l-1)} 
\partial_i \mathop w\limits^n +R_i^{-1} \mathop M \limits^n {}_{i3}^{(l)} 
,\\
& {\bf x}\in G,\quad i=1,2,\quad 0\leqslant l\leqslant M+1,\quad n=1,2,3,\ldots, 
\end{aligned}
\end{equation} \tag{17} \]
где 
\[ \begin{equation}
\begin{aligned}
& 
\rho \equiv \rho _0 \omega _0 +\sum\limits_{k=1}^K \rho_k \omega _k,\quad 
\omega _0 =1-\sum\limits_{k=1}^K \omega _k, \quad 
M_{ij}^{(l)} (t, {\bf x} )\equiv \int _{-h}^h \sigma_{ij} (t, {\bf r} )x_3^l dx_3, 
\\ 
& 
lM_{33}^{(l-1)} (t, {\bf x} )\equiv l\int _{-h}^h \sigma _{33} (t, {\bf r} )x_3^{l-1} dx_3 = 
\frac{h^l}{2}
\Bigl[ (\sigma _{33}^{(+)} + \sigma _{33}^{(-)} ) \bigl(1- (-1 )^l \bigr)+\frac{l}{l+1} (\sigma _{33}^{(+)} -\sigma _{33}^{(-)} ) \bigl(1+ (-1 )^l\bigr) 
\Bigr],\quad i,j=\overline {1,3},\quad 0\leqslant l\leqslant M+1; 
\end{aligned}
\end{equation} \tag{18} \]
$\rho _0$, $\rho _k$ — объемные плотности связующей матрицы и волокон $k$-го семейства; $\tau =\rm{const}>0$ — шаг по времени. Внешние объемные силы в равенствах (17) не учитываются. Осредненные напряжения $\sigma _{ij}$ в выражениях (18) связаны с напряжениями в фазах композиции $\sigma _{ij}^{(k)}$ ($i,j=\overline {1,3}$, $0\leqslant k\leqslant K$) структурными соотношениями из [23] (по правилу простой смеси).

По формулам (18) с учетом предположений (15) в текущий момент времени $t_n$ можно определить силовые факторы $M_{ij}^{(l)}$ и внешние силы $\sigma _{33}^{(\pm )}$ в правых частях уравнений (17), поэтому, учитывая необходимые граничные условия [22], из равенств (17) по явной схеме можем вычислить искомые функции $\mathop w\limits^{n+1}$ и 
$\mathop u \limits^{n+1}{\!}_i^{(l)}$ при $t=t_{n+1}$.

Для численного интегрирования теплофизической составляющей исследуемой задачи также применим явную схему шагов по времени, но на двухточечном шаблоне $ \{ t_n , t_{n+1} \}$. При этом с учетом выражений (14), (16) и обозначений типа (15) двумерные уравнения теплового баланса пологой КМ-оболочки примут вид [23]:
\[ \begin{equation}
\begin{aligned}
&
\frac{\rho }{\tau } (
\mathop U \limits^{n+1}{}^{(m)} -
\mathop U \limits^n{}^{(m)} )=-\partial _1 \mathop Q \limits^n {}_1^{(m)} 
-\partial _2 \mathop Q \limits^n {}_2^{(m)} -
\mathop {\bar {Q}}\limits^n {}_3^{(m)} +
\mathop W \limits^n{}^{(m)} , \\ 
& {\bf x} \in G,\quad 0\leqslant m\leqslant L-2,\quad n=0, 1, 2,\ldots; 
\end{aligned}
\end{equation} \tag{19} \]
\[ \begin{equation} 
\begin{aligned}
& 
-\sum\limits_{l=0}^L (-1 )^lh^{l-1} (l\lambda _{33}^{(-)} +h\alpha ^{(-)} )\Theta _l (t,{\bf x} )=\alpha ^{(-)} (\Theta _\infty ^{(-)} -\Theta ^0 )+q_\infty ^{(-)} (t, {\bf x} ), 
\\ 
&
\sum\limits_{l=0}^L h^{l-1} (l\lambda _{33}^{(+)} + h\alpha ^{(+)} )\Theta _l (t, {\bf x} )=\alpha ^{(+)} ( \Theta _\infty ^{(+)} -\Theta ^0 ) - q_\infty ^{(+)} (t, {\rm {\bf x}} ),\quad {\bf x}\in G,\; t\geqslant t_0 ; 
\end{aligned}
\end{equation} \tag{20} \]
\[ \begin{equation}
C_0 \sum\limits_{i=0}^L H (i+m )\Theta _i + 
 \frac{C_1 }{2}\sum\limits_{i=0}^L \sum\limits_{j=0}^L H \ (i+j+m )\Theta _i \Theta _j +\frac{C_2}{3}\sum\limits_{i=0}^L \sum\limits_{j=0}^L \sum\limits_{l=0}^L H (i+j+l+m )
 \Theta _i \Theta _j \Theta _l = {} \\
 {}= U^{(m)} (t, {\bf x} ), \quad {\bf x} \in G,\; t\geqslant t_0,\; 0\leqslant m\leqslant L-2, 
\end{equation} \tag{21} \]

где
\[ \begin{equation}
\begin{aligned}
& H ( s )\equiv \frac{h^{s+1}}{s+1} [ 1- (-1 )^{s+1}],\quad
 Q_i^{(m)} (t,{\bf x} ) \equiv \int_{-h}^h q_i (t, {\bf r} ) x_3^m dx_3 \quad (i=\overline {1,3}), 
 \\ 
& \bar {Q}_3^{(m)} (t, {\bf x} )\equiv \int _{-h}^h \partial _3 q_3 (t, {\bf r} ) x_3^m dx_3 =
h^m [q_3^{(+)} - (-1 )^mq_3^{(-)} ]-mQ_3^{(m-1)} ( t, {\bf x} ), \\ 
& W^{(m)} (t, {\bf x} )\equiv \int _{-h}^h \sigma _{ij} \dot {\varepsilon }_{ij} x_3^m dx_3,
\quad
C_l ( {\bf x} )\equiv \frac{1}{\rho }\sum\limits_{k=0}^K c_l^{(k)} \rho _k \omega _k ( {\bf x} ) \quad (l=0,1,2), \\ 
&\lambda _{33}^{(\pm )} \equiv \lambda _{33} \big|_{\Theta =\Theta (t, {\bf x},\pm h )},\quad q_3^{(\pm 
)} \equiv q_\infty ^{(\pm )}; 
\end{aligned}
\!\!\!
\end{equation} \tag{22} \]
$\lambda _{33}$ — коэффициент теплопроводности композиции в поперечном направлении $Ox_3$, определяемый по структурным формулам из [28]; $\alpha ^{(\pm )}$ — коэффициент теплоотдачи на верхней $(+)$ и нижней $(-)$ лицевой поверхности оболочки; $\Theta _\infty ^{(\pm )} $ — температура окружающей среды со стороны той же поверхности; $c_l^{(k)}$ — коэффициент разложения удельной теплоемкости материала $k$-й фазы композиции $c_k (\Theta -\Theta ^0 )$ по формуле (в случае термочувствительности):
\[ \begin{equation}
c_k (\Theta -\Theta ^0 )=c_0^{(k)} +c_1^{(k)} (\Theta -\Theta ^0 )+
c_2^{(k)} (\Theta -\Theta ^0 )^2,\quad 0\leqslant k\leqslant K.
\end{equation} \tag{23} \]

Равенства (20) — тепловые граничные условия на лицевых поверхностях искривленной КМ-панели, преобразованные с учетом представления температуры в виде (14). Соотношение (21) выражает двумерные функции $U^{(m)}$ (см. (16)) через коэффициенты разложения (14) с учетом равенства (23). Уравнения (20), (21) справедливы для любого момента времени $t$.

Используя соотношения (22) и предположения (15), в текущий момент времени $t_n$ можно вычислить правые части в равенствах (19). Учитывая при этом необходимые граничные (заданные на контуре $\Gamma$, ограничивающем область $G$) и начальные тепловые условия [23], получаем из (19) по явной схеме значения $\mathop U\limits^{n+1}{\!}^{(m)}$ при $t=t_{n+1}$. Далее, в следующий момент времени $t_{n+1}$ из равенств (20), (21) (где правые части уже известны) при учете (22) определяем коэффициенты разложения для температуры $\mathop \Theta \limits^{n+1} {}_l ( {\bf x} )$, $0\leqslant l\leqslant L$ (см. (14)). В случае термочувствительности компонентов композиции ($c_1^{(k)} \ne 0$ и/или $c_2^{(k)} \ne 0$ ($0\leqslant k\leqslant K$) в равенстве (23)) система (20), (21) является нелинейной. Для ее линеаризации может быть использован метод переменных теплофизических параметров, аналогичный методу переменных параметров упругости. В остальных деталях разработанный численный метод (17), (19) при учете соотношений (2), (5)–(13), (18), (22) и структурных формул, полученных в [28] (см. выражения (20) и (22)), реализуется так же, как подробно описано в [22, 23]. При этом в соотношениях (5), (7)–(9), (11) и (12) всем величинам следует приписать верхний индекс $n$ (аналогично (17) и (19)), означающий номер текущего дискретного момента времени $t_n$. (Величины же, входящие в равенства (10), от $n$ не зависят.) Кроме того, в соотношениях (7) и (8) скорость температуры $\dot {\Theta }$ необходимо выразить разностной формулой, вытекающей из формулы трапеций [23]:
\[ \begin{equation}
\dot {\Theta }\equiv \mathop {\dot {\Theta }}\limits^n =
\frac{2}{\tau }\bigl( \mathop \Theta \limits^n - \mathop \Theta \limits^{n-1/2} \bigr),
\quad 
\mathop \Theta \limits^{n-1/2} \equiv \mathop \Theta \limits^{n-1} +\frac{\tau }{2} \mathop {\dot {\Theta }}\limits^{n-1},
\quad n=1, 2, 3,\dots .
\end{equation} \tag{24} \]
Согласно предположениям (15), при $t=t_n$ правые части в равенствах (24) известны, т. е. известна и скорость температуры $\dot {\Theta }=\dot {\Theta }\left( {t_n ,{\rm {\bf r}}} \right)$.

В работе [23] показано, что в уравнениях (17), (19) шаг по времени $\tau$ нужно задавать так, чтобы выполнялось необходимое условие устойчивости Куранта – Фридрихса – Леви [5].

2. Обсуждение результатов расчетов

Рассматривается термоупругопластический динамический изгиб относительно тонких пологих цилиндрических КМ-оболочек толщиной $2h=2$ см, имеющих в плане удлиненную прямоугольную форму ($G$: $ |x_1 |\leqslant a$, $ |x_2 |\leqslant b$, $a=3b$; $1/R_1 \equiv 0$, $R_2 \equiv R=\rm{const}$; $b=50$ см). Стрела подъема искривленных панелей $f$ над продольными (в направлении оси $Ox_1$) кромками равна 10 см. Радиус цилиндрической оболочки $R$ при этом определяется формулой [22]: $R= (b^2+f^2 )/ (2f )$, $0\leqslant f<0.4b$. По кромкам конструкции имеют жесткое закрепление: $w=0$, $u_i^{(m)} =0$, ${\bf x}\in \Gamma $, $t\geqslant t_0$ (см. (16), (17)). До начального момента времени $t=t_0=0$ панели покоятся ($w=0$, $u_i^{(m)} =0$, ${\bf x}\in G$, $t=t_0$, $i=1, 2$, $0\leqslant m\leqslant M+1$) в естественном состоянии при $\Theta =\Theta ^0=20$°C (${\bf x}\in G$, $ |x_3 |\leqslant h$, $t=t_0$). Начиная с момента времени $t_0$ оболочки подвергаются воздействию избыточного давления $p( t )$, порожденного воздушной взрывной волной [26]:
\[ \begin{equation}
\begin{aligned}
& \sigma _{33}^{(+)} -\sigma _{33}^{(-)} \equiv p ( t )=
\begin{cases}
 p_{\max} \cdot t/t_{\max},& 0\leqslant t\leqslant t_{\max}, \\ 
 p_{\max} \cdot \exp [-\alpha (t-t_{\max} ) ], & t>t_{\max}, 
\end{cases} 
\\ 
&\alpha =-\ln (0.01)/ (t_{\min} -t_{\max} )>0,\quad t_{\min} \gg t_{\max},\quad p (t_{\min} )=0.01p_{\max}, 
\end{aligned}
\end{equation} \tag{25} \]
где смысл величин $p_{\max}$, $t_{\max}$, $t_{\min}$ и $\alpha$ вполне очевиден и описан в [22, 23]. 
В соответствии с экспериментальными данными [26] в расчетах задано $t_{\max}=0.1$ мс и $t_{\min}=2$ мс. Из соотношений (25) получаем: при $p_{\max}>0$ конструкции нагружены со стороны нижней (вогнутой) лицевой поверхности (при этом условно считаем $\sigma _{33}^{(+)} \equiv 0$), а при $p_{\max} <0$ — со стороны верхней (выпуклой) поверхности (при этом принимаем $\sigma _{33}^{(-)} \equiv 0$).

Теплообмен оболочек с окружающей средой через их лицевые поверхности реализуется в условиях естественной конвекции ($\alpha^{(\pm )}=30$ Вт/(м$^2{}\cdot{}$K [29], $q_\infty^{(\pm)} \equiv 0$) при температуре воздуха, равной температуре естественного состояния конструкций: $\Theta _\infty ^{(\pm )} =\Theta ^0$ (см. (20)). На торцевых поверхностях цилиндрических панелей задана температура, которая поддерживается равной $\Theta^0$.

Оболочки выполнены из эпоксисвязующего [17], армированного стеклянными волокнами [16] (стеклопластик), или магниевого сплава ВТ65 [30], усиленного стальной проволокой У8А [16] (металлокомпозиция). Упругопластическое поведение компонентов композиции при постоянной температуре $\Theta$ и активном нагружении определяется билинейной диаграммой растяжения – сжатия:
\[ \begin{equation*}
\sigma =\begin{cases}
\qquad\qquad \quad 
E^{(k)}\varepsilon, & |\varepsilon |\leqslant \varepsilon _s^{(k)} =\sigma _s^{(k)} /E^{(k)}, 
\\ 
\operatorname{sign}(\varepsilon)\sigma_s^{(k)} +E_s^{(k)} 
\bigl(\varepsilon -\operatorname{sign} (\varepsilon)\varepsilon_s^{(k)} \bigr), & | \varepsilon |
>\varepsilon_s^{(k)},\; 0\leqslant k\leqslant K, 
\end{cases}
\end{equation*} \]
где $\sigma $, $\varepsilon $ — осевое напряжение и соответствующая ему деформация; $E^{(k)} =E^{(k)}( \Theta )$, $E_s^{(k)}=E_s^{(k)} ( \Theta)$ — модули Юнга и линейного упрочнения материала $k$-й фазы композиции; $\sigma _s^{(k)} =\sigma _s^{(k)} ( \Theta)$ — предел текучести этого же компонента. Физико-механические характеристики материалов фаз композиций представлены в таблице, где $\nu $ — коэффициент Пуассона, $\lambda $ — коэффициент теплопроводности, $\alpha $ — коэффициент линейного теплового расширения, $c$ — удельная теплоемкость, а в скобках приведено значение температуры ($\Theta$, °C), при которой была определена данная характеристика. В расчетах использованы линейные аппроксимации по $\Theta$ для всех физико-механических характеристик по данным, указанным в таблице.

Физико-механические характеристики материалов композиций [16, 17, 30]
[Physico-mechanical characteristics of the components of composite [16, 17, 30]

Характеристика материала [Material characteristic]	Эпоксисвязующее [Epoxy]	Стекловолокно [Fiberglass]	Магниевый сплав ВТ65(Mg) [VT65(Mg) magnesium alloy]	Стальная проволока У8А [U8A steel wire]
$\rho$, kg/m$^3$	1210.0 (20°C)	2520.0 (20°C)	1800.0 (20°C)	7800.0 (20°C)
	1208.0 (40°C)	2519.6 (80°C)	1796.2 (100°C)	7791.8 (100°C)
$E$, GPa	2.8 (20°C)	86.8 (20°C)	43.0 (20°C)	210.0 (20°C)
	2.3 (40°C)	86.3 (80°C)	38.5 (100°C)	195.0 (100°C)
$\nu$	0.33 (20°C)	0.25 (20°C)	0.330 (20°C)	0.3 (20°C)
	0.333 (40°C)	0.254 (80°C)	0.334 (100°C)	0.305 (100°C)
$\sigma_s$, MPa	20 (20°C)	4500 (20°C)	267 (20°C)	3968 (20°C)
	18 (40°C)	4400 (80°C)	219 (100°C)	3971 (200°C)
$E_s$, GPa	1.114 (20°C)	6.230 (20°C)	0.379 (20°C)	6.973 (20°C)
	0.783 (40°C)	5.168 (80°C)	0.367 (100°C)	5.014 (200°C)
$\lambda$, W/(m${}\cdot{}$K)	0.243 (20°C)	0.89 (20°C)	117.23 (20°C)	42.7 (20°C)
	0.240 (40°C)	0.86 (80°C)	121.42 (100°C)	41.7 (100°C)
$\alpha\cdot 10^6$, K$^{-1}$	68.1 (20°C)	2.5 (20°C)	20.9 (20°C)	12.3 (20°C)
	70.3 (40°C)	2.6 (80°C)	22.6 (100°C)	13.2 (100°C)
$c$, kJ/(kg${}\cdot{}$K)	1.54 (20°C)	0.800 (20°C)	1.032 (20°C)	0.485 (20°C)
	1.60 (40°C)	0.839 (80°C)	1.054 (100°C)	0.488 (100°C)

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

Пологие оболочки армированы двумя ($K=2$) семействами непрерывных волокон, ориентированных по направлениям $Ox_1$ и $Ox_2$ (см. рис. 1), с плотностями $\omega _1 =0.1$ и $\omega _2 =0.3$ соответственно. Углы армирования при этом имеют следующие значения (см. рис. 2): $\theta _1 =\theta _2 =\pi /2$, $\varphi _1 =0$ и $\varphi _2 =\pi /2$.

В [23] были проведены расчеты для гибких КМ-пластин с рассматриваемой структурой армирования при тех же характерных размерах тонкостенных конструкций ($h$, $a$, $b$ и $f=0$). При этом в соотношениях (2) и (3) принималось значение $M=0$ (использовалась теория изгиба Амбарцумяна [1]), а в разложении температуры (14) варьировалось значение $L\geqslant 2$. Было показано, что в динамически изгибаемых КМ-пластинах температурные поля следует определять по формуле (14) при задании в них $L=6$ или $L=7$.

В работе же [22] было продемонстрировано, что в случае изотермического изгибного динамического деформирования пологих КМ-оболочек для адекватного описания их неупругого механического поведения в кинематических соотношениях (2) и (3) нужно принимать не $M=0$, а $M=6$ или $M=7$ (т. е. следует применять уточненную теорию изгиба).

Как уже отмечалось ранее, в данной работе исследуется влияние использования уточненных соотношений (2) и (3) (при $M=7$) на результаты расчета термомеханического отклика динамически изгибаемых искривленных КМ-панелей при упругопластическом деформировании материалов компонентов их композиций.

В связи с этим на рис. 3 изображены зависимости максимальных значений температуры $\Theta_m (t;M)=\mathop {\max }\limits_{\bf r} \Theta (t, {\bf r}; M )$ в рассматриваемых пологих КМ-оболочках. При этом на основании результатов, полученных в [23], в формуле (14) было принято $L=7$. Кривые 1 и 2 рассчитаны соответственно при $M=0$ (теория Амбарцумяна) и ${M=7}$ (уточненная теория изгиба) для стеклопластиковой панели, нагруженной снизу (рис. 3, a) при $p_{\max} =8$ МПа (см. (25)) и сверху (рис. 3, b) при $p_{\max } =-8$ МПа, а также металлокомпозитной конструкции, нагруженной сверху при $p_{\max } =-30$ МПа (рис. 3, c). Для случая нагружения металлокомпозитной оболочки снизу ($p_{\max} =30$ МПа) зависимости $\Theta _m (t;M)$ здесь не приводятся, так как кривые 1 и 2 при этом на рассматриваемом интервале времени $0\leqslant t\leqslant 24$ мс практически совпадают (незначительное различие наблюдается лишь при $t>22$ мс).

Рис. 3. Зависимость от времени максимального значения температуры в стеклопластиковой пологой оболочке, нагруженной снизу (a) и сверху (b), а также в металлокомпозитной панели, нагруженной сверху (c)
[Figure 3. Time dependence of the maximum temperature value in a fiberglass shallow shell loaded from below (a) and from above (b), as well as in a metal-composite panel loaded from above (c)]

Сравнение кривых 1 и 2 на рис. 3 показывает, что температурные поля, рассчитанные по простейшей неклассической теории изгиба искривленной КМ-панели — теории Амбарцумяна (см. кривые 1) — и по уточненной теории изгиба (см. кривые 2), с течением времени все более различаются. Следовательно, есть основания предполагать, что теория Амбарцумяна, широко используемая в расчетной практике, в случае динамического термоупругопластического деформирования тонкостенных КМ-конструкций типа оболочек может оказаться непригодной для проведения таких расчетов. Окончательное же решение этого вопроса возможно только после сравнения с результатами соответствующих экспериментов, которые автору неизвестны. (Здесь следует отметить, что дополнительные расчеты показали: в случаях квазистатического нагружения и динамического термоупругого деформирования пологих КМ-оболочек результаты, полученные по теории Амбарцумяна и по уточненной теории изгиба, практически совпадают.)

Сделанный выше вывод подтверждают и расчеты поперечных колебаний рассматриваемых панелей. Так, на рис. 4 представлены осцилляции прогиба центральной точки стеклопластиковой пологой КМ-оболочки ($w_0 ( t ) \equiv w (t, 0, 0)$), определенные при $p_{\max } =8$ МПа (рис. 4, a) и $p_{\max } =-8$ МПа (рис. 4, b). Кривые 1 и 2 на рис. 4 получены при тех же значениях $M$, что и на рис. 3, а кривые 2$'$ на рис. 4 изображены для сравнения: они рассчитаны по уточненной теории изгиба и соответствуют изотермическому упругопластическому случаю, когда изменение температуры в пологой КМ-оболочке не учитывается. Поведение кривых 1 и 2 на рис. 4 свидетельствует о том, что с увеличением времени зависимости $w_0 (t)$, рассчитанные по теории Амбарцумяна (кривые 1) и по уточненной теории (кривые 2), существенно различаются. Аналогичное сопоставление кривых 2 и 2$'$ демонстрирует, что с увеличением $t$ зависимости $w_0 ( t )$, определенные с учетом (кривые 2) и без учета (кривые 2$'$) температурного отклика в искривленной стеклопластиковой панели, все более отличаются друг от друга. Согласно этим результатам, в еще большей степени различаются напряженно-деформированные состояния в компонентах композиции таких тонкостенных конструкций, рассчитанные по теории Амбарцумяна ($M=0$) и уточненной теории изгиба ($M=7$), а также с учетом и без учета температурных полей, возникающих в них. Соответствующие зависимости здесь не представлены.

Рис. 4. Осцилляции прогиба центральной точки стеклопластиковой пологой оболочки, нагруженной снизу (a) и сверху (b), рассчитанные по разным теориям изгиба
[Figure 4. Deflection oscillations of the central point of a fiberglass shallow shell loaded from below (a) and from above (b) calculated using different bending theories]

На рис. 5 изображены аналогичные кривые 1, 2 и 2$'$, рассчитанные для металлокомпозитной пологой оболочки при $p_{\max } =30$ МПа (рис. 5, a) и $p_{\max } =-30$ МПа (рис. 5, b). Как видно из рис. 5, b, при нагружении искривленной Mg–У8А-панели сверху уже на рассматриваемом малом интервале времени ($0\leqslant t\leqslant 24$ мс) кривые 1 и 2, а также 2 и 2$'$ значительно различаются между собой. Кривые 1 и 2 на рис. 5, a (нагружение снизу) визуально незначительно различаются между собой лишь при $t \approx 23.5$ мс. Объяснятся это тем, что, как отмечалось выше, при нагружении металлокомпозитной пологой оболочки снизу зависимости $\Theta_m (t;M)$, рассчитанные по теории Амбарцумяна (кривая 1) и по уточненной теории изгиба (кривая 2), заметно различаются только при $t>22$ мс. Однако сравнение кривых 2 и 2$'$ на рис. 5, a свидетельствует о том, что даже в этом случае нагружения расчеты Mg–У8А-панели, выполненные по уточненной теории изгиба с учетом (кривая 2) и без учета (кривая 2$'$) теплового отклика в ней, приводят к значительно различающимся зависимостям $w_0 ( t )$.

Несмотря на то, что на рис. 5, a кривые 1 и 2 визуально почти не различаются, деформированные состояния компонентов композиции, рассчитанные по теории Амбарцумяна и по уточненной теории изгиба, в этом случае на рассматриваемом интервале времени различаются весьма заметно. Так, на рис. 6 изображены осцилляции максимальных значений интенсивности деформаций $\varepsilon _\ast ^{(k)}$ ($\varepsilon_m^{(k)} (t;M)=\mathop {\max }\limits_{\bf r} \varepsilon_\ast ^{(k)} (t, {\bf r}; M )$) связующего материала (рис. 6, a; $k=0$) и арматуры второго семейства (рис. 6, b; $k=2$), которая испытывает наибольшее деформирование в металлокомпозитной панели, нагруженной снизу. Кривые 1 и 2 на рис. 6 получены при тех же условиях, что и аналогичные кривые на рис. 5, a. Кривые 1 и 2 на рис. 6 (особенно на рис. 6, b) показывают, что расчет деформированного состояния материалов фаз композиции даже в этом случае необходимо проводить по уточненной теории изгиба (кривые 2). Действительно, поведение кривой 1 на рис. 6, b (расчет по теории Амбарцумяна) свидетельствует о том, что наибольшее значение интенсивности деформаций в волокнах второго семейства $\varepsilon _{\max}^{(2)} =\mathop {\max }\limits_t \varepsilon _m^{(2)} ( t )$ достигается при $t=19.1$ мс и равно 4.99 %, а согласно поведению кривой 2 на этом рисунке (уточненная теория изгиба), значение $\varepsilon _{\max}^{(2)} =4.78$ % достигается при $t=6.7$ мс.

Рис. 5. Осцилляции прогиба центральной точки металлокомпозитной пологой оболочки, нагруженной снизу (a) и сверху (b), рассчитанные по разным теориям изгиба
[Figure 5. Deflection oscillations of the central point of a metal-composite shallow shell loaded from below (a) and from above (b) calculated using different bending theories]

Согласно поведению кривой 2 на рис. 3, c, температура в металлокомпозитной панели, нагруженной сверху, в отдельных точках в некоторые моменты времени может достигать порядка 150–170°C. При нагружении же такой пологой КМ-оболочки снизу температура не превышает 70°C. (Как уже отмечалось, соответствующие зависимости $\Theta _m (t;M)$ здесь не приводятся.) Целесообразно предположить, что высокий уровень нагрева динамически изгибаемой металлокомпозитной конструкции обусловлен тем, что в проведенных расчетах предполагалось: теплообмен КМ-панели с окружающей средой через лицевые поверхности реализуется при естественной конвекции воздуха ($\alpha^{(\pm )}=30$ Вт/(м$^2{}\cdot{}$K)). Однако при нагружении искривленной панели взрывной волной она осциллирует в поперечном направлении с высокой частотой (см. кривые на рис. 5), поэтому теплообмен через ее лицевые поверхности следовало бы моделировать не при естественной, а при вынужденной конвекции. Согласно этому допущению, были выполнены дополнительные расчеты при значении коэффициентов теплоотдачи $\alpha^{(\pm )}=300$ Вт/(м$^2{}\cdot{}$K) (см. равенства (20)), которое характеризует вынужденную конвекцию для газов [29]. Определенные таким образом значения $\Theta_m (t;M)$ отличаются от полученных ранее, т. е. при естественной конвекции (см. рис. 3), лишь в четвертой значащей цифре. А значит, увеличение коэффициентов теплоотдачи на лицевых поверхностях изучаемых искривленных КМ-панелей в 10 раз практически не оказывает влияния на результаты расчетов их неизотермического упругопластического поведения на рассматриваемых временных интервалах ($0\leqslant t \leqslant 50$ мс). Следовательно, на таких, достаточно малых, интервалах времени теплофизические процессы в пологих КМ-оболочках можно считать адиабатическими, т. е. можно не учитывать кондуктивную составляющую процесса теплопереноса в них.

Рис. 6. Осцилляции максимальных значений интенсивности деформаций связующего материала (a) и волокон второго семейства (b) в металлокомпозитной пологой оболочке, нагруженной снизу, рассчитанные по разным теориям изгиба
[Figure 6. Oscillations of the maximum values of the intensity of strains of the binder material (a) and fibers of the second family (b) in a metal-composite shallow shell loaded from below, calculated according to different bending theories]

Заключение

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

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

Выполненные расчеты показали, что при нагружении поперечной взрывной волной относительно тонкие искривленные панели из стеклопластика в отдельных точках кратковременно (менее 1 мс) могут дополнительно нагреваться на 14–34°C (аналогичные по структуре и характерным размерам стеклопластиковые пластины нагреваются не более чем на 14°C [23]), а металлокомпозитные пологие оболочки в отдельных точках способны дополнительно нагреваться на 50°C при их нагружении со стороны вогнутой лицевой поверхности и на 150°C при нагружении со стороны выпуклой лицевой поверхности (аналогичные металлокомпозитные гибкие пластины нагреваются не выше, чем на 30–32°C [23]). Более интенсивное тепловыделение и, как результат, более активное воздействие температурных полей на механический отклик всегда наблюдается при динамическом нагружении искривленных КМ-панелей со стороны их выпуклой лицевой поверхности.

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

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

About the authors

Andrey P. Yankovskii

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

Author for correspondence.
Email: lab4nemir@rambler.ru
ORCID iD: 0000-0002-2602-8357
SPIN-code: 9972-3050
Scopus Author ID: 7003288442
http://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