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

## Abstract

Modern mechanical engineering sets the tasks of calculating thin-walled structures that combine lightness and economy on the one hand and high strength and reliability on the other. In this regard, the use of anisotropic materials and plastics seems justified. The problems of the theory of plates and shells belong to the class of boundary value problems, the analytical solution of which, due to various circumstances (nonlinearity of differential equations, complexity of geometry and boundary conditions, etc.), cannot be
determined. Numerical methods help to solve this problem. Among numerical methods, undeservedly little attention is paid to the boundary element method. In this regard, the further development of indirect method of compensating loads for solving problems of the anisotropic plates and shells theory based on the application
of exact fundamental solutions is relevant.
The paper considers the application of the indirect boundary element method for solving of an anisotropic plates and shells nonlinear deformation problem. Since the kernels of the system of singular integral equations to which the solution of the problem is reduced are expressed in terms of the fundamental solution and its
derivatives, first of all, the article provides a method for determining the fundamental solutions to the problem of bending and the plane stress state of an anisotropic plate. The displacement vector is determined from the solution of linear equations system describing the bending and plane stress state of an anisotropic plate. The solution of the system is performed by the method of compensating loads, according to which the area representing the plan of the shallow shell is supplemented to an infinite plane, and on the contour that limits the area, compensating loads are applied to the infinite plate. Integral equations of indirect BEM are given. In this paper, the study of nonlinear deformation of anisotropic plates and shallow shells is carried out using the “deflection – load” dependencies. The deflection at a given point on the median surface of the shell was taken as the leading parameter.

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) = L1 F (x, y) = G × F =

G(x ξ, y η) F (ξ, η) dξdη, ξ, η ,

0

где L1

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

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) L1

– ассоциированный к 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�2n1

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)) или кратными.

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

α2n1

2n2

1 α1 ... 1

2n1

2n2

d = α2 α2 ... 1 .

... ... ... ...

α2n1

2n2

2n α2n ... 1

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

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

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

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

i 2n

j 2n2

ψ(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)i1 δ(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 = wi1 + ∆w

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

i

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

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

где β = w

i

w

i1

u(i) = βu(i1), v(i) = βv(i1), w(i) = βw(i1), (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 x4 + 4D16 x3∂y+ 2(D12 + 2D66) x2∂y2 + 4D26 x∂y3 + D22 ∂y4

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

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

### Пример 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...

точности ε = 103. Итерационный процесс сходился за 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

## Выводы

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

×

### P. G. Velikanov

Kazan (Volga Region) Federal University

Author for correspondence.
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

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

Master's degree student

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

## References

1. Galimov K.Z. Fundamentals of the nonlinear theory of thin shells. Kazan: Izd-vo Kazanskogo un-ta, 1975, 326 p. (In Russ.)
2. Volmir A.S. Nonlinear dynamics of plates and shells. Moscow: Nauka, 1972, 432 p. Available at: https://dwg.ru/lib/1965. (In Russ.)
3. Lekhnitsky S.G. Anisotropic plates. Moscow: OGIZ-Gostekhizdat, 1947, 355 p. Available at: https://lib-bkm.ru/12625. (In Russ.)
4. Kornishin M.S. Nonlinear problems of the plates and shallow shells theory and methods of their solution. Moscow: Nauka, 1964, 192 p. (In Russ.)
5. Hormander (должна быть О с двумя точками наверху) L. The Analysis of Linear Partial Differential Operators I Distribution Theory and Fourier Analysis. Moscow: Mir, 1986, 464 p. Available at: http://www.bookre.org/reader?file=442772. (In Russ.)
6. Shanz M., Antes H. A Boundary Integral Formulation for the Dynamic Behavior of a Timoshenko Beam. Electronic Journal of Boundary Elements, 2002, vol. BETEQ, no. 3, pp. 348–359. Available at: https://www.researchgate.net/publication/228746551_A_boundary_integral_formulation_for_the_dynamic_behavior_of_a_Timoshenko_beam.
7. Wu B.C., Altiero N.J. A new numerical method for the analysis of anisotropic thin plate bending problems. Computer Methods in Applied Mechanics and Engineering, 1981, vol. 25, issue 3, pp. 343–353. DOI: https://doi.org/10.1016/0045-7825(81)90037-2.
8. Shi G., Bezine G. A general boundary integral formulation for the anisotropic plate bending problems. Journal of Composite Materials, 1988, vol. 22, pp. 694–716. DOI: http://doi.org/10.1177/002199838802200801.
9. Albuquerque E.L., Sollero P., Venturini W.S., Aliabadi M.H. Boundary element analysis of anisotropic Kirchhoff plates. International Journal of Solids and Structures, July 2006, vol. 43, issues 14–15, pp. 4029–4046. DOI: http://doi.org/10.1016/j.ijsolstr.2006.03.027.
10. Guryanov I.N., Artyukhin Yu.P. Fundamental solution of the plane problem and the bending plates problem of an anisotropic body. Kazan: KGU, 1994, 33 p. (In Russ.)
11. Levada V.S. Construction of a fundamental solution for the bending an anisotropic plate problem. Pridneprovskii nauchnyi vestnik, 1996, no. 4, p. 8. (In Russ.)
12. Velikanov P.G., Artyukhin Yu.P., Kukanov N.I. Bending of an anisotropic plate by the method of boundary elements. Actual problems of continuum mechanics-2020. Kazan: Kazanskii universitet – izd-vo Akademii nauk RT, 2020, pp. 105–111. (In Russ.)
13. Artyukhin Yu.P., Gribov A.P. Solving problems of nonlinear deformation of plates and shallow shells by the method of boundary elements. Kazan: Fen, 2002, 199 p. Available at: https://booksee.org/book/438702. (In Russ.)
14. Gribov A.P., Velikanov P.G. Application of the Fourier transformation to obtain a fundamental solution to the bending an orthotropic plate problem. In: Proceedings of the All-Russian Scientific Conference (26-28 May 2004). Part 3. Matem. Mod. Kraev. Zadachi. Samara: SamGTU, 2004, pp. 67–71. Available at: https://elibrary.ru/item.asp?id=7556596; http://mi.mathnet.ru/mmkz191 (In Russ.)
15. Crouch S.I., Starfield S. Boundary element in solid mechanics. Moscow: Mir, 1987, 328 p. Available at: https://booksee.org/book/438609. (In Russ.)
16. Vorovich I.N. Mathematical problems of the nonlinear theory of shallow shells. Moscow: Nauka, 1989, 376 p. Available at: https://booksee.org/book/449596. (In Russ.)
17. Grigolyuk E.I., Shalashilin V.I. Problems of nonlinear deformation. Moscow: Nauka, 1988, 232 p. Available at: https://booksee.org/book/438580. (In Russ.)
18. Rvachev V.L., Kurpa L.V., Khomenko M.M., Bolotina A.Yu. Application of the theory of R-functions to the calculation of flexible orthotropic flat shells of complex shape in plan. In: Proceedings of the XV All-Union Conference on the Theory of Shells and Plates, vol. 1. Kazan: Izdatel’stvo Kazanskogo universiteta, 1990, pp. 342–347. (In Russ.)
19. Chia C.Y. Nonlinear analysis of plates. New York: Mc. Graw-Hill International Book Co., 1980, XIV, 422 p.

Copyright (c) 2022 Velikanov P.G., Khalitova D.M.

You consent to our cookies if you continue to use our website.