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

Abstract

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

Full Text

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

Важной задачей при анализе закритического поведения композитных пластин является обоснование эффективной и адекватной кинематической модели для точного моделирования глобального и локального поведения многослойных структур при действии плоскостных нагрузок [11, 12]. Классическая теория пластин Кирхгофа, успешно применяемая для тонких пластин, не учитывает влияние деформации поперечного сдвига, которые значимы для пластин с низкой поперечной жесткостью на сдвиг. В результате анализ на основе этой простой теории завышает параметры критического состояния пластины и характеристики динамического отклика. В теории деформации сдвига первого порядка Миндлина поперечные сдвиги учтены в предположении, что нормаль к срединной поверхности пластины остается прямой при изгибе, но не обязательно перпендикулярной ей [13]. Эта теория широко используются для инженерного анализа тонких и умеренно толстых композитных структур из-за ее простоты и низких вычислительных затрат. Однако эта теория содержит дополнительный параметр – поправочный коэффициент для сведения к нулю напряжений сдвига на верхней и нижней поверхностях панели. От этого недостатка свободны теории деформации сдвига более высокого порядка, в которых предполагают нелинейное изменение поперечных касательных напряжений по толщине пластины. Эти модели не требуют введения поправочного коэффициента на сдвиг и обладают достаточной точностью. Теории деформации сдвига более высокого порядка в полиномиальной форме основаны на разложении перемещений в плоскости как квадратичных или кубических функций координаты толщины, что обеспечивает удовлетворение требований отсутствия поперечных напряжений сдвига на верхней и нижней поверхностях ламината. Нгуен [14] разработал новую унифицированную обобщенную форму теории деформаций сдвига высокого порядка для слоистых композитных пластин, в которой все существующие модели сдвиговой деформации математически представлены уникальной полиномиальной формулировкой, которой вполне достаточно для описания любого нелинейного распределения поперечных касательных напряжений. Известны также некоторые неполиномиальные теории деформации сдвига высокого порядка с тригонометрическими [2, 3], экспоненциальными [4] и гиперболическими [9] функциями, которые достаточно точно аппроксимируют поле касательных напряжений и граничные условия, но значительно усложняют вычислительные алгоритмы. Существуют также зигзагообразные теории [4], которые основаны на полных квадратичных или кубических полях с кусочно-линейным зигзагообразным представлением перемещений в плоскости и автоматически удовлетворяющие условия непрерывности межслойных и поперечных касательных напряжений и условиям на поверхности.

Одним из преимуществ аналитических методов является закрытая форма представления расчетных моделей для различных режимов потери устойчивости, и, следовательно, возможности последующей параметрической оптимизации с помощью известных алгоритмов. Однако точные аналитические решения доступны только для некоторых простых задач при определенных нагрузках и граничных условиях [15]. Поэтому получили развитие эффективные процедуры приближенных и численных решений задач закритического деформирования композитных панелей. К таким методам относятся метод Галеркина [16], метод Релея–Ритца [6–8], дифференциально-квадратурный метод [10], бессеточные методы [9], метод дискретной сингулярной свертки, расширенный метод Канторовича, метод конечных полос, метод конечных элементов [3] и др. Обзор литературы показывает, что метод Рэлея–Ритца является одним из наиболее успешных подходов к исследованиям поведения композитных пластин после потери устойчивости, обеспечивающим достаточную точность и высокую вычислительную эффективность. Важнейшим этапом этого энергетического метода является выбор подходящей функции перемещения, чтобы правильно описать нелинейное поведение пластины и удовлетворить граничным условиям. Критическая нагрузка потери устойчивости определяется путем решения характеристического уравнения в соответствии с принципом минимальной потенциальной энергии. Метод Рэлея–Ритца часто используется для анализа больших прогибов композитных панелей и тонкостенных конструкций на основе теории сдвига первого порядка [4], однако недостатком метода Ритца является его негибкость по отношению к граничным условиям [15]. Решением этой проблемы может быть применение нелинейной модели В. Оливери [17], в которой для аппроксимации перемещений предлагается схема интерполяции pb‑2 и стандартные ортогональные полиномы Лагранжа вместе с так называемой граничной функцией, способной гарантировать выполнение геометрических граничных условий.

Большинство исследовательских работ посвящено анализу закритического состояния прямоугольных или косоугольных композитных панелей, и лишь ограниченные решения применены для трапециевидных или общих четырехугольных анизотропных пластин. В целом анализ косоугольных и трапециевидных пластин сложнее, чем анализ прямоугольных [18]. Это связано с особенностями нелинейного деформирования непрямоугольных пластин с произвольными граничными условиями – влиянием угла скоса сторон пластины [19]; неравномерностью распределение напряжения до и после потери устойчивости [20]; сингулярностью напряжений в тупых углах панели [21] и др. Применение метода конечных элементов и других сетчатых методов ограничено недостаточной точностью при искажении сетки из-за нелинейной связи между изопараметрическими координатами (x, h) и декартовыми координатами (х, у) в случае непрямоугольной расчетной области [22].

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

Основные соотношения. Рассмотрим четырехугольную k-слойную композитную пластину из ортотропного материала со слоями постоянной толщины c различной структурой, определяемой углом армирования ±φ. Геометрия панели (рис. 1) характеризуется размерами продольных λ₁ ≠ λ₂ и поперечных сторон b₁ ≠ b₂, толщиной h, углами скоса χ₁, χ₂ и углом q между левой и правой кромками пластины.

Рис. 1. Геометрические параметры и системы координат панели и элементарного композитного слоя.

Введем косоугольную систему координат x, h, направления осей (контравариантный базис) заданы таким образом, что ось h совпадает с направлением оси y исходной системы декартовых координат, а направление оси ξ составляет угол χ_₁ с направлением оси x.

Это позволяет описать метрики неортогональных четырехугольных объектов в косоугольном базисе, в котором можно преобразовать тензоры деформаций и напряжений, связанные с локальной системой армирования элементарных ортотропных слоев композитного материала. В косоугольном базисе ξ, h компоненты ковариантных и контравариантных тензоров деформаций и напряжений определяются известными соотношениями Лурье [23]. Напомним, что контравариантная матрица, преобразующая компоненты вектора, должна быть обратной по отношению к матрице, преобразующей базисные векторы, а ковариантный вектор преобразуется той же матрицей, что и изменение базовой матрицы. Тогда основные соотношения механики композиционных материалов в косоугольном базисе получим в виде:

• – условия равновесия

$\begin{array}{c}{N}_{\xi ,\xi }+N_{\xi \eta ,\eta }=0,N_{\xi \eta ,\xi }+N_{\eta ,\eta }=\mathrm{0,}\\ M_{\xi ,\xi \xi }+2M_{\xi \eta ,\xi \eta }+M_{\eta ,\eta \eta }+N_{\xi }{w}_{,\xi \xi }+2N_{\xi \eta }{w}_{,\xi \eta }+N_{\eta }{w}_{,\eta \eta }=0.\end{array}$                                                        (1)

• – нелинейные геометрические соотношения

${\epsilon }_{\xi }=u_{\xi ,\xi }+\frac{{(w_{,\xi })}^2}{2};{\epsilon }_{\eta }=u_{\eta ,\eta }+\frac{{(w_{,\eta })}^2}{2};{\epsilon }_{\xi \eta }=u_{\xi ,\eta }+u_{\eta ,\xi }+w_{,\xi }{w}_{,\eta };$

• – физические соотношения

$N_{\xi \eta }=N_{\xi }{N}_{\eta }{N}_{\xi \eta }={{\displaystyle \int }}_{-h/2}^{h/2}{\sigma }_{\xi }{\sigma }_{\eta }{\tau }_{\xi \eta }dz={\left[B\right]}^{\xi \eta }{\epsilon }_{\xi \eta }+{\left[C\right]}^{\xi \eta }{\zeta }_{\xi \eta },$

$M_{\xi \eta }={{\displaystyle \int }}_{-h/2}^{h/2}{\sigma }^{\left(p\right)}{}_{\xi \eta }zdz={\left[C\right]}^{\xi \eta }{\epsilon }_{\xi \eta }+{\left[D\right]}^{\xi \eta }{\zeta }_{\xi \eta },$

где ${[B]}^{\xi \eta },\text{ [}C{\text{]}}^{\xi \eta },\text{ [}D{\text{]}}^{\xi \eta }$ – матрицы мембранных, смешанных и изгибных жесткостей соответственно.

Зависимости, связывающие деформации и напряжения в косоугольном базисе с деформациями и напряжениями элементарного слоя имеют вид [6]

$\begin{array}{c}{\left[A\right]}^{\xi \eta }={\left[R_{\sigma }^{\xi \eta }\right]}^{-1}{\left[A\right]}_{\alpha \beta }=\left[R_{\epsilon }^{\xi \eta }\right]=\\ =\left[\begin{array}{c}\sec \chi 00\\ \sin \chi \text{tg}\chi \cos \chi -2\sin \chi \\ -\text{tg}\chi 01\end{array}\right]\left[\begin{array}{c}{\sec }^2\chi {\text{\hspace{1em}tg}}^2\chi -\text{tg}\chi \sec \chi \\ 010\\ 0-2\text{tg}\chi \sec \chi \end{array}\right].\end{array}$

Здесь ${[A]}^{\xi \eta }$ и ${[A]}_{\alpha \beta }$ – матрицы обобщенных жесткостей в косоугольном базисе и в локальной системе, связанной с направлением армирования [2].

$\begin{array}{c}{e}_{\alpha \beta }={\epsilon }_{\alpha \beta }+z{\zeta }_{\alpha \beta }=[R_{\alpha \beta }^{\epsilon }][R_{\epsilon }^{\xi \eta }]e_{\xi \eta };{\sigma }_{\alpha \beta }=[R_{\alpha \beta }^{\sigma }][R_{\sigma }^{\xi \eta }]{\sigma }_{\xi \eta };\\ R_{\alpha \beta }^{\sigma }=\left[\begin{array}{c}{\cos }^2\left(\gamma \pm \phi \right) {\sin }^2\left(\gamma \pm \phi \right) \sin 2\left(\gamma \pm \phi \right)\\ {\sin }^2\left(\gamma \pm \phi \right) {\cos }^2\left(\gamma \pm \phi \right) -\sin 2\left(\gamma \pm \phi \right)\\ -0.5\sin 2\left(\gamma \pm \phi \right) 0.5\sin 2\left(\gamma \pm \phi \right) \cos 2\left(\gamma \pm \phi \right)\end{array}\right];\\ [R_{\sigma }^{\xi \eta }]=\left[\begin{array}{c}\cos \chi 0 0\\ \sin \chi \text{tg}\chi \sec \chi 2\text{tg}\chi \\ \sin \chi 0 1\end{array}\right].\end{array}$

Анализ потери устойчивости композитной панели. Полагаем, что форма закритического прогиба будет подобна первой собственной форме, которая соответствует минимальному собственному значению. Решение задачи на собственные значения четырехугольной панели выполним с помощью метода Релея–Ритца, для чего функцию перемещений представим в виде двойного ряда, образованного функциями Крылова [6, 9]

$w={\displaystyle \sum }_{m=1}^{\infty }{\displaystyle \sum }_{n=1}^{\infty }{A}_{mn}{\overline{w}}_{1m}\left(\xi \right){\overline{w}}_{2n}\left(\eta \right).$                                                                             (2)

Для определения минимальных собственных значений и принадлежащих им собственных векторов и форм, полную потенциальную энергию E панели площадью S представим в форме Брайана [3]

$\begin{array}{c}E=\frac{1}{2}{\displaystyle {\iint }_S\left\{D_{11}^{\xi \eta }{\left(w_{,\xi \xi }\right)}^2+2D_{12}^{\xi \eta }\left(w_{,\xi \xi }\right)\left(w_{,\eta \eta }\right)+D_{22}^{\xi \eta }{\left(w_{,\eta \eta }\right)}^2+\right.4D_{33}^{\xi \eta }{\left(w_{,\xi \eta }\right)}^2-}\\ \left.-N_{\xi }{\left(w_{,\xi }\right)}^2-N_{\eta }{\left(w_{,\eta }\right)}^2-2N_{\xi \eta }\left(w_{,\xi }\right)\left(w_{,\eta }\right)\right\}d\xi d\eta .\end{array}$             (3)

Подставляя в функционал энергии E прогиб (2), вычислим интеграл (3) методом Симпсона [24], согласно которому область интегрирования разбивается на i× j $(i\in \{\mathrm{0,}\mathrm{...,}{i}_{\max }\},j\in \{\mathrm{0,}\mathrm{...,}{j}_{\max }\})$ ячеек сетки. Минимизируя получающееся после численного интегрирования выражение по каждому элементу неизвестного вектора амплитуд $⌊A⌋,$ получим однородную систему линейных алгебраических уравнений

$(\left[\Omega \right]-{\lambda }_{\xi }[t_{\xi }]-{\lambda }_{\eta }[t_{\eta }]-{\lambda }_{\xi \eta }[t_{\xi \eta }])⌊A⌋=\mathrm{0,}$                                                                     (4)

где $[\Omega ],[t_{\xi }],[t_{\eta }],[t_{\xi \eta }]$ – квадратные матрицы [6]; ${\lambda }_{\xi },{\lambda }_{\eta },{\lambda }_{\xi \eta }$ – векторы собственных значений.

Условие равенства нулю определителя однородной системы уравнений (4)

$\det (\left[\Omega \right]-{\lambda }_{\xi }[t_{\xi }]-{\lambda }_{\eta }[t_{\eta }]-{\lambda }_{\xi \eta }[t_{\xi \eta }])=\mathrm{0,}$

позволяет определить векторы собственных значений $\left|{\lambda }_{\xi }\right|$, $\left|{\lambda }_{\eta }\right|$ и $\left|{\lambda }_{\xi \eta }\right|$  с помощью двойного QR/QL-алгоритма Френсиса [24], а также несколько первых собственных векторов $\left|q_1\right|, \left|q_2\right|, \left|q_3\right|$ …, соответствующих собственным значениям ${\lambda }_1< {\lambda }_2< {\lambda }_3$ …. Затем вычисляем необходимые для аппроксимации закритического прогиба собственные формы в виде ${\overline{w}}_i={\left\{{{\displaystyle \sum }}_p{q}_p{\tilde{w}}_1(\xi ){\tilde{w}}_2(\eta )\right\}}_i,$ где р – число членов ряда (3).

Закритическое деформирование ортотропных четырехугольных панелей с произвольными граничными условиями. Перейдем к закритической стадии деформирования и представим закритический прогиб панели в виде произведения неизвестной амплитуды A₀ и собственной формы

$w=A_0{{\displaystyle \sum }}_p{q}_p{\tilde{w}}_1\left(\xi \right){\tilde{w}}_2\left(\eta \right).$                                                                                    (5)

Здесь ${\tilde{w}}_1\left(\xi \right),{\tilde{w}}_2\left(\eta \right)$ – функции, описывающие собственные формы в продольном и поперечном направлениях.

Перемещения в каждом узле сетки Симпсона представим в виде полиномов второй степени

$\begin{array}{c}{u}_{\xi }=C_1^{\left(i,j\right)}\eta +C_2^{\left(i,j\right)}\xi +C_3^{\left(i,j\right)}\frac{{\xi }^2}{2l_{sr}} ; u_{\eta }=C_4^{\left(i,j\right)}\xi +C_5^{\left(i,j\right)}\eta +C_6^{\left(i,j\right)}\frac{{\eta }^2}{2}{b}_{sr},\\ l_{sr}=\raisebox{1ex}{$(l_1+l_2)$}\!\left/ \!\raisebox{-1ex}{$2$}\right.; b_{sr}=\raisebox{1ex}{$(b_1+b_2)$}\!\left/ \!\raisebox{-1ex}{$2$}\right.,\end{array}$

где $C_1^{(i,j)},C_2^{(i,j)},C_3^{(i,j)},C_4^{(i,j)},C_5^{(i,j)},C_6^{(i,j)}$ – относительные параметры прогиба панели в узлах сетки $q_p{\tilde{w}}_1(i){\tilde{w}}_2(j),$ соответствующие граничным условиям

$\begin{array}{c}\xi =0\text{:} N_{\eta }{l}_1=T_{\eta }{l}_1; N_{\xi \eta }{l}_1=T_{\xi \eta }{l}_1; \eta =0:N_{\xi }{b}_1=T_{\xi }{b}_1; N_{\xi \eta }{b}_1=T_{\xi \eta }{b}_1;\\ \xi =l_2\text{:} N_{\xi }{b}_2=T_{\xi }{b}_2-T_{\xi \eta }{b}_{2}tg\theta ; N_{\xi \eta }{b}_2=T_{\xi \eta }{b}_2-T_{\xi }{b}_{2}tg\theta ;\\ \eta =b_2\text{:} N_{\eta }{l}_2=T_{\eta }{l}_2-T_{\xi \eta }{l}_{2}tg{\chi }_{sr}; N_{\xi \eta }{l}_2=T_{\xi \eta }{l}_2-T_{\eta }{l}_{2}tg{\chi }_{sr}; {\chi }_{sr}=\raisebox{1ex}{$\left({\chi }_1+{\chi }_2\right)$}\!\left/ \!\raisebox{-1ex}{$2$}\right..\end{array}$

Подставив перемещения и собственную форму (5) в уравнения равновесия, определим оставшиеся $C_3^{(i,j)},C_6^{(i,j)}$

$\begin{array}{c}{C}_3^{\left(i,j\right)}=-\frac{F_1^{\left(i,j\right)}}{B_{11}^{\xi \eta }};C_6^{\left(i,j\right)}=-\frac{F_2^{\left(i,j\right)}}{B_{22}^{\xi \eta }},\\ F_1^{\left(i,j\right)}=B_{11}^{\xi \eta }{w}_{,\xi }{w}_{,\xi \xi }+B_{33}^{\xi \eta }{w}_{,\xi }{w}_{,\eta \eta }+(B_{12}^{\xi \eta }+B_{33}^{\xi \eta })w_{,\eta }{w}_{,\xi \eta }=\raisebox{1ex}{$B_{11}^{\xi \eta }{A}_0^2{f}_1(i,j)$}\!\left/ \!\raisebox{-1ex}{$2$}\right.,\\ F_2^{\left(i,j\right)}=(B_{21}^{\xi \eta }+B_{33}^{\xi \eta })w_{,\xi }{w}_{,\xi \eta }+B_{22}^{\xi \eta }{w}_{,\eta }{w}_{,\eta \eta }+B_{33}^{\xi \eta }{w}_{,\eta }{w}_{,\xi \xi }=\raisebox{1ex}{$B_{22}^{\xi \eta }{A}_0^2{f}_2(i,j)$}\!\left/ \!\raisebox{-1ex}{$2$}\right.,\end{array}$ 

$\begin{array}{c}{f}_1(i,j)={\sum }_p{q}_p{\tilde{w}}_{\mathrm{1,}\xi }{\tilde{w}}_2{\sum }_r{q}_r{\tilde{w}}_{\mathrm{1,}\xi \xi }{\tilde{w}}_2+{\sum }_p{q}_p{\tilde{w}}_{\mathrm{1,}\xi }{\tilde{w}}_2{\sum }_r{q}_r{\tilde{w}}_1{\tilde{w}}_{\mathrm{2,}\eta \eta }\cdot B_{11}^{\xi \eta }/B_{11}^{\xi \eta }+\\ +{\sum }_p{q}_p{\tilde{w}}_1{\tilde{w}}_{\mathrm{2,}\eta }{\sum }_r{q}_r{\tilde{w}}_{\mathrm{1,}\xi }{\tilde{w}}_{\mathrm{2,}\eta }\cdot (B_{12}^{\xi \eta }+B_{33}^{\xi \eta })/B_{11}^{\xi \eta },\\ f_2\left(i,j\right)={\sum }_p{q}_p{\tilde{w}}_1{\tilde{w}}_{\mathrm{2,}\eta }{\sum }_r{q}_r{\tilde{w}}_1{\tilde{w}}_{\mathrm{2,}\eta \eta }+{\sum }_p{q}_p{\tilde{w}}_1{\tilde{w}}_{\mathrm{2,}\eta }{\sum }_r{q}_r{\tilde{w}}_{\mathrm{1,}\xi \xi }{\tilde{w}}_2\cdot B_{33}^{\xi \eta }/B_{22}^{\xi \eta }+\\ +{\sum }_p{q}_p{\tilde{w}}_{\mathrm{1,}\xi }{\tilde{w}}_2{\sum }_r{q}_r{\tilde{w}}_{\mathrm{1,}\xi }{\tilde{w}}_{\mathrm{2,}\eta }\cdot (B_{12}^{\xi \eta }+B_{33}^{\xi \eta })/B_{22}^{\xi \eta }.\end{array}$ 

Определим деформации панели

Полные деформации панели равны

$\begin{array}{c}\begin{array}{c}{e}_{\xi }={\epsilon }_{\xi }+z{\zeta }_{,\xi \xi }={\epsilon }_{\xi }-z{A}_0{\sum }_p{q}_p{\tilde{w}}_{\mathrm{1,}\xi \xi }{\tilde{w}}_2; \\ e_{\eta }={\epsilon }_{\eta }+z{\zeta }_{,\eta \eta }={\epsilon }_{\eta }-z{A}_0{\sum }_p{q}_p{\tilde{w}}_1{\tilde{w}}_{\mathrm{2,}\eta \eta };\\ e_{\xi \eta }={\epsilon }_{\xi \eta }+2z{\zeta }_{,\xi \eta }={\epsilon }_{\xi \eta }-2z{A}_0{\sum }_p{q}_p{\tilde{w}}_{\mathrm{1,}\xi }{\tilde{w}}_{\mathrm{2,}\eta }.\end{array}\end{array}$ 

Полная потенциальная энергия U закритического деформирования панели c учетом найденных деформаций примет вид

$\begin{array}{c}U=\frac{1}{2}{\displaystyle \underset{S}{\iint }\{A_0^2{\zeta }^T{[D]}^{\xi \eta }\zeta +A_0^4{\tilde{\epsilon }}^T{[B]}^{\xi \eta }\tilde{\epsilon }+A_0^2({\tilde{N}}_{\xi }{\tilde{\epsilon }}_{\xi }+{\tilde{N}}_{\eta }{\tilde{\epsilon }}_{\eta }+{\tilde{N}}_{\xi \eta }{\tilde{\epsilon }}_{\xi \eta })\}d\xi d\eta -}\\ -{\displaystyle \underset{S}{\iint }[T_{\xi }{u}_{\xi ,\xi }+T_{\eta }{u}_{\eta ,\eta }+2T_{\xi \eta }(u_{\xi ,\eta }+u_{\eta ,\xi })]d\xi d\eta ;}\end{array}$                             (6)

$\begin{array}{c}{\tilde{N}}_{\xi }=T_{\xi }-T_{\xi \eta }\frac{i}{i_{\max }}tg\theta ; {\tilde{N}}_{\eta }=T_{\eta }-T_{\xi \eta }\frac{j}{j_{\max }}tg{\chi }_{sr};\\ {\tilde{N}}_{\xi \eta }=T_{\xi \eta }-T_{\xi }\frac{i}{i_{\max }}tg\theta -T_{\eta }\frac{j}{j_{\max }}tg{\chi }_{sr}.\end{array}$

Исходя из минимума полной потенциальной энергии U в виде (6), получим уравнение третьей степени относительно амплитуды прогиба A₀

$\frac{\partial U}{\partial A_0}=-A_0({\upsilon }_{\xi }{T}_{\xi }+{\upsilon }_{\eta }{T}_{\eta }+{\upsilon }_{\xi \eta }{T}_{\xi \eta }-N^{*})+A_0^3{R}^{*}=\mathrm{0,}$                                                           (7)

$\begin{array}{c}{\upsilon }_{\xi }={\displaystyle {\iint }_S\left[{\left({{\displaystyle \sum }}_p{q}_p{\tilde{w}}_{\mathrm{1,}\xi }{\tilde{w}}_2\right)}^2-2\frac{i}{i_{\max }}tg\theta {{\displaystyle \sum }}_p{q}_p{\tilde{w}}_{\mathrm{1,}\xi }{\tilde{w}}_2{{\displaystyle \sum }}_r{q}_r{\tilde{w}}_1{\tilde{w}}_{\mathrm{2,}\eta }\right]d\xi d\eta },\\ {\upsilon }_{\eta }={\displaystyle {\iint }_S\left[{\left({{\displaystyle \sum }}_p{q}_p{\tilde{w}}_1{\tilde{w}}_{\mathrm{2,}\eta }\right)}^2-2\frac{j}{j_{\max }{}_{}}tg{\chi }_{sr}{{\displaystyle \sum }}_p{q}_p{\tilde{w}}_{\mathrm{1,}\xi }{\tilde{w}}_2{{\displaystyle \sum }}_r{q}_r{\tilde{w}}_1{\tilde{w}}_{\mathrm{2,}\eta }\right]d\xi d\eta },\\ {\upsilon }_{\xi \eta }={\displaystyle {\iint }_S\left\{2{{\displaystyle \sum }}_p{q}_p{\tilde{w}}_{\mathrm{1,}\xi }{\tilde{w}}_2{{\displaystyle \sum }}_r{q}_r{\tilde{w}}_1{\tilde{w}}_{\mathrm{2,}\eta }-\frac{i}{i_{\max }}tg\theta {{\displaystyle \sum }}_p{q}_p{\tilde{w}}_{\mathrm{1,}\xi }{\tilde{w}}_2{{\displaystyle \sum }}_r{q}_r{\tilde{w}}_{\mathrm{1,}\xi }{\tilde{w}}_2-\right.}\\ \left.-\frac{j}{j_{\max }{}_{}}tg{\chi }_{sr}{{\displaystyle \sum }}_p{q}_p{\tilde{w}}_1{\tilde{w}}_{\mathrm{2,}\eta }{{\displaystyle \sum }}_r{q}_r{\tilde{w}}_1{\tilde{w}}_{\mathrm{2,}\eta }\right\}d\xi d\eta ,\end{array}$

$\begin{array}{c}{N}^{*}={\displaystyle {\iint }_S\left\{D_{11}^{\xi \eta }{{\displaystyle \sum }}_p{q}_p{\tilde{w}}_{\mathrm{1,}\xi \xi }{\tilde{w}}_2{{\displaystyle \sum }}_r{q}_r{\tilde{w}}_{\mathrm{1,}\xi \xi }{\tilde{w}}_2+2D_{12}^{\xi \eta }{{\displaystyle \sum }}_p{q}_p{\tilde{w}}_{\mathrm{1,}\xi \xi }{\tilde{w}}_2{{\displaystyle \sum }}_r{q}_r{\tilde{w}}_1{\tilde{w}}_{\mathrm{2,}\eta \eta }+\right.}\\ \left.+D_{22}^{\xi \eta }{{\displaystyle \sum }}_p{q}_p{\tilde{w}}_1{\tilde{w}}_{\mathrm{2,}\eta \eta }{{\displaystyle \sum }}_r{q}_r{\tilde{w}}_1{\tilde{w}}_{\mathrm{2,}\eta \eta }+4D_{33}^{\xi \eta }{{\displaystyle \sum }}_p{q}_p{\tilde{w}}_{\mathrm{1,}\xi }{\tilde{w}}_2{{\displaystyle \sum }}_r{q}_r{\tilde{w}}_1{\tilde{w}}_{\mathrm{2,}\eta }\right\}d\xi d\eta ,\\ R^{*}={\displaystyle {\iint }_S\left[(B_{11}^{\xi \eta }{\tilde{\epsilon }}_{\xi }+B_{12}^{\xi \eta }{\tilde{\epsilon }}_{\eta }){\tilde{\epsilon }}_{\xi }+(B_{21}^{\xi \eta }{\tilde{\epsilon }}_{\xi }+B_{22}^{\xi \eta }{\tilde{\epsilon }}_{\eta }){\tilde{\epsilon }}_{\eta }+B_{33}^{\xi \eta }{({\tilde{\epsilon }}_{\xi \eta })}^2\right]d\xi d\eta }.\end{array}$

Корни уравнения (7) равны

$A_1=\mathrm{0,}{A}_{\mathrm{2,3}}=\sqrt{\frac{({\upsilon }_{\xi }{T}_{\xi }+{\upsilon }_{\eta }{T}_{\eta }+{\upsilon }_{\xi \eta }{T}_{\xi \eta }-N^{*})}{R^{*}}}.$                                                             (8)

Корень A₁ = 0 не имеет смысла, т.к. в этом случае перемещения равны нулю и панель остается плоской. Значения двух других корней определяют величину закритического прогиба в форме, подобной первой собственной форме. Направление прогиба, соизмеримого с толщиной панели, для симметричной укладки слоев не может быть определено.

Предельное состояние четырехугольной композитной панели. Предельное состояние композитной панели в закритической области можно установить по разрушению наиболее нагруженного слоя на основании одного из критериев прочности анизотропных материалов [25]. Соотношения (5)–(8) полностью определяют компоненты тензора деформаций слоистой панели $\left\{e_{\xi }, e_{\eta }, e_{\xi \eta }\right\}$ в косоугольном базисе и позволяют рассчитать напряженное состояние ${⌊{\sigma }_1^{(i)},{\sigma }_2^{(i)},{\tau }_{12}^{(i)}⌋}^T$ каждого i-го слоя в локальной системе α, β связанной с направлением армирования слоя (индексы 1, 2 соответствуют осям α, β). Текущему напряженному состоянию соответствует значение коэффициента нагруженности $F_p^{(i)}(T_{\xi },T_{\eta },T_{\xi \eta })\le 1$ [25], рассчитанное на основании сопоставления компонент тензора напряжений с предельными характеристиками прочности композита для различных критериев разрушения:

• – критерий Цая-Ву $F_p^{\left(i\right)}={\left(\frac{{\sigma }_1^{\left(i\right)}}{{\sigma }_{1\text{в}}}\right)}^2-\left(\frac{{\sigma }_1^{\left(i\right)}}{{\sigma }_{1\text{в}}}\frac{{\sigma }_2^{\left(i\right)}}{{\sigma }_{2\text{в}}}\right)+{\left(\frac{{\sigma }_2^{\left(i\right)}}{{\sigma }_{2\text{в}}}\right)}^2+{\left(\frac{{\tau }_{12}^{\left(i\right)}}{{\tau }_{12\text{в}}}\right)}^2;$
• – критерий максимальных напряжений $F_p^{\left(i\right)}=\max \left\{\frac{{\sigma }_1^{\left(i\right)}}{{\sigma }_{1\text{в}}},\frac{{\sigma }_2^{\left(i\right)}}{{\sigma }_{2\text{в}}},\frac{{\tau }_{12}^{\left(i\right)}}{{\tau }_{12\text{в}}}\right\};$                    (9)
• – критерий сплошности материала $F_p^{\left(i\right)}={\left(\frac{{\sigma }_2^{\left(i\right)}}{{\sigma }_{2\text{в}}}\right)}^2+{\left(\frac{{\tau }_{12}^{\left(i\right)}}{{\tau }_{12\text{в}}}\right)}^2.$

Выполняя пошаговую процедуру пропорционального увеличения компонентов тензора контурных усилий $T_Z=f_Z\{T_{\xi }^{*},T_{\eta }^{*},T_{\xi \eta }^{*}\},$ для каждого шага вычисляют напряжения ${\sigma }_1^{\left(i\right)}, {\sigma }_2^{\left(i\right)}, {\sigma }_{11}^{\left(i\right)}$и значения коэффициента нагруженности $F_p^{\left(i\right)}$ по критериям (12). Предельная нагрузка соответствует выполнению условия прочности в виде $F_p^{\left(i\right)}(T_{\xi },T_{\eta },T_{\xi \eta })=1$.

На рис. 2 приведены результаты расчета критической T* и предельной T^разр нагрузок при сжатии и сдвиге в зависимости от граничных условий и структуры армирования слоев панели из углепластика, характеристики механических свойств которого указаны в табл. 1.

Таблица 1. Характеристики механических свойств углепластика 

Характеристика	вдоль волокон	поперек волокон
Модуль упругости, ГПа	110	17.5
Коэффициент Пуассона	0.260	0.041
Модуль сдвига в плоскости слоя, ГПа	7.5
Предел прочности при растяжении, МПа	1350	75
Предел прочности при сжатии, МПа	860	240
Предел прочности при сдвиге в плоскости слоя, МПа	56

Панель состоит из 20 ортотропных слоев со структурой армирования ± φ и направлением армирования γ = 20° и имеет четырехугольную форму с размерами сторон l₁ = l₂/cos χ_₁, l₂ = 2 м, b₁ = 3 м, b₂ = 2 м и углами скоса χ₁ = 30° и θ = 10°. Рассмотрены два варианта закрепления панели: (а) – шарнирное опирание по всем сторонам, (б) – смешанный вариант– левый и нижний края защемлены, правый – шарнирно оперт, верхний свободен.

Рис. 2. Расчетные значения критических ${\mathit{T}}_{\mathbf{\xi }}^{\mathbf{*}}$ (а) и предельных усилий ${\mathit{T}}_{\mathbf{разр}}^{\mathbf{\xi }}$ (б) при сжатии: 1 – шарнирное закрепление; 2 – смешанное закрепление. 

Расчеты показали, что предельная нагрузка углепластиковой панели значительно выше критических усилий потери устойчивости как при сжатии в осевом направлении (рис. 2а, б), так и при сдвиге (рис. 3а, б). Существенно отличается также характер зависимостей критических и предельных усилий от структуры армирования при сжатии для шарнирного закрепления. Критическая нагрузка потери устойчивости при шарнирном закреплении панели имеет наибольшее значение при φ ≈ ±45°, тогда как предельная нагрузка при таком армировании минимальна. При сжатии предпочтительна продольная укладка слоев, при сдвиге – поперечная (рис. 2 и 3).

Рис. 3. Влияние схемы закрепления панели на величину критических${\mathit{T}}_{\mathbf{\xi }}^{\mathbf{*}}$ – (1) и разрушающих ${\mathit{T}}_{\mathbf{р}\mathbf{а}\mathbf{з}\mathbf{р}}^{\mathbf{\xi }}$ – (2) нагрузок при сдвиге: (а) – шарнирное опирание; (б) – смешанное опирание.

Установлено, что влияние граничных условий имеет место только для предельных сжимающих усилий (рис. 3). При действии на панель касательных усилий (рис. 4) отмечено незначительное влияние схемы закрепления и совпадение величин предельных усилий в пределах точности расчетов во всем диапазоне углов укладки слоев. Предельные усилия при сжатии и сдвиге, найденные по критерию сплошности материала, заведомо выше усилий, найденных по другим критериям (9) (рис. 4а, б).

Рис. 4. Сравнение предельных усилий при сжатии 1 и сдвиге 2 по критерию сплошности материала: (а) – шарнирное опирание; (б) – смешанное опирание.

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

Финансирование. Данная работа финансировалась за счет средств бюджета Института машиноведения им. А. А. Благонравова РАН и Московского авиационного института. Никаких дополнительных грантов на проведение или руководство данным конкретным исследованием получено не было.

Конфликт интересов. Авторы заявляют об отсутствии конфликта интересов.