Using RBF-FD for calculation of hydroelastic vibrations of axisymmetric orthotropic shells of rotation

Толық мәтін

Осесимметричными оболочками вращения моделируют топливные баки ракет-носителей с жидкостным ракетным двигателем при решении задачи обеспечения продольной устойчивости (англ. POGO) [1]. Наиболее популярным и универсальным методом анализа тонкостенных конструкций является метод конечных элементов (МКЭ), для учета жидкости часто используют метод граничных элементов (МГЭ). Целью данной работы является развитие и программная реализация альтернативного МКЭ и МГЭ подхода к расчету тонкостенных оболочечных конструкций, взаимодействующих с жидкостью.

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

В трехмерном пространстве рассмотрим криволинейную поверхность оболочки, заданную радиус-вектором:

$\begin{array}{l}\overrightarrow{r}\left(s_1,s_2\right)=x_1\left(s_1,s_2\right)\overrightarrow{i_1}+x_2\left(s_1,s_2\right)\overrightarrow{i_2}+\\ +x_3\left(s_1,s_2\right)\overrightarrow{i_3}=x_k\left(s_1,s_2\right)\overrightarrow{i_k},\end{array}$ (1)

где $\overrightarrow{i_k}$ – базисные векторы глобальной системы координат, по повторяющимся индексам ведется суммирование от 1 до 3. Поверхность разделяется координатными линиями, бесконечно малый элемент, образуемый линиями s₁, s₂, s₁ + ds₁, s₂ + ds₂, в точке O этого элемента поставим локальную систему координат с тремя единичными векторами ${\overrightarrow{e}}_1$, ${\overrightarrow{e}}_2$, ${\overrightarrow{e}}_3$ (рис. 1).

Рис. 1. Деформирование и равновесие малого элемента оболочки.

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

${\overrightarrow{e}}_1=\frac{\partial \overrightarrow{r}}{\partial s_1}, {\overrightarrow{e}}_2=\frac{\partial \overrightarrow{r}}{\partial s_2}, {\overrightarrow{e}}_3={\overrightarrow{e}}_1\times {\overrightarrow{e}}_2;$  

$\left\{\begin{array}{c}{\overrightarrow{e}}_1\\ {\overrightarrow{e}}_2\\ {\overrightarrow{e}}_2\end{array}\right\}=\left[\begin{array}{ccc}{\beta }_{11}& {\beta }_{12}& {\beta }_{13}\\ {\beta }_{21}& {\beta }_{22}& {\beta }_{23}\\ {\beta }_{31}& {\beta }_{32}& {\beta }_{33}\end{array}\right]\left\{\begin{array}{c}{\overrightarrow{i}}_1\\ {\overrightarrow{i}}_2\\ {\overrightarrow{i}}_3\end{array}\right\}\iff {\overrightarrow{e}}_i={\beta }_{ik}{\overrightarrow{i}}_k,$ (2)

где β – матрица начальной геометрии.

Поворот малого элемента определим следующим образом:

$\left\{\begin{array}{c}{\overrightarrow{i}}_1^{*}\\ {\overrightarrow{i}}_2^{*}\\ {\overrightarrow{i}}_3^{*}\end{array}\right\}=\left[\begin{array}{ccc}{\lambda }_{11}& {\lambda }_{12}& {\lambda }_{13}\\ {\lambda }_{21}& {\lambda }_{22}& {\lambda }_{23}\\ {\lambda }_{31}& {\lambda }_{32}& {\lambda }_{33}\end{array}\right]\left\{\begin{array}{c}{\overrightarrow{i}}_1\\ {\overrightarrow{i}}_2\\ {\overrightarrow{i}}_3\end{array}\right\}\iff {\overrightarrow{i}}_n^{*}={\lambda }_{nk}{\overrightarrow{i}}_k,$

$\begin{array}{c}{\lambda }_{ii}=1-\frac{1-\cos \omega }{{\omega }^2}\left({\omega }_j^2+{\omega }_k^2\right),\\ {\lambda }_{ij}=\frac{\sin \omega }{\omega }{\omega }_k+\frac{1-\cos \omega }{{\omega }^2}{\omega }_i{\omega }_j,\\ {\lambda }_{ik}=-\frac{\sin \omega }{\omega }{\omega }_k+\frac{1-\cos \omega }{{\omega }^2}{\omega }_i{\omega }_j,\\ i=\mathrm{1,}\mathrm{2,}3;j=\mathrm{2,}\mathrm{3,}1;k=\mathrm{3,}\mathrm{1,}2;\end{array}$ (3)

где ω_k – проекции вектора конечного поворота (вектора Эйлера) на оси глобальной системы координат.

В результате перемещения параллельно самому себе как жесткого целого точка O займет новое положение O^*, длины отрезков ds₁ = | OA₁ |, ds₂ = | OA₂ | изменятся и станут равными $\left|O^{*}{A}_1^{*}\right|=d{s}_1^{*}$, $\left|O^{*}{A}_2^{*}\right|=d{s}_2^{*}$. Обозначим ${\epsilon }_1=\left(d{s}_1^{*}-d{s}_1\right)\text{/}d{s}_1$, ${\epsilon }_2=\left(d{s}_2^{*}-d{s}_2\right)\text{/}d{s}_2$ – относительные удлинения сторон малого элемента вдоль направлений s₁ и s₂ соответственно. При деформировании векторы локального базиса ${\overrightarrow{e}}_1,{\overrightarrow{e}}_2$ переходят в ${\overrightarrow{e}}_1^{*},{\overrightarrow{e}}_2^{*}$. С учетом (2) и (3) повернутые орты будут иметь следующие выражения:

${\overrightarrow{e}}_i^{*}={\beta }_{ik}{\overrightarrow{i}}_k^{*}={\beta }_{ik}{\lambda }_{kn}{\overrightarrow{i}}_n.$ (4)

При наличии сдвига в плоскости ортогональность ${\overrightarrow{e}}_1^{*}$ и ${\overrightarrow{e}}_2^{*}$ нарушается, и они переходят в ${\overrightarrow{e}}_1^{**}$, ${\overrightarrow{e}}_2^{**}$. Обозначим угол сдвига χ, тогда:

${\overrightarrow{e}}_k^{**}={\alpha }_{kj}{\overrightarrow{e}}_j^{*},$

где

$\alpha =\left(\begin{array}{cc}\cos \chi & \sin \chi \\ \sin \chi & \cos \chi \end{array}\right),$

$\begin{array}{c}d{\overrightarrow{r}}^{*}={\overrightarrow{e}}_1^{**}d{s}_1^{*}=\left(1+{\epsilon }_1\right){\alpha }_{1j}{\beta }_{jk}{\lambda }_{kn}\overrightarrow{i_n}d{s}_1,\\ d{\overrightarrow{r}}^{*}={\overrightarrow{e}}_2^{**}d{s}_2^{*}=\left(1+{\epsilon }_2\right){\alpha }_{2j}{\beta }_{jk}{\lambda }_{kn}\overrightarrow{i_n}d{s}_2,j=\mathrm{1,}2.\end{array}$ (5)

Поскольку $\overrightarrow{U}={\overrightarrow{r}}^{\ast }-\overrightarrow{r}$, то

$\begin{array}{c}\frac{\partial \overrightarrow{U}}{\partial s_1}=\left[\left(1+{\epsilon }_1\right)\left({\alpha }_{1j}{\beta }_{jk}{\lambda }_{kn}+\text{}{S}_{1n}\right)-{\beta }_{1n}\right]{\overrightarrow{i}}_n,\\ \frac{\partial \overrightarrow{U}}{\partial s_2}=\left[\left(1+{\epsilon }_2\right)\left({\alpha }_{2j}{\beta }_{jk}{\lambda }_{kn}+\text{}{S}_{2n}\right)-{\beta }_{2n}\right]{\overrightarrow{i}}_n,\end{array}$  (6)

где

$$\begin{gathered} S_{1n}=\frac{T_{1j}}{\kappa G_{13}{h}^{*}}{\lambda }_{kj}{\beta }_{3k}{\beta }_{3p}{\lambda }_{pn} \\ {S}_{2n}=\frac{T_{2j}}{\kappa G_{23}{h}^{*}}{\lambda }_{kj}{\beta }_{3k}{\beta }_{3p}{\lambda }_{pn} \end{gathered}$$ ,

суть дополнительные углы деформирования за счет поперечного сдвига, G₁₃ – модуль сдвига ортотропного материала, h^* – толщина деформированной оболочки, ê – корректирующий коэффициент сдвига прямоугольного сечения (ê ≈ π² / 12), ${\overrightarrow{T}}_1=T_{1j}{\overrightarrow{i}}_j$, ${\overrightarrow{T}}_2=T_{2j}{\overrightarrow{i}}_j$ – векторы внутренних погонных сил на площадках, ортогональных координатным линиям.

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

Для деформированной оболочки имеем $d{s}_{{}_1}^{\ast }=A_1^{\ast }d{\alpha }_1$, $d{s}_{{}_2}^{\ast }=A_2^{\ast }d{\alpha }_2$. Уравнения динамического равновесия сил (рис. 1), записанные для деформированного состояния, имеют вид

$\frac{1}{A_1^{*}{A}_2^{*}}\frac{\partial \left(A_2^{*}{\overrightarrow{T}}_1\right)}{\partial {\alpha }_1}+\frac{1}{A_1^{*}{A}_2^{*}}\frac{\partial \left(A_1^{*}{\overrightarrow{T}}_2\right)}{\partial {\alpha }_2}+\overrightarrow{q}=-{\overrightarrow{q}}_{è},$ (7)

$\begin{array}{c}\frac{1}{A_1^{*}{A}_2^{*}}\frac{\partial \left(A_2^{*}{\overrightarrow{M}}_1\right)}{\partial {\alpha }_1}+\frac{1}{A_1^{*}{A}_2^{*}}\frac{\partial \left(A_1^{*}{\overrightarrow{M}}_2\right)}{\partial {\alpha }_2}+\frac{1}{A_1^{*}}\frac{\partial {\overrightarrow{r}}^{*}}{\partial {\alpha }_1}\times {\overrightarrow{T}}_1+\\ +\frac{1}{A_2^{*}}\frac{\partial {\overrightarrow{r}}^{*}}{\partial {\alpha }_2}\times {\overrightarrow{T}}_2+\overrightarrow{m}=-{\overrightarrow{m}}_{è},\end{array}$ (8)

где ${\overrightarrow{T}}_1=T_{1j}{\overrightarrow{i}}_j$, ${\overrightarrow{T}}_2=T_{2j}{\overrightarrow{i}}_j$ – векторы внутренних погонных сил; ${\overrightarrow{M}}_1=M_{1j}{\overrightarrow{i}}_j$, ${\overrightarrow{M}}_2=M_{2j}{\overrightarrow{i}}_j$ – векторы внутренних моментов; $\overrightarrow{q}=q_j{\overrightarrow{i}}_j$ – вектор внешней распределенной нагрузки; $\overrightarrow{m}=m_j{\overrightarrow{i}}_j$ – вектор внешнего распределенного момента, ${\overrightarrow{q}}_{è}$ – сила инерции, ${\overrightarrow{m}}_{è}$ – инерционный момент. Здесь необходимо отметить, что инерционная сила и момент, в силу сохранения массы элементарного объема, зависят только от ускорения, поэтому в формуле (7) сила инерции

${\overrightarrow{q}}_{è}=-\frac{\rho h\ddot{\overrightarrow{U}}}{\left(1+{\epsilon }_1\right)\left(1+{\epsilon }_2\right)}$,

а в (8)

${\overrightarrow{m}}_{è}=-\frac{\rho h^3\left({\ddot{\Omega }}_1{\overrightarrow{e}}_1^{*}+{\ddot{\Omega }}_2{\overrightarrow{e}}_2^{*}\right)}{12\left(1+{\epsilon }_1\right)\left(1+{\epsilon }_2\right)}$,

где ${\ddot{\Omega }}_{\mathrm{1,2}}$ – угловые ускорения вокруг главный осей (точки обозначают дифференцирование по времени t).

При осесимметричном деформировании оболочки вращения (рис. 2) координатная сетка остается ортогональной, т.е. параметр сдвига χ = 0. Введем в рассмотрение цилиндрическую систему координат с ортами

$\begin{array}{c}{\overrightarrow{i}}_r=\cos \phi {\overrightarrow{i}}_1+\sin \phi {\overrightarrow{i}}_2,\\ {\overrightarrow{i}}_{\phi }=-\sin \phi {\overrightarrow{i}}_1+\cos \phi {\overrightarrow{i}}_2,\\ {\overrightarrow{i}}_z={\overrightarrow{i}}_3.\end{array}$ (9)

Радиус-вектор срединной поверхности будет иметь следующий вид:

$\overrightarrow{r}(s,\phi )=r\left(s\right){\overrightarrow{i}}_r+z\left(s\right){\overrightarrow{i}}_z,$ (10)

где r(s), z(s) – параметрические уравнения меридиана оболочки.

Рис. 2. Осесимметричная оболочка вращения с жидкостью.

Введем в рассмотрение безразмерные величины следующим образом:

$\overline{s}=\frac{s}{\mathcal{l}}, \overline{h}=\frac{h}{\mathcal{l}}, \overline{r}=\frac{r}{\mathcal{l}}, \overline{z}=\frac{z}{\mathcal{l}}, \overline{U}=\frac{U}{\mathcal{l}},$    

$\overline{T}=\frac{\left(1-{\nu }^2\right)}{E\overline{h}\mathcal{l}}T, \overline{M}=\frac{12\left(1-{\nu }^2\right)}{E{\overline{h}}^2{\mathcal{l}}^2}M,$  (11)

где $\mathcal{l}$ – характерный геометрический размер оболочки (например, длина меридиана, радиус), $\nu =\sqrt{{\nu }_{12}{\nu }_{21}}$, $E=\sqrt{E_1{E}_2}$. Введем безразмерный параметр длины меридиана ξ, т.е. s = s(ξ). Определим параметры Ламе следующими соотношениями: A₁ = r, A₂ = s_,ξ, и dα₁ = dφ, dα₂ = dξ, здесь и далее нижний индекс после запятой означает дифференцирование по этой переменной.

Пропустим промежуточные выкладки. Разрешающая система 6 геометрически нелинейных дифференциальных уравнений первого порядка относительно 6 функций – ${\overline{U}}_r,{\overline{U}}_z,\gamma ,{\overline{T}}_{2r},{\overline{T}}_{2z},{\overline{M}}_{2\phi }$, описывающих динамическое напряженно-деформированное состояние осесимметричной ортотропной оболочки при больших продольных деформациях, перемещениях и поворотах, с учетом поперечного сдвига имеет следующий вид:

1) $\frac{d{\overline{U}}_r}{d\xi }=\left(1+{\epsilon }_2\right)\left({\overline{d}}_r+{\overline{d}}_z\theta \right)-{\overline{r}}_{,\xi },$

2) $\frac{d{\overline{U}}_z}{d\xi }=\left(1+{\epsilon }_2\right)\left({\overline{d}}_z-{\overline{d}}_r\theta \right)-{\overline{z}}_{,\xi },$                               (12)

3) $\frac{d\gamma }{d\xi }=-\frac{{\overline{s}}_{,\xi }}{\overline{h}{f}_h^3}\sqrt{\frac{E_1}{E_2}}{\overline{M}}_{2\phi }-{\nu }_{12}{\overline{s}}_{,\xi }\Delta {\overline{k}}_1,$

$\Delta {\overline{k}}_1=\frac{1}{{\overline{s}}_{,\xi }\overline{r}}\left[{\overline{z}}_{,\xi }\left(\cos \gamma -1\right)+{\overline{r}}_{,\xi }\sin \gamma \right], \Delta {\overline{k}}_2=\frac{{\gamma }_{,\xi }}{{\overline{s}}_{,\xi }};$              (13)

4) $\frac{d{\overline{T}}_{2r}}{d\xi }=-A{\overline{T}}_{2r}+\frac{{\overline{s}}_{,\xi }}{\overline{r}}B{\overline{T}}_{1\phi }-C\frac{{\overline{s}}_{,\xi }\left(1+{\epsilon }_2\right)}{\overline{h}}{q}_r+D{\ddot{\overline{U}}}_r,$

5) $\frac{d{\overline{T}}_{2z}}{d\xi }=-A{\overline{T}}_{2z}-C\frac{{\overline{s}}_{,\xi }\left(1+{\epsilon }_2\right)}{\overline{h}}{q}_z+D{\ddot{\overline{U}}}_z,$                  (14)

6) $\begin{array}{l}\frac{d{\overline{M}}_{2\phi }}{d\xi }=-A{\overline{M}}_{2\phi }-\frac{{\overline{r}}_{,\xi }}{\overline{r}}B{\overline{M}}_{1s}-\\ -\frac{12}{\overline{h}}\left(1+{\epsilon }_2\right)\left({\overline{T}}_{2r}{\overline{d}}_z-{\overline{T}}_{2z}{\overline{d}}_r\right)-D\overline{h}\ddot{\gamma },\end{array}$                       (15)

дополнительные выражения и обозначения:

$$\begin{gathered} {\overline{d}}_r=\left({\overline{r}}_{,\xi }\cos \gamma -{\overline{z}}_{,\xi }\sin \gamma \right), \\ {\overline{d}}_z=\left({\overline{r}}_{,\xi }\sin \gamma +{\overline{z}}_{,\xi }\cos \gamma \right), \\ \theta =\frac{{\overline{d}}_z{\overline{T}}_{2r}-{\overline{d}}_r{\overline{T}}_{2z}}{\kappa G_{23}{\overline{s}}_{,\xi }{f}_{h}C}, A=\left(\frac{{\overline{r}}_{,\xi }}{\overline{r}}+\frac{{\epsilon }_{\mathrm{1,}\xi }}{\left(1+{\epsilon }_1\right)}\right), \\ {\epsilon }_{\mathrm{1,}\xi }=\left(\frac{{\overline{r}}_{,\xi }}{{\overline{r}}^2}{\overline{U}}_r-\frac{1}{\overline{r}}\frac{d{\overline{U}}_r}{d\xi }\right), B=\frac{\left(1+{\epsilon }_2\right)}{\left(1+{\epsilon }_1\right)}, \\ C=\frac{{\displaystyle 1-{\nu }_{12}{\nu }_{21}}}{{\displaystyle \sqrt{E_1{E}_2}}}, D=C\frac{{\displaystyle {\overline{s}}_{,\xi }\rho {\mathcal{l}}^2}}{{\displaystyle \left(1+{\epsilon }_1\right)}}, \\ {\epsilon }_1=\frac{1}{\overline{r}}{\overline{U}}_r, {\epsilon }_2=\exp \left({\widehat{\epsilon }}_2\right)-\mathrm{1,} \end{gathered}$$

$$\begin{gathered} {\widehat{\epsilon }}_2=\sqrt{\frac{E_1}{E_2}}\frac{{\overline{T}}_{2r}{\overline{d}}_r+{\overline{T}}_{2z}{\overline{d}}_z}{{\overline{s}}_{,\xi }{f}_h\left({\widehat{\epsilon }}_2\right)}-{\nu }_{12}{\widehat{\epsilon }}_1, \\ {f}_h\left({\widehat{\epsilon }}_2\right)=\exp \left[\frac{{\nu }_{12}{\nu }_{23}+{\nu }_{13}}{{\nu }_{12}{\nu }_{21}-1}{\widehat{\epsilon }}_1+\frac{{\nu }_{21}{\nu }_{13}+{\nu }_{23}}{{\nu }_{12}{\nu }_{21}-1}{\widehat{\epsilon }}_2\right], \\ {\widehat{\epsilon }}_1=\ln \left(1+{\epsilon }_1\right), \\ {\overline{T}}_{1\phi }=\sqrt{\frac{E_1}{E_2}}{f}_h\left({\widehat{\epsilon }}_1+{\nu }_{21}{\widehat{\epsilon }}_2\right), \end{gathered}$$

${\overline{M}}_{1s}=\sqrt{\frac{E_1}{E_2}}\overline{h}{f}_h^3\left(\Delta {\overline{k}}_1+{\nu }_{21}\Delta {\overline{k}}_2\right),$  (16)

где f_h = h^* / h – функция изменения толщины оболочки при деформировании, которая вычисляется итерационно.

Рассмотрим расчет малых колебаний осесимметричного бака с идеальной несжимаемой жидкостью. Для такой модели жидкости ее малые перемещения определяются потенциалом перемещений φ [4], а перемещения стенки бака определяются системой дифференциальных уравнений движения осесимметричной оболочки вращения (12)–(15), которые в сокращенном виде можно записать следующим образом:

$\frac{dX}{d\xi }=F\left(\xi ,X,\ddot{X}\right),$ (17)

где $X\left(\xi ,t\right)={\left\{\begin{array}{cccccc}{U}_r& U_z& \gamma & T_{2r}& T_{2z}& M_{2\phi }\end{array}\right\}}^T$ – вектор разрешающих функций. В случае малых колебаний решение разыскивается в следующем виде:

$X\left(\xi ,t\right)=\Delta X\left(\xi \right)\exp \left(i\omega t\right),$ (18)

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

$\frac{d\Delta X}{d\xi }=\left[\frac{\partial F}{\partial X}-{\omega }^2\frac{\partial F}{\partial \ddot{X}}\right]\Delta X.$ (19)

В рамках допущений для потенциала малых перемещений жидкости φ получим краевую задачу 4:

∆φ = 0 – в объеме жидкости, (20)

φ = 0 – на свободной поверхности жидкости, (21)

∂φ / ∂n = 0 – на оси симметрии и/или неподвижной части бака, (22)

$\overrightarrow{U}\overrightarrow{n}=\partial \phi \text{/}\partial n$, ${\overrightarrow{T}}_2\overrightarrow{n}=-{\rho }_{\ae }\phi {\omega }^2$ – на смоченной части оболочки, (23)

где $\overrightarrow{n}$ – внешняя нормаль.

Для аппроксимации производной по координатам в уравнении применяется шаблон центральной разности либо интерполяционный многочлен Лагранжа по четырем точкам для повышения точности аппроксимации при малом числе разбиений меридиана [5].

Для решения краевой задачи (20)–(23) применим метод конечных разностей, а для генерирования весовых коэффициентов разностных схем на произвольном трафарете будем использовать радиальные базисные функции (РБФ). В англоязычной литературе такой подход называется Radial Basis Function Finite Difference method (RBF-FD). В данной работе используется полигармоническая РБФ [6–8]:

ϕ(r) – r^2m ln(r), (24)

где m – степень четной полигармонической РБФ.

Использование РБФ позволяет вычислять весовые коэффициенты конечных разностей на произвольной (нерегулярной) сетке, имеющей сгущение там, где это требуется (например, в области взаимодействия жидкости и оболочки), а в других местах иметь разреженную структуру для экономии ресурсов ЭВМ.

Быстрый алгоритм генерирования нерегулярной сетки узлов в прямоугольной области предложен в [9]. Этот алгоритм был слегка модифицирован для генерирования сетки из предварительно дискретизированной выпуклой граничной области с произвольной дискретизацией. Примеры результата работы этого алгоритма приведены на рис. 3.

Рис. 3. Примеры заполнения расчетной области узлами.

На рис. 3 кроме граничных и внутренних точек области показаны законтурные (вспомогательные) узлы (в англ. ghost nodes). Они вводятся для уравнивания числа уравнений и неизвестных. Таким образом, уравнение ∆φ = 0 должно выполняться во всех узлах границы и внутренней области, а законтурные узлы участвуют в формировании трафаретов конечных разностей. На рис. 3 приведен пример трафарета, составленного из 22 точек. Такие трафареты строятся для каждого узла контура и внутри области, и по ним вычисляются весовые коэффициенты конечных разностной аппроксимации требуемых производных.

Рассмотрим алгоритм вычисления весовых коэффициентов для разностных схем. Интерполяционный сплайн на основе РБФ записывается в стандартном виде с дополнительным полиномом первой степени [10], гарантирующим уникальность решения [6], а также условиями ортогональности весовых коэффициентов следующим образом:

$S\left(x,y\right)=\sum _{k=1}^n{\lambda }_k\varphi \left(x-x_k,y-y_k\right)+{\gamma }_1+{\gamma }_{2}x+{\gamma }_{3}y,$

$\sum _{k=1}^n{\lambda }_k=\sum _{k=1}^n{\lambda }_k{x}_k=\sum _{k=1}^n{\lambda }_k{y}_k=\mathrm{0,}$ (25)

где r_k = (x_k, y_k) – узлы интерполяции, n – число узлов интерполяции, λ_k – весовые коэффициенты РБФ, γ_{1, 2, 3} – коэффициенты дополнительного многочлена. Сплайн (25) может быть расширен до пространственного случая, а также до произвольного количества размерностей.

Пусть интерполируемая функция имеет узловые значения f_k, тогда из СЛАУ (25) для определения весовых коэффициентов будет иметь вид:

$\underset{\left[A\right]}{\underbracket{\left[\begin{array}{ccccccc}\varphi \left({\overrightarrow{r}}_1-{\overrightarrow{r}}_1\right)& \varphi \left({\overrightarrow{r}}_1-{\overrightarrow{r}}_2\right)& \mathrm{...}& \varphi \left({\overrightarrow{r}}_1-{\overrightarrow{r}}_n\right)& 1& x_1& y_1\\ \varphi \left({\overrightarrow{r}}_2-{\overrightarrow{r}}_1\right)& \varphi \left({\overrightarrow{r}}_2-{\overrightarrow{r}}_2\right)& \mathrm{...}& \varphi \left({\overrightarrow{r}}_2-{\overrightarrow{r}}_n\right)& 1& x_2& y_2\\ \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}\\ \varphi \left({\overrightarrow{r}}_n-{\overrightarrow{r}}_1\right)& \varphi \left({\overrightarrow{r}}_n-{\overrightarrow{r}}_2\right)& \mathrm{...}& \varphi \left({\overrightarrow{r}}_n-{\overrightarrow{r}}_n\right)& 1& x_n& y_n\\ 1& 1& \mathrm{...}& 1& 0& 0& 0\\ x_1& x_2& \mathrm{...}& x_n& 0& 0& 0\\ y_1& y_2& \mathrm{...}& y_n& 0& 0& 0\end{array}\right]}}$ $\left\{\begin{array}{c}{\lambda }_1\\ {\lambda }_2\\ \mathrm{...}\\ {\lambda }_n\\ {\gamma }_1\\ {\gamma }_2\\ {\gamma }_3\end{array}\right\}=\left\{\begin{array}{c}{f}_1\\ f_2\\ \mathrm{...}\\ f_n\\ 0\\ 0\\ 0\end{array}\right\}.$  (26)

СЛАУ (26) имеет уникальное решение, если узлы интерполяции не повторяются (это условие обеспечивается генератором распределения узлов), а ее матрица [A] не меняет свой вид при аппроксимации произвольного дифференциального оператора L [10], что легко показать. СЛАУ для определения весовых коэффициентов w_j разложения дифференциальной оператора по узловым значениям функции f_j, в которой матрица совпадает с (26):

$\left[A\right]\cdot \left\{\begin{array}{c}{w}_1\\ w_2\\ \mathrm{...}\\ w_n\\ w_{n+1}\\ w_{n+2}\\ w_{n+3}\end{array}\right\}={\left.\left\{\begin{array}{c}L\varphi \left(x-x_1,y-y_1\right)\\ L\varphi \left(x-x_2,y-y_2\right)\\ \mathrm{...}\\ L\varphi \left(x-x_n,y-y_n\right)\\ L1\\ Lx\\ Ly\end{array}\right\}\right|}_{\begin{array}{l}x\text{ = }{x}_c\\ y\text{ = }{y}_c\end{array}}.$ (27)

Расчетная область, занятая жидкостью, разбивалась неравномерной сеткой, имеющей сгущение на границах и разреженность внутри (рис. 3). Это дает существенную экономию размерности задачи и незначительно сказывается на точности решения. Алгоритм решения уравнений ортотропной осесимметричной оболочки совместно с жидкостью добавлен в функциональные возможности программы DARSYS [5].

Рассмотрим несколько примеров расчета статического деформирования. Цилиндрической оболочки под действием внутреннего распределенного давления P (рис. 4). Радиус R = 0.5 м, толщина h = 0.02 м, высота Н = 0.7 м, модуль упругости материала E = 2 ∙ 10⁷ Па, коэффициент Пуассона v = 0.3.

Рис. 4. Цилиндрическая оболочка. Результаты расчета.

На рис. 4 приведены зависимости полного перемещения верхней точки меридиана от числа разбиений, а также деформированные конфигурации, рассчитанные в ANSYS разными конечными элементами и по формулам (12)–(16) в DARSYS при значении давления P = 1.1 ∙ 10⁵ Па при свободном опирании и с защемленным основанием на рис. 5. Кроме того, на рис. 5 показаны деформированные конфигурации меридиана с h = 0.2 м, чтобы показать влияние учета поперечного сдвига.

Рис. 5. Защемленная цилиндрическая оболочка. Результаты расчета.

Рассмотрим задачу о деформировании эллиптической оболочки с полуосями α = 0.4 м, b = 2a под действием внутреннего давления P (рис. 6). Толщина и материал взяты из предыдущего примера. На рис. 6 приведены зависимости полного перемещения точки меридиана, лежащей на малой полуоси ${\overrightarrow{i}}_r$ (ξ = 0) от числа разбиений, рассчитанные в ANSYS и по формулам (12)–(16) при значении давления P = 3.3 ∙ 10⁵ Па и деформированная конфигурация. Из рисунков видно, что результаты хорошо согласуются.

Рис. 6. Эллиптическая оболочка. Результаты расчета.

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

Рис. 7. Цилиндрическая (а) и эллиптическая (б) оболочки, зависимость перемещения от давления.

Рассмотрим задачу о деформировании конической оболочки под действием внутреннего распределенного давления P, радиус окружности дна R = 0.5 м, угол конуса φ = π / 6 (рис. 8), толщину и параметры материала взяты из предыдущих примеров. На рис. 8 приведена деформированная конфигурация и зависимости полного перемещения вершины конуса (ξ = 0) от числа разбиений, рассчитанные в ANSYS и по формулам (12)–(16) при значении давления P = 3.5 ∙ 10⁵ Па.

Рис. 8. Коническая оболочка. Результаты расчета.

Рассмотрим задачу о деформировании составной оболочки (рис. 9) под действием внутреннего распределенного давления P, радиус цилиндра R = 0.5 м, высота цилиндра H = 0.8 м, длина образующей конуса L = 0.3 м, φ = π / 6, толщина и материал взяты из предыдущих примеров. На рис. 9 приведена деформированная конфигурация и зависимости полного перемещения верхней точки меридиана составной оболочки (ξ = 0) от числа разбиений, рассчитанные в ANSYS и по формулам (12)–(16) при значении давления P = 1.8 ∙ 10⁵ Па. На рис. 10 представлены зависимости перемещения от давления для цилиндрической (а) и эллиптической (б) оболочки.

Рис. 9. Составная оболочка. Результаты расчета.

Рис. 10. Коническая (а) и составная (б) оболочки, зависимость перемещения от давления.

Рассмотрим тестовую задачу из статьи [11]: полусферический открытый бак заполнен водой, шарнирно закреплен по радиусу. Радиус R = 5.08 м, толщина 0.0254 м, модуль упругости Е = 7 ∙ 10¹⁰ Па, коэффициент Пуассона v = 0.3, плотность материала ρ = 2770 кг/м³.

В табл. 1 приведены три низшие частоты осесимметричных гидроупругих колебаний, рассчитанные различными методами в разных программах. Три первые столбца с частотами взяты из статьи [12], в которой нет информации о сходимости в зависимости от дискретизации модели. В табл. 1 в четвертом столбце приведены частоты, которые практически не меняются при увеличении числа конечных и граничных элементов, при расчете в программе В.Е. Левина [4] по МКЭ-МГЭ. В пятом столбце приведены частоты, рассчитанные в ANSYS Workbench 14.5 при ~50 разбиениях вдоль меридиана.

Таблица 1. Частоты колебаний полусферической оболочки с водой, Гц

Из табл. 1 видно, что расхождение в первой частоте между разными методами достигает 16%. Для анализа сходимости были проведены расчеты с разной дискретизацией модели. На рис. 11 приведены графики сходимости рассчитанных частот в ANSYS и по разработанной методике. Из рисунка видно, что графики пересекаются с ростом числа разбиений. Сходимости частот в ANSYS Workbench 14.5 (Shell181+Fluid80) получить не удалось, а частоты, рассчитанные по предлагаемой методике, имеют выраженную асимптотику при увеличении числа разбиений меридиана. На рис. 11б приведены зависимости низшей частоты от дискретизации меридиана для различных размеров трафарета конечных разностей n = 15, 22, 50 для степени полигармонической РБФ m = 2. Повышение степени m не приводит к существенным изменениям, как и увеличение размера трафарета n.

Рис. 11. Сходимость частот гидроупругих колебаний полусферической оболочки с водой.

В табл. 2 приведены частоты собственных колебаний составной оболочки из изотропного материала (рис. 9), заполненной водой до уровня 0.5 м. Были проведены расчеты частот собственных осесимметричных колебаний в программе [4] по МКЭ-МГЭ (~60 элементов), в ANSYS Workbench 14.5 моделировалась четверть модели с условиями симметрии (~119 тыс. элементов) и по предложенной методике в программе DARSYS (300 разбиений меридиана), результаты сведены в табл. 2.

Таблица 2. Частоты колебаний составной оболочки с водой, Гц

Из табл. 2 видно, что наибольшее различие приходится на низшую частоту и достигает 11%, с ростом номера тона разница падает до 0.5%. Уточнить результаты в ANSYS путем увеличения числа КЭ не удалось.

Рассмотрим тороидальный бак, наполовину заполненный водой, образующая окружность имеет радиус 0.5 м, внутренний радиус – 0.5 м, внешний – 1.5 м. Толщина оболочки 0.02 м, изотропный материал: модуль упругости E = 2 ∙ 10⁷ Па, коэффициент Пуассона v = 0.3, плотность ρ = 7850 кг/м³. Тор шарнирно закреплен по внешнему радиусу. На рис. 12 приведены частоты и формы колебаний, рассчитанные в программе [4] (140 КЭ, 140 ГЭ вдоль меридиана) и по предлагаемой методике в программе DARSYS [5] (140 разбиений меридиана, m = 2, n = 22). В табл. 3 приведено сравнение первых 6 частот осесимметричных колебаний. Из табл. 3 видно, что частоты хорошо согласуются между собой.

Рис. 12. Формы колебаний тороидального бака (МКЭ-МГЭ [4] (а), DARSYS (б), сходимость частот (в)).

Таблица 3. Частоты колебаний тороидального бака с водой, Гц

На рис. 12 приведены зависимости 3 низших частот от числа разбиений меридиана, рассчитанные в DARSYS. Из рис. 12 видно, что частоты стремятся к горизонтальным прямым с ростом числа разбиений без заметных осцилляций.

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

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

С. Nguyen

Хат алмасуға жауапты Автор.
Email: mckq1985@gmail.com
Ресей, Novosibirsk

D. Shelevaya

Email: mckq1985@gmail.com
Ресей, Novosibirsk; Novosibirsk

D. Krasnorutsky

Email: mckq1985@gmail.com
Ресей, Novosibirsk; Novosibirsk

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

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

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

Әрекет

1. JATS XML

Жүктеу

2. Fig. 1. Deformation and equilibrium of a small shell element.

Жүктеу (137KB)

Метадеректер

3. Fig. 2. Axisymmetric shell of revolution with liquid.

Жүктеу (195KB)

Метадеректер

4. Fig. 3. Examples of filling the calculation area with nodes.

Жүктеу (1MB)

Метадеректер

5. Fig. 4. Cylindrical shell. Calculation results.

Жүктеу (368KB)

Метадеректер

6. Fig. 5. Clamped cylindrical shell. Calculation results.

Жүктеу (522KB)

Метадеректер

7. Fig. 6. Elliptical shell. Calculation results.

Жүктеу (350KB)

Метадеректер

8. Fig. 7. Cylindrical (a) and elliptical (b) shells, dependence of displacement on pressure.

Жүктеу (297KB)

Метадеректер

9. Fig. 8. Conical shell. Calculation results.

Жүктеу (304KB)

Метадеректер

10. Fig. 9. Composite shell. Calculation results.

Жүктеу (391KB)

Метадеректер

11. Fig. 10. Conical (a) and composite (b) shells, dependence of displacement on pressure.

Жүктеу (275KB)

Метадеректер

12. Fig. 11. Convergence of frequencies of hydroelastic oscillations of a hemispherical shell with water.

Жүктеу (283KB)

Метадеректер

13. Fig. 12. Forms of oscillations of a toroidal tank (FEM-MGE [4] (a), DARSYS (b), frequency convergence (c)).

Жүктеу (281KB)

Метадеректер

Presented by Academician of the RAS B.D. Annin