SOLUTIONS OF BOUNDARY VALUE PROBLEMS FOR ANISOTROPIC PLATES AND SHELLS BY BOUNDARY ELEMENTS METHOD

Full Text

1. Предварительные сведения

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

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

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

   Для решения поставленной задачи методом граничных элементов (МГЭ) необходимо предварительно определить фундаментальное решение задачи о плоском напряженном состоянии и изгибе анизотропной пластины.

   Для определения фундаментального решения задачи о плоском напряженном состоянии и изгибе анизотропной пластины и затем с помощью НМГЭ решения этих задач была использована универсальная методика расчета.

   Для отработки методики были решены несколько задач нелинейного деформирования анизотропных пластин и оболочек НМГЭ.

    

2. Постановка задачи

   Рассмотрим малые деформации тонкой линейно-упругой пологой оболочки, деформирование которой описывается моделью, основанной на гипотезах Кирхгофа — Лява в рамках теории среднего изгиба [1].

   Уравнения равновесия элемента пологой оболочки в усилиях имеют вид [2]

   ∂Tx + ∂Txy + p

   ∂Txy

   = 0;

   ∂Ty

   ∂x

   ∂2Mx

   ∂y x

   ∂2My

   +

   ∂x

   ∂2Mxy

   + py = 0,

   ∂y

   

   ∂x2 +

   

   

   ∂y2 + 2 ∂x∂y + kxTx + ky Ty + (1)

   ∂ ( ∂w

   

   

   + ∂x Tx ∂x

    

   + Txy

   ∂w )

   

   +

   ∂y

   ∂ ( ∂w

   

   

   ∂y Ty ∂y

    

   + Txy

   ∂w )

    

   

   ∂x

    

   + pz = 0,

   где kx, ky — главные кривизны оболочки в направлении осей x, y; px, py , pz — компоненты вектора внешней нагрузки, ось z направлена по нормали к центру кривизны.

   Усилия и моменты определяются через перемещения по формулам [3; 4]

   Tx = B11 εx + B12 εy + B16 γxy ; Ty = B12 εx + B22 εy + B26 γxy ; Txy = B16 εx + B26 εy + B66 γxy ;

    

    

   Mx = −

   ∂u

   1 (

   ∂w

   )2

   ∂v

   1 (

   ∂w

   )2

   ∂u

   ∂v

   ∂w ∂w

   ∂x

   — kxw +

   2

   ∂x

   ; εy =

   ∂y

   — ky w +

   2

   ∂y

   ; γxy =

   ∂y

   +

   ∂x

   +

   ; ∂x ∂y

    

   εx =

   D11 ∂ w + D

    

   ( 2

   ∂x2 12

   + 2D

    

   ∂2 w

   ∂y2

    

   ∂2 w )

   16 ∂x∂y

   (

    

   ; My = −

   ∂2w

    

   D12 ∂ w + D

    

   ( 2

   ∂x2 22

   ∂2w

    

   + 2D

    

   ∂2 w

   ∂y2

   ∂2w )

    

   ;

    

   ∂2 w )

   26 ∂x∂y

   (2)

   Mxy = −

   D16

   

   ∂x2 + D26

   

   

   ∂y2 + 2D66 ∂x∂y ,

   Великанов П.Г., Халитова Д.М. Решение задач нелинейного деформирования анизотропных пластин и оболочек...

   50 Velikanov P.G., Khalitova D.M. Solutions of boundary value problems for anisotropic plates and shelles by boundary...

    

   где Bij , Dij (i, j = 1, 2) – коэффициенты анизотропии [3]; (u, v, w) – компоненты вектора перемещения точки срединной поверхности; (εx, εy , γxy ) – тангенциальные деформации срединной поверхности.

   Если рассматривать пластины и оболочки постоянной толщины h = h0, то, подставив (2) в (1), получим систему нелинейных дифференциальных уравнений равновесия гибкой пологой анизотропной оболочки вида:

   L11 u(x, y) + L12 v(x, y) = l1(w(x, y)) − px;

   L21 u(x, y) + L22 v(x, y) = l2(w(x, y)) − py ; (3)

   L w(x, y) = l3(u(x, y), v(x, y), w(x, y)) + pz ,

   2

    

   где L11 = B11 ∂ + B

   2

    

   ∂ + 2B

   2

    

   2

    

   ∂ ; L12 = L21 = B16 ∂ + B

   2

    

   ∂ + (B

   2

    

   + B ) ∂ ;

   

   ∂x2

   2

   

   66 ∂y2

   2

   

   16 ∂x∂y

   2

   

   ∂x2

   

   26 ∂y2

   12 66

   

   ∂x∂y

   L22 = B66 ∂ 2 + B22 ∂ 2 + 2B26 ∂

    

   L = D11 ∂

    

   

   

   

   ∂x

   4

   + 4D

   ∂y ∂x∂y ;

   4

    

   4

    

   ∂ + 2(D12 + 2D66) ∂

    

   + 4D

    

   4

    

   ∂ + D

    

   4

    

   ∂ — линейные дифференциальные

   

   ∂x4

   

   16 ∂x3 ∂y

   

   ∂x2 ∂y2

   

   26 ∂x∂y3

   

   22 ∂y4

   

   операторы; li (i = 1, 3) — нелинейные дифференциальные операторы.

   В (3) были выделены слева линейные операторы задач о плоском напряженном состоянии и изгибе анизотропной пластины.

    

3. Определение фундаментального решения задачи изгиба и матрицы фундаментального решения задачи о плоском напряженном состоянии анизотропной пластины

   Пусть дана система линейных дифференциальных уравнений:

   L0 U (x, y) = F (x, y), (4)

   

   

   где L0 = [Llm ( d , d )] — заданный линейный дифференциальный оператор;

   dx dy

   U (x, y) = [ul(x, y)]T — векторная функция, подлежащая определению;

   F (x, y) = [fl(x, y)]T — заданная векторная функция правых частей; l, m = 1, N . Решение системы (4) U (x, y) представимо в виде свертки:

   ∫

   0

    

   U (x, y) = L−1 F (x, y) = G × F =

   Ω

   G(x − ξ, y − η) F (ξ, η) dξdη, ξ, η ∈ Ω,

   0

    

   где L−1

   — интегральный оператор, ядром которого является матрица фундаментального решения

   G(x − ξ, y − η) системы линейных дифференциальных уравнений; Ω — область определения дифференциального оператора L0. Матрица фундаментального решения G(x − ξ, y − η) определяется

   из выражения вида:

   L0 G(x − ξ, y − η) = δ(x − ξ, y − η) I, (5)

   где δ(x − ξ, y − η) — дельта-функция Дирака; I — единичная матрица размерностью N × N . Фундаментальные решения определяются с точностью до решения однородной системы уравнений

   L0 G(x − ξ, y − η) = 0 и, кроме того, являются обобщенными функциями. Из (5) видно, что фундаментальное решение G(x − ξ, y − η) зависит только от свойств дифференциального оператора L0.

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

   Для того чтобы свести задачу по поиску компонент матрицы фундаментального решения задачи о плоском напряженном состоянии анизотропной пластины к дифференциальному уравнению, подобному задаче изгиба анизотропной пластины, воспользуемся методикой, предложенной Л. Хермандером в [5; 6]. Согласно этой методике фундаментальное решение ищем в виде (via the ansatz [6]):

   0

    

   G(x − ξ, y − η) = L∗ ϑ(x − ξ, y − η), (6)

   0

    

   где L∗

   0

    

   = det(L0) L−1

   – ассоциированный к L0 дифференциальный оператор, элементами

   0

    

   которого являются алгебраические дополнения оператора LT

   (matrix of cofactors [6]). Таким образом,

   первоначальный и ассоциированный к нему дифференциальные операторы имеют вид

    

   L0 =

   ( L11 L12

   L21 L22

   )

   0

    

   , L∗ =

   ( L22 −L21 ) .

   −L12 L11

   Вестник Самарского университета. Естественнонаучная серия. 2021. Том 27, № 2. С. 48–61

   Vestnik of Samara University. Natural Science Series. 2021, vol. 27, no. 2, pp. 48–61 51

    

   В результате подстановки (6) в (5) приходим к уравнению относительно скалярной функции

   ϑ(x − ξ, y − η):

   det(L0) ϑ(x − ξ, y − η) = δ(x − ξ, y − η), (7)

   где det(L0) – определитель первоначального дифференциального оператора.

   Таким образом, определение матрицы фундаментального решения системы линейных дифференциальных уравнений данным методом предполагает: 1) вычисление компонент ассоциированного дифференциального оператора по методике, аналогичной методике вычисления компонент обратной матрицы; 2) решение уравнения (7), подобного уравнению задачи изгиба;

   0

    

   3) в соответствии с найденным ассоциированным дифференциальным оператором L∗

   проведение всех

   необходимых вычислений производных функции ϑ(x − ξ, y − η) по (6).

   Таким образом, необходимо найти фундаментальные решения двух подобных дифференциальных

   уравнений:

   [ ∂4 ∂4

   ∂4 ∂4

   ∂4 ]

   

   

   

   

   

   D11 ∂x4 + 4D16 ∂x3∂y + 2(D12 + 2D66) ∂x2∂y2 + 4D26 ∂x∂y3 + D22 ∂y4

   = δ(x − ξ, y − η),

   G(x − ξ, y − η) =

   [ ∂4 ∂4

   ∂4 ∂4

   ∂4 ]

   

   

   

   

   

   a22 ∂x4 − 2a26 ∂x3∂y + (2a12 + a66) ∂x2∂y2 − 2a16 ∂x∂y3 + a11 ∂y4

   ϑ(x − ξ, y − η) = δ(x − ξ, y − η),

   где aij (i, j = 1, 2) – упругие постоянные (коэффициенты деформации) [3].

   Таким образом, эти дифференциальные уравнения отличаются лишь коэффициентами в дифференциальных операторах, и их решение может быть найдено по аналогичной методике, например, фундаментальные решения (для разных вариантов корней характеристического уравнения) для задачи изгиба анизотропной пластины при действии нормальной к поверхности нагрузки q(x, y) имеют вид [3]. Без ограничения общности и в целях экономии места далее предполагаем, что точка приложения единичных сосредоточенных нагрузок, которые моделируются обобщенной δ-функцией Дирака, находится в начале координат:

    

   где

    

   ∂4 ∂4

   L[w(x, y)] = q(x, y),

    

   ∂4

    

   ∂4 ∂4

   

   

   

   

   

   L = D11 ∂x4 + 4D16 ∂x3∂y + 2(D12 + 2D66) ∂x2∂y2 + 4D26 ∂x∂y3 + D22 ∂y4

   — эллиптический дифференциальный оператор; x и y — декартовы прямоугольные координаты срединной поверхности; Dij — коэффициенты анизотропии [3]; w(x, y) — искомая функция прогиба пластины.

   Дифференциальное уравнение задачи изгиба бесконечной анизотропной пластины при действии нормальной к поверхности и приложенной в начале системы координат единичной сосредоточенной силы имеет вид:

   L [G(x, y)] = δ(x, y),

   где G(x, y) — фундаментальное решение, δ(x, y) — дельта-функция Дирака.

   Рассмотрим характеристическое уравнение:

   D11 µ4 + 4D16 µ3 + 2(D12 + 2D66) µ2 + 4D26 µ + D22 = 0.

   В настоящее время известны следующие фундаментальные решения задачи изгиба анизотропных пластин [7–12] (см. таблицу).

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

   Для нахождения нового точного фундаментального решения задачи изгиба анизотропных пластин (и, как следствие, матрицы фундаментального решения задачи о плоском напряженном состоянии анизотропных пластин) воспользуемся методом Леви, предварительно сделав краткое описание метода.

   Пусть имеется следующее дифференциальное уравнение эллиптического типа:

   Λ(u) + ∆(u) = 0,

   Великанов П.Г., Халитова Д.М. Решение задач нелинейного деформирования анизотропных пластин и оболочек...

   52 Velikanov P.G., Khalitova D.M. Solutions of boundary value problems for anisotropic plates and shelles by boundary...

    

    

   Информация об известных фундаментальных решениях задачи изгиба анизотропных пластин

   Information about known fundamental solutions to the problem of bending anisotropic plates

   Таблица

   Table

    

   Статьи	Решение для комплексно сопряженных корней характеристического уравнения	Решение для комплексно сопряженных кратных корней характеристического уравнения	Метод получения фундаментального решения	Примечание
   [7–9]	+	+	Не приведен	Решение приведено в полярной системе координат, что неудобно
   [10]	+	+	Интегральное преобразование Фурье	Хотя был использован один метод, но получены решения, отличающиеся на решение однородного уравнения
   [11]	+	+	Интегральное преобразование Фурье
   [12]	+	+	Метод последовательного интегрирования	-

    

    

   где операторы Λ и ∆ имеют вид:

    

   ∂2n ∑

    

   ∂l+m

   Λ = ∑

   l+m=2n

   alm(x, y)

    

   

   ∂xl∂ym

   , ∆ =

    

   l+m�2n−1

   blm(x, y)

   

   ∂xl∂ym .

   Обозначим через aj корни характеристического уравнения

   ∑

    

   l+m=2n

   alm(x, y)αl = 0.

   Корней будет ровно 2n, и все они будут комплексно сопряженными (для реальных коэффициентов анизотропии материалов при изгибе Dij (i, j = 1, 2)), действительными (для реальных упругих постоянных (коэффициентов деформации) материала в условиях плоского напряженного состояния aij (i, j = 1, 2)) или кратными.

   Составим из этих корней определитель Вандермонда

    α2n−1

   2n−2 

   1 α1 ... 1

   2n−1

   2n−2

   d =  α2 α2 ... 1  .

   

    

   

    

    ... ... ... ... 

   α2n−1

   2n−2

   2n α2n ... 1

   Алгебраические дополнения членов первого столбца обозначим dj .

   Введем в рассмотрение выражение

   σi(x − ξ, y − η) = αi(x − ξ) + (y − η)

   и составим функцию

    

   i 2n

   ∑

    

    

   j 2n−2

   ψ(x − ξ, y − η) = d

    

   j=1

   (−1) dj σj (x − ξ, y − η),

   которая является действительной и однозначной, имеющей особенность в точке ζ(ξ, η), что видно из самой функции ψ.

   Вестник Самарского университета. Естественнонаучная серия. 2021. Том 27, № 2. С. 48–61

   Vestnik of Samara University. Natural Science Series. 2021, vol. 27, no. 2, pp. 48–61 53

    

   Искомое фундаментальное решение задачи изгиба анизотропных пластин может быть получено после анализа особенностей функции ψ в виде:

   1

   

   G(x − ξ, y − η) = 2π(2n − 2)!a

    

   2n,0

   ψ(x − ξ, y − η).

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

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

   Для проверки правильности полученного фундаментального решения задачи изгиба и матрицы фундаментального решения задачи о плоском напряженном состоянии анизотропной пластины были использованы две методики: с помощью формулы дифференцирования обобщенных функций [13]

   ∂Φ ( ∂Φ )

   

   

   =

   ∫

   + (−1)i−1 δ(x1, ..., xn)

    

   Φ dx1...dxi

    

   1dxi+1...dxn,

   ∂xi

   ∂xi G −

   Γ

   

   ∂x

    

   где ( ∂Φ )

   i G

    

   — обычная производная от функции Φ; Γ —– граница области G (один из контуров, внутри

   которого находится особенность) и с помощью проверки равновесия пластины, ограниченной кривой,

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

   Найденные в представленной вашему вниманию работе фундаментальные решения были сравнены с решениями, ранее полученными другими методами. Так оказалось, что с точностью до решения однородного уравнения найденные фундаментальные решения совпадают с решениями, полученными в статьях [10–12]. Для случая ортотропного материала полученные фундаментальные решения с точностью до решения однородного уравнения совпадают с результатами, полученными в [13; 14], а для изотропного материала – с результатом, полученным в [10–15].

    

4. Непрямой метод граничных элементов (метод компенсирующих нагрузок) для решения задач нелинейного деформирования анизотропных пластин и оболочек

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

   На контуре Γ, ограничивающем оболочку, рассматриваются следующие граничные условия:

   ∂w

   w = 0,

   

   ∂n = 0, u = 0, v = 0 − жесткая заделка,

   w = 0, Mn = 0, u = 0, v = 0 − шарнирное закрепление, где Mn = Mx cos2 α + My sin2 α + Mxy sin 2α — изгибающий момент на контуре Γ.

   Процесс последовательных приближений для решения системы (3) принимается в виде

   uk+1 = u(k) + αu(u − u(k)); vk+1 = v(k) + αv (v − v(k)); wk+1 = w(k) + αw (w − w(k)); (k = 0, 1, 2, ...), (8)

   где αu, αv , αw (0 < αu :( 1, 0 < αv :( 1, 0 < αw :( 1) – параметры, обеспечивающие сходимость итерационного процесса (параметры релаксации). Вектор U (u, v, w) определяется из решения системы линейных уравнений, описывающих изгиб и растяжение анизотропной пластины

   L11 u(x, y) + L12 v(x, y) = l1(w(k)) − px;

   L21 u(x, y) + L22 v(x, y) = l2(w(k)) − py ; (9)

   L w(x, y) = l3(u(k), v(k), w(k)) + pz .

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

   Великанов П.Г., Халитова Д.М. Решение задач нелинейного деформирования анизотропных пластин и оболочек...

   54 Velikanov P.G., Khalitova D.M. Solutions of boundary value problems for anisotropic plates and shelles by boundary...

    

   По методу компенсирующих нагрузок решение (9) ищется в виде:

   ∫

   u(t) =

   Γ

   ∫

   v(t) =

   Γ

   [G11(t, ζ) φ1(ζ) + G12(t, ζ) φ2(ζ)] ds(ζ) + ur (t);

    

   [G21(t, ζ) φ1(ζ) + G22(t, ζ) φ2(ζ)] ds(ζ) + vr (t); (10)

   ∫

   w(t) =

   Γ

   [

   G(t, ζ) q(ζ) −

   ∂G(t, ζ) ]

   m(t, ζ)

   ∂n(ζ)

   ds(ζ) + wr (t).

   Здесь t(x, y) ∈ Ω, ζ(ξ, η) ∈ Γ; ds – элемент длины контура Γ; Gij (t, ζ), (i, j = 1, 2) – компоненты

   матрицы фундаментального решения задачи о плоском напряженном состоянии анизотропной пластины;

   G(t, ζ) – фундаментальное решение задачи изгиба анизотропной пластины; n(ζ) – внешняя нормаль к контуру Γ; φ1(ζ), φ2(ζ), q(ζ), m(ζ) – компоненты вектора компенсирующих нагрузок на контуре Γ (φ1, φ2 – усилия в срединной поверхности, направленные вдоль координатных осей x, y; q, m – усилие нормальное срединной поверхности и момент вокруг касательной к контуру Γ).

   Функции ur (t), vr (t), wr (t) определяются соотношениями

   ∫

   ur (t) =

   Ω

   [G11(t, ζ) l1(w(k)(ζ)) + G12(t, ζ) l2(w(k)(ζ))]

    

   dΩ(ζ);

   vr (t) =

   ∫ [ ]

   G21(t, ζ) l1(w(k)(ζ)) + G22(t, ζ) l2(w(k)(ζ))

   Ω

   ∫

    

   dΩ(ζ); (11)

   l3

    

   wr (t) =

   Ω

   G(t, ζ) [

   (u(k)(ζ), v(k)(ζ), w(k)(ζ)) + p(ζ)] dΩ(ζ).

   Разрешающая система сингулярных интегральных уравнений с неизвестными компенсирующими нагрузками получается при подстановке (10) в граничные условия на контуре Γ и после анализа предельных значений потенциалов имеет вид:

   1. жесткая заделка

      ∫ ∫

      G(t, ζ)q(ζ)ds(ζ) −

      Γ Γ

      ∂G(t, ζ)

      ∂n1

      m(ζ)ds(ζ) + wr (t) = 0,

      ∫ ∂G(t, ζ)

       

      

      ∂n

      Γ

      ∫

      ∫

      q(ζ)ds(ζ) −

      Γ

      ∂2G(t, ζ)

       

      

      ∂n∂n1

       

      m(ζ)ds(ζ) +

      ∂wr (t)

       

      

      ∂n

       

      = 0,

       

      Γ

      ∫

       

      Γ

   2. шарнирное закрепление

   [G11(t, ζ)φ1(ζ) + G12(t, ζ)φ2(ζ)] ds(ζ) + ur (t) = 0,

    

   [G21(t, ζ)φ1(ζ) + G22(t, ζ)φ2(ζ)] ds(ζ) + vr (t) = 0;

   ∫ ∫

   G(t, ζ) q(ζ) ds(ζ) −

   Γ Γ

   ∂G(t, ζ)

   ∂n1

   m(ζ) ds(ζ) + wr (t) = 0,

   m(t) ∫

   

   +

   2

   Γ

   ∫

   L1(G(t, ζ)) q(ζ) ds(ζ) −

   Γ

   ∫

   ( ∂G(t, ζ))

   

   L1

   ∂n1

   n

    

   m ds(ζ) + M r (t) = 0,

   [G11(t, ζ) φ1(ζ) + G12(t, ζ) φ2(ζ)] ds(ζ) + ur (t) = 0,

   Γ

   ∫

   [G21(t, ζ) φ1(ζ) + G22(t, ζ) φ2(ζ)] ds(ζ) + vr (t) = 0.

   Γ

   Перемещения и усилия в области Γ и на контуре Γ определяются соотношениями вида (10).

   Ядра граничных интегральных уравнений содержат особенности типа ln r, 1/r, 1/r2 при r → 0.

   Интегралы с особенностью типа 1/r2 определяются в смысле конечных значений по Адамару.

   Вестник Самарского университета. Естественнонаучная серия. 2021. Том 27, № 2. С. 48–61

   Vestnik of Samara University. Natural Science Series. 2021, vol. 27, no. 2, pp. 48–61 55

    

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

   Функции ur , vr , wr определяются интегральными операторами со слабой особенностью (интегрирование проводится по срединной поверхности Ω). Для вычисления этих интегралов срединная поверхность разбивается на треугольники или секторы (рис. 1). Далее каждый треугольник разбивается на отдельные элементы (ячейки). На срединную поверхность Ω наносится сеточная область, узлы которой являются центрами тяжести отдельных ячеек. Интегралы по отдельным ячейкам можно представить в виде

   ∫

   J = k(ti, ζ) φ(ζ) dΩ(ζ), ζ ∈ ∆j (j = 1, 2, ..., m),

   ∆j

   α

    

   где ti — центр тяжести i-й ячейки; ∆j — площадь j-й ячейки; m — число ячеек в области интегрирования; k(ti, ζ) — ядро интегрального оператора; φ(ζ) — плотность, которая удовлетворяет условию Гельдера

   |φ(ζ) − φ(t)| :( C |ζ − t|

   ; 0 < α :( 1; C = const > 0.

    

   

    

   Рис. 1. Методика вычисления частных решений Fig. 1. Technique for calculating particular solutions

    

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

   Площадь дугового сектора и положение центра тяжести определяются следующими соотношениями:

   

   

   

   

   

   1 1 sin α (2 + cos2 α ) − α cos α

   

   2

    

   F = R2 (α − sin α); yc = 2R 3

   2 2 2 2 .

   α − sin α

   Если i ̸= j, то интеграл J вычисляется по формуле

   J = k(ti, tj ) φ(tj ) ∆j .

   При i = j эти интегралы имеют слабую особенность и вычисляются численно по квадратурным формулам Гаусса.

   Условие окончания итерационного процесса (9) принимается в следующем виде:

   11 11

   11U − U (k)11

   

   11 11

   < ε,

   ||

    

   где ||U 2

   11U (k)11

   = ∑[u2(ti) + v2(ti) + w2(ti)] – принятая норма; ε – малая положительная величина; U (k) =

   i

   = (u(k), v(k), w(k)).

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

   Великанов П.Г., Халитова Д.М. Решение задач нелинейного деформирования анизотропных пластин и оболочек...

   56 Velikanov P.G., Khalitova D.M. Solutions of boundary value problems for anisotropic plates and shelles by boundary...

    

   значений выбранного параметра. Эффективность алгоритма зависит от способа выбора этого параметра. Вопросом применения метода продолжения по параметру посвящены работы [16; 17].

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

   оболочек проводится с помощью зависимостей “прогиб – нагрузка”. За ведущий параметр принимался прогиб w∗ в заданной точке t∗(x∗, y∗) срединной поверхности оболочки.

   В данном случае к системе нелинейных уравнений добавляется уравнение w(t∗) = w∗ и нагрузка p

   считается неизвестной величиной. Строится итерационный процесс с шагом по прогибу ∆w∗ в заданной

   точке t∗, т. е. решаются нелинейные задачи для значений прогибов w∗, w∗, ..., w∗

   в точке t∗. При этом

   w∗ ∗ ∗

   0 1 N

   i = wi−1 + ∆w

   (i = 1, 2, ..., N ).

   i

    

   За начальное приближение к решению при значении прогиба w∗

   в точке t∗ можно принять

    

   где β = w

    

   ∗

   i

   

   w∗

   i−1

   u(i) = βu(i−1), v(i) = βv(i−1), w(i) = βw(i−1), (i = 1, 2, ..., N ),

   Задание вышеприведенного начального приближения является эффективным, если деформация

   оболочки происходит без резкого изменения формы.

   0

    

   При значении в точке t∗ прогиба w∗

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

   u(0) = 0, v(0) = 0, w(0) = 0.

   Это соответствует решению на первой итерации задачи изгиба и растяжения анизотропной пластины в линейной постановке.

    

5. Сопоставление ортотропного и анизотропного материала

Рассмотрим ортотропный материал (рис. 2), главные направления упругости (физические оси x′ и y′) которого не совпадают с его геометрическими осями (x и y). Между ними угол φ.

 

Рис. 2. Oртотропный материал Fig. 2. Orthotropic material

 

В результате дифференциальный оператор ортотропного материала

∂4 ∂4 ∂4

L = D1 ∂x4 + 2D3 ∂x2∂y2 + D2 ∂y4 − ортотропный оператор задачи изгиба;

в вышеприведенных обстоятельствах после пересчета коэффициентов жесткости (Di(i = 1, 3))

′

 

D11 = D1 cos4 φ + 2D3 sin2 φ cos2 φ + D2 sin4 φ,

′

 

′

 

D22 = D1 sin4 φ + 2D3 sin2 φ cos2 φ + D2 cos4 φ, D66 = Dk + (D1 + D2 − 2D3) sin2 φ cos2 φ,

1

D

 

γ1 = ′

22

[D2ν1 + (D1 + D2 − 2D3) sin2

φ cos2

φ],

′ 1 2 2

D16 = 2 (D2 sin

φ − D1 cos

φ + D3 cos 2φ) sin 2φ,

′ 1 2 2

D26 = 2 (D2 cos

φ − D1 sin

φ − D3 cos 2φ) sin 2φ,

Вестник Самарского университета. Естественнонаучная серия. 2021. Том 27, № 2. С. 48–61

Vestnik of Samara University. Natural Science Series. 2021, vol. 27, no. 2, pp. 48–61 57

 

примет вид, аналогичный анизотропному материалу

′ ∂4

′ ∂4

′ ′ ∂4

′ ∂4

′ ∂4

L = D11 ∂x′ 4 + 4D16 ∂x′ 3∂y′ + 2(D12 + 2D66) ∂x′ 2∂y′ 2 + 4D26 ∂x′ ∂y′ 3 + D22 ∂y′ 4

• анизотропный оператор

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

   

  Пример 1

  Круглая пластина (рис. 3) с жестко заделанными краями находится под действием равномерно распределенной поперечной нагрузки интенсивности q (φ = 300).

  На рис. 3 представлены зависимости “максимальный прогиб — нагрузка” в центре пластины, полученные НМГЭ, на базе метода R-функций [18] и результаты работы [19]. Сплошная линия соответствует решению [18], пунктирная (с длинным пунктиром) — решению [19], а пунктирная (с коротким пунктиром и точками) — решению НМГЭ.

  При решении задачи НМГЭ контур пластины был равномерно разбит на 20 одинаковых по длине дуговых элементов. Релаксационные параметры и невязка приняты в виде:

  αu = αv = αw = 0.2; ε = 0.001. Итерационный процесс сходился за 4–13 итераций.

   

  Исследование проведено на примере пластин, изготовленных из материалов со следующими механическими характеристиками [18]:

  1. боропластик (кривая BR): E1/E2 = 10;

     G˜/E2 = 1/3; ν1 = 0.22;

  2. графит (кривая GR): E1/E2 = 40;

     G˜/E2 = 0.6; ν1 = 0.25;

  3. стеклопластик (кривая GL): E1/E2 = 3;

  Размеры пластины: R = 0.1(м); h = 0.01(м).

  G˜/E2 = 0.6; ν1 = 0.25;.

   

  

   

  Рис. 3. Результаты вычислений для круглой пластины Fig. 3. Calculation results for a round plate

   

  Пример 2

  Квадратная в плане пологая сферическая оболочка (рис. 4) с шарнирно закрепленными краями с размером стороны 2a находится под действием равномерно распределенного нормального давления интенсивности p, направленного к центру кривизны оболочки (φ = 300).

  На рис. 4 приведены зависимости “максимальный прогиб — нагрузка” для оболочки с параметром

  2 4

  

  

  

  Rh

   

  кривизны kx = ky = k = (2a)

  

  

  Eh4

   

  = 16. Параметр нагрузки p = p (2a)

  . Контур разбивался на 28 элементов.

  Длина стороны: a = 0.1(м). Параметры релаксации принимались равными: αu = αv = 0.3; αw = 0.1 при

  Великанов П.Г., Халитова Д.М. Решение задач нелинейного деформирования анизотропных пластин и оболочек...

  58 Velikanov P.G., Khalitova D.M. Solutions of boundary value problems for anisotropic plates and shelles by boundary...

   

  точности ε = 10−3. Итерационный процесс сходился за 5–68 итераций. Механические характеристики материала аналогичны приведенным в примере 1.

  

  

  

  

  h

   

  На рис. 4 приведены зависимости q − w, где w = w

• относительный прогиб в центре оболочек с

кривизнами k = 0; 4; 8; 16. При k = 0; 4 зависимости q − w носят монотонный характер. При k = 8; 16,

как видно из рис. 4, оболочка теряет устойчивость хлопком. Например, верхняя и нижняя критические

 

нагрузки для оболочки с кривизной k = 16 соответственно равны

q = 115.73; q = 40.4.

 

 

Рис. 4. Результаты вычислений для квадратной в плане пологой сферической оболочки Fig. 4. Calculation results for a shallow spherical shell that is square in plan

 

Выводы

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

 

About the authors

P. G. Velikanov

Kazan (Volga Region) Federal University

Email: pvelikanov@mail.ru
ORCID iD: 0000-0003-0845-2880

Candidate of Physical and Mathematical Sciences, assistant professor of the Department of Theoretical Mechanics

Russian Federation, 18, Kremlevskaya Street, Kazan, 420008, Russian Federation

D. M. Khalitova

Kazan (Volga Region) Federal University

Author for correspondence.
Email: diana982000@gmail.com
ORCID iD: 0000-0002-2239-9222

Master's degree student

Russian Federation, 18, Kremlevskaya Street, Kazan, 420008, Russian Federation

References

Supplementary files