USE OF THE INDIRECT BOUNDARY ELEMENTS METHOD FOR ISOTROPIC PLATES ON AN ELASTIC WINKLER BASE AND PASTERNAK — VLASOV BASE

Full Text

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

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

   теории малого изгиба [1].

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

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

    

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

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

   αt

   D∇2∇2W = pz (x, y) − p(W ) − (1 + ν)

   ∆T,

   h

   где ∇2 — оператор Лапласа; W — прогиб точки срединной поверхности пластины; pz — интенсивность

   нормального давления, действующего на пластину; p — реактивное давление, которое зависит от

   упругого основания, на котором лежит пластина; ν — коэффициент Пуассона (коэффициент поперечной деформации); αt — коэффициент линейного температурного расширения; h — толщина пластины; ∆T — изменение температуры между поверхностями пластины.

   Рассмотрим упругое основание, которое характеризуется следующей функцией:

   p(W ) = K1W + K2W 2 + K3W 3 − g1

   ∂2W

    

   

   ∂x2

   — g2

   ∂2W

   

   ∂y2 ,

   где Ki (i = 1, 2, 3) и gi (i = 1, 2) — коэффициенты упругого основания. В таблице приведено несколько возможных вариантов оснований.

    

    

   Наиболее распространенные варианты упругих оснований The most common options for elastic foundations

   Таблица

   Table

    

   K1 ̸= 0 ; K2 = K3 = g1 = g2 = 0	Линейное изотропное основание Винклера (Циммермана — Фусса)
   K1 ̸= 0; K2 = K3 = 0; g1 = g2 = g ̸= 0	Линейное основание Пастернака — Власова
   K1 ̸= 0; K2 = K3 = 0; g1 ̸= g2 ̸= 0	Линейное ортотропное основание
   K1 ̸= 0; K2 ̸= 0; K3 ̸= 0	Нелинейное основание

    

   В дальнейшем будем рассматривать линейные упругие основания Винклера и Пастернака — Власова как самые распространенные.

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

   Vestnik of Samara University. Natural Science Series. 2021, vol. 27, no. 2, pp. 33–47 35

    

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

   Пусть дана система линейных дифференциальных уравнений в частных производных [5]:

   L0U (x, y) = F (x, y), (1)

   

   

   где L0 = [Llm( ∂ , ∂ )] — исходный линейный дифференциальный оператор в частных производных,

   ∂x ∂y

   

   

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

   Решение системы (1) U (x, y) представимо в виде свертки:

   ∫

   0

    

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

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

   Ω

   0

    

   где L−1

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

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

     дифференциального оператора L0 .

     Матрица фундаментального решения G(x − ξ, y − η) определяется из выражения вида:

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

     L0G(x − ξ, y − η) = 0 и, кроме того, являются обобщенными функциями.

     Из (2) видно, что фундаментальное решение G(x − ξ, y − η) зависит только от свойств

     дифференциального оператора L0.

     В дальнейшем фундаментальное решение G(x−ξ, y −η) будем обозначать G(t, ζ) = G(t−ζ), где t(x, y)

     и ζ(ξ, η).

     Элементы матрицы G(t, ζ) определятся из решения систем дифференциальных уравнений вида:

      

     где δlm — символ Кронекера.

     N

     

     ∑ Lmk Gkl(t, ζ) = δlmδ(t − ζ), (l, m = 1, N ), (3)

     k=1

     Решение систем (3) осуществляется, например, методом двумерного интегрального преобразования

     Фурье [1–4]. Без ограничения общности и в целях экономии места далее предполагаем, что точка приложения единичных сосредоточенных нагрузок, которые моделируются обобщенной δ-функцией Дирака, находится в начале координат.

     

     Двумерная трансформанта (образ, изображение) f (ξ, η) функции f (x, y) определяется выражением вида:

     ∫

      

     +∞ +∞

     1 ∫

     

     F [f (x, y)] = f (ξ, η) =

     2π

     

     f (x, y)ei(ξx+ηy)dxdy, i = √−1.

     −∞ −∞

     

     Соотношение между трансформантами производных от f (x, y) и самой трансформантой f (ξ, η) может быть определено следующим выражением:

     ∂(n+m)f (x, y)

     

     

     F [ ∂xn∂ym

     ](ξ, η) = (−iξ)n(−iη)m f (ξ, η).

     Трансформанта двумерной обобщенной δ-функции Дирака, которая, в силу присущих ей свойств, представима в виде δ(x, y) = δ(x)δ(y), определяется из следующего соотношения: F [δ(x, y)] =

     +∞ +∞

     =

      

     1 ∫ ∫

     

     2π

     

     2π

      

     δ(x, y)ei(ξx+ηy)dxdy = 1 .

     −∞ −∞

     Трансформанты для производных δ-функции Дирака примут вид:

     

     − −

      

     F [δ(m+n)(x, y)] = 1 ( iξ)m( iη)n.

     2π

     Системы (3) в пространстве трансформант примут вид:

     N

     

     ∑ Lmk (−iξ, −iη)Gkl(ξ, η) =

     k=1

     

     1

     

     2π δlm, (l, m = 1, N ).

     Великанов П.Г., Куканов Н.И., Халитова Д.М. Использование непрямого метода граничных элементов...

     36 Velikanov P.G., Kukanov N.I., Khalitova D.M. Use of the indirect boundary elements method for isotropic plates...

      

     Формула обращения двумерного интегрального преобразования Фурье (прообраз, оригинал) имеет

     +∞ +∞

     

     2π

      

     вид: f (x, y) = 1 ∫

     ∫ f (ξ, η)e−i(ξx+ηy)dξdη.

     −∞ −∞

     При вычислении несобственных интегралов будем пользоваться свойством:

     +∞ +∞

     ∂(n+m)f (x, y)

     

     1 ∫ ∫ n m

     i(ξx+ηy)

     

     

     ∂xn∂ym = 2π

      

     −∞ −∞

     (−iξ)

     (−iη)

     f (ξ, η)e−

     dξdη.

     Для вычисления интегралов будем использовать теорию вычетов из комплексного анализа: если

     Q(z) — функция, где z — комплексное число, удовлетворяющая условиям:

     1. аналитическая при Im z ): 0 всюду, за исключением числа полюсов, лежащих на действительной оси;

     2. не имеет полюсов на действительной оси;

     3. при z → ∞ функция zQ(z) → 0 равномерно при 0 � arg z � π;

     4. интегралы

     ∞

      

     ∫ Q(x)dx и

     0

     0

     ∫ Q(x)dx являются сходящимися, то

     −∞

     +∞

     ∫

      

     Кроме того,

     Q(x)dx = 2πi ∑

     −∞ Imz�0

      

     ∫∞

     Res Q(z).

     Q(x) cos mx dx = πi ∑

     0 Imz�0

     Res Q(z) eim z ;

     ∞

      

     ∫ Q(x) sin mx dx = π ∑

     0 Imz�0

      

     Res Q(z) eim z.

     Если z = a — полюс функции Q(z) порядка m, то вычет функции Q(z) определяется по формуле:

     1

     

     Res Q(a) =

      

     lim

     

     dm−1

     [(z − a)m Q(z)].

     (m − 1)! z→a dzm−1

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

     Рассмотрим вычисление несобственных интегралов следующих типов:

     +∞ +∞

     1 ∫

     

     I1(x, y) = 2π

     ∫ 1

     

     (ξ2 + η2) e−

     i(ξx+ηy)

     dξdη,

     −∞ −∞

      

     +∞ +∞

     1 ∫

     

     I2(x, y) = 2π

     ∫ 1

     

     (ξ2 + η2 + λ2) e−

     i(ξx+ηy)dξdη.

     −∞ −∞

     Интеграл I1(x, y) определим, воспользовавшись свойством интегралов от четных и нечетных функций на интервале, симметричном относительно начала координат:

     +∞ +∞

     ∫∞ ∫∞ ∫∞

     1 ∫

     ∫ 1

     i(ξx+ηy) 2

     cos ξxdξ

     1 ηx

     I1(x, y) = 2π

     −∞ −∞

     (ξ2 + η2) e−

     dξdη =

     π

     0

     cos ηydη

     0

     (ξ2 + η2) = η e−

     0

     cos ηydη,

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

     ∫∞ cos ξx

     eizx

     eiz x 1 π 1

     

     

     1

      

     ∑ 1 1 ηx

     (ξ2 + η2) dξ = πi

     Res = πi

     (z2 + η2) 2z

     = e− .

     2 η

     0 Imz�0

     1 1z1 =iη

     Внешний интеграл вычислим с помощью теории обобщенных функций [4]. Рассмотрим обобщенную функцию вида:

     Θ(η) = ln η+ =

     { 0, η � 0,

     ln η, η > 0.

     

     η

      

     Производная этой функции есть обобщенная функция Θ′(η) = 1 , η > 0.

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

     Vestnik of Samara University. Natural Science Series. 2021, vol. 27, no. 2, pp. 33–47 37

      

     Согласно правилу дифференцирования обобщенных функций, имеем

     ∫

     (Θ′, φ) = −(Θ, φ′) = −

     φ′(η)Θ(η)d(η),

      

     где φ(η) – основная функция (бесконечное число раз дифференцируемая и в бесконечно удаленной точке сама функция и ее бесконечное число производных стремятся к нулю).

     Здесь в качестве основной функции выступает функция φ(η) = e−ηx cos ηy. Учитывая все вышесказанное, получим:

      

     I1(x, y) =

     ∫∞ 1

     e−ηx

     η

     0

     ∫∞

     cos ηydη = x

     0

     1 2

      

     e−ηx

      

     2

      

     cos ηy ln η dη + y

     ∞

      

     ∫ e−ηx

     0

      

     sin ηy ln ηdη =

     

     = − 2 ln(x

     + y ) − C = −(ln r + C),

     где входящие интегралы являются табличными; r2 = x2 + y2, C = 0.5772157 ... – постоянная Эйлера.

     Рассмотрим интеграл I2(x, y):

     +∞ +∞

     ∫∞ ∫∞

     1 ∫

     

     I2(x, y) = 2π

     ∫ 1

     

     (ξ2 + η2 + λ2) e−

     

     i(ξx+ηy)dξdη = 2

     π

     cos ηydη

     cos ξxdξ

     

     (ξ2 + η2 + λ2) .

     −∞ −∞ 0 0

     Аналогично рассмотренному выше

     ∫∞ cos ξx π 1

     отсюда

     

     −x√η2 +λ2

     

     (ξ2 + η2 + λ2) dξ = 2 √η2 + λ2 e ,

     0

      

     I2(x, y) =

     

     ∞

      

     ∫ 1 −x√η2 +λ2

     √η2 + λ2 e

     0

      

     cos ηy dη = K0(λr),

     где K0(λr) — модифицированная функция Бесселя второго рода нулевого порядка; этот интеграл также является табличным.

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

      

     3.1. Нахождение фундаментального решения с помощью метода последовательного интегрирования

     Запишем, например, уравнение изгиба изотропной пластины на упругом основании Винклера в виде:

     ∂4 ∂4

     ( + 2

     4

      

     + ∂ + χ4)W (x, y) = q(x, y)

     q(x, y)

     ( L W (x, y) = ),

     

     ∂x4

     

     ∂x2∂y2

     

     

     ∂y4

     

     D 0 D

     

     

     3

      

     где χ4 = K1 ; D = Eh

   • изгибная (цилиндрическая) жесткость пластины.

   D 12(1−ν2 )

   Представим дифференциальный оператор изгиба пластины на упругом основании L0 в виде двух

   последовательных операторов Гельмгольца:

   ∂4 ∂4 ∂4 4

   ∂2 ∂2

   2 ∂2 ∂2 2

   

   

   

   

   

   

   

   L0 = ( ∂x4 + 2 ∂x2∂y2 + ∂y4 + χ ) = ( ∂x2 + ∂y2 − iχ )( ∂x2 + ∂y2 + iχ ).

   Фундаментальное решение ищем в виде:

    

   L0G1(x, y) =

   δ(x, y)

   

   .

   D

   После введения тривиального обозначения получим

   ∂2

   

   ( ∂x2 +

   ∂2

    

   

   ∂y2

   − iχ2) G˜1(x, y) =

   δ(x, y)

   

   .

   D

   Далее к этому дифференциальному уравнению последовательно применяем двумерное интегральное преобразование Фурье, выделение трансформанты G˜1(x, y), применение формулы обращения двумерного интегрального преобразования Фурье. В результате G˜1(x, y) примет вид:

   

   1

    

   G˜1(x, y) = − 2πD K0(

   

   √iχr).

   Великанов П.Г., Куканов Н.И., Халитова Д.М. Использование непрямого метода граничных элементов...

   38 Velikanov P.G., Kukanov N.I., Khalitova D.M. Use of the indirect boundary elements method for isotropic plates...

    

   Далее по известной величине G˜1(x, y) восстанавливаем G1(x, y):

   ∫∞

   i G1(x, y) = − 2D

   0

   

   

   rK0(

   0

    

   √iχ(r0)H(2)(√

   iχ(r − r0))dr0,

   0

    

   где H(2)

   — функция Ханкеля 0-го порядка 2-го рода.

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

   ∂4 ∂4 ∂4 4

   ∂2 ∂2

   2 ∂2 ∂2 2

   

   

   

   

   

   

   

   L0 = ( ∂x4 + 2 ∂x2∂y2 + ∂y4 + χ ) = ( ∂x2 + ∂y2 + iχ )( ∂x2 + ∂y2 − iχ ),

   ∂2

   

   ( ∂x2

   ∂2

   

   + ∂y2

   + iχ2)G˜2(x, y) =

   i

   δ(x, y)

   

   ,

   D

   

   G˜2(x, y) =

   ∫∞

   

   0

    

   H(2)(√iχr),

   4D

   i G2(x, y) = − 2D

   0

   

   0

    

   rH(2)(√

   iχ(r0)K0(

   √

   

   iχ(r − r0))dr0.

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

   l2 r

    

   

   где χ4 = 1

   

   

   G(x, y) = G1(x, y) + G2(x, y) = − 2πD kei0( l ),

   

   

   = K1 ; kei0( r ) — функция Томпсона-Кельвина.

   l4 D l

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

   пластины на упругом основании Винклера рассмотрим еще 2 варианта его нахождения:

   а) трансформанта Фурье фундаментального решения задачи изгиба изотропной пластины на упругом основании Винклера имеет вид:

   1

   

   G(ξ, η) =

   1 i 1

   

   

   = (

   1

   

   − ).

   2πD (ξ2 + η2)2 + χ4

   4πDχ2

   ξ2 + η2 + iχ2

   ξ2 + η2 − iχ2

   Восстанавливая его оригинал с точностью до решения однородного уравнения, получим

   G(x, y) =

   i

   

   4πDχ2

   

   K0(√iχr) −

   1

   

   8Dχ2

   

   0

    

   H(2)(√iχr) = −

   l2

   

   2πD

   

   r kei0( l );

   б) трансформанта Фурье фундаментального решения задачи изгиба изотропной пластины на упругом основании Винклера имеет вид:

   1

   

   G(ξ, η) =

   1 1 1

   

   

   = − Im .

   2πD (ξ2 + η2)2 + χ4

   Восстанавливая оригинал, получим:

   2πDχ2

   ξ2 + η2 + iχ2

   1

   

   √ l2 r

   G(x, y) = − 2πDχ2 Im(K0(χr

   

   i)) = − 2πD kei0( l ).

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

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

   ∂Φ

   

   ∂xi

   ∂Φ

   

   = (

   ∂xi

   ∫

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

   Γ

    

   Φdx1 . . . dxi

    

   −1dx

    

   i+1

    

   . . . dxn,

   

   ∂xi

    

   где ( ∂Φ )G — обычная производная от функции Φ, Γ — граница области G (один из контуров, внутри которого находится особенность) и с помощью проверки равновесия пластины на упругом основании Винклера, ограниченной кривой, при действии на нее единичной нагрузки. Результаты проверок подтверждают правильность найденного фундаментального решения задачи изгиба изотропной пластины на упругом основании Винклера.

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

   Vestnik of Samara University. Natural Science Series. 2021, vol. 27, no. 2, pp. 33–47 39

    

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

    

   L0W (x, y) =

   pz (x, y) , D

   где L0 = ∇2 −2p2∆+χ4 — дифференциальный оператор; kz, kt — параметры упругого основания (первый и второй коэффициенты постели или коэффициенты сжатия и сдвига соответственно);

   2

   

   2p2 = g

   kth

   

   =

   

   

   

   , χ4 = K1 = kz = 1 .

   D 4D

   D D l4

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

   

   4πs2 D

    

   других, которые получаются из них как частные случаи: a) p > χ ⇒ G(x, y) = − 1

   

   [K0(√p2 + s2r) −

   

   

   − K0(√p2 − s2r)], s2 = √p4 − χ4;

   2

   

   

   b) p < χ ⇒ G(x, y) = − 1

   

   

   kei0(r/l, φ), β2 = √χ4 − p4, φ = 1 arctanβ ;

   2πβ2 D

   

   4πDχ2

    

   c) p = χ ⇒ G(x, y) = − 1

   

   4πp2 D

    

   d) χ = 0 ⇒ G(x, y) = − 1

    

   0

    

   kei′ (r/l, 0);

   

   [ln r + K0(√2pr)];

   2 p2

   

   2

    

   e) p = 0 ⇒ G(x, y) = − l

   kei0(r/l);

   2πD

   

   8πD

    

   f) p = χ = 0 ⇒ G(x, y) = 1

    

   r2ln r.

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

   Для реально существующих материалов пластин, лежащих на упругом основании Пастернака — Власова, имеет место случай b):

   1

   

   G(x, y) = − 2πβ2D kei0(r/l, φ),

   где kei0(x, φ) — обобщенная функция Томпсона-Кельвина (мнимая часть обобщенной модифицированной функции Бесселя второго рода нулевого порядка K0).

    

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

   Рассмотрим тонкую линейно-упругую пластину, срединная поверхность которой занимает область Ω+, ограниченную гладким кусочно-ляпуновским контуром Γ (рис. 1). Дополним область Ω+ областью Ω− до бесконечной области. По контуру Γ к бесконечной пластине приложим компенсирующие нагрузки q(ζ), m(ζ). Нагрузки q(ζ) и m(ζ) – распределенные по контуру Γ усилие, нормальное к поверхности пластины, и момент вокруг касательной к контуру Γ соответственно.

    

   

    

   Рис. 1. Тонкая линейно-упругая пластина Fig. 1. Thin linear elastic plate

   Великанов П.Г., Куканов Н.И., Халитова Д.М. Использование непрямого метода граничных элементов...

   40 Velikanov P.G., Kukanov N.I., Khalitova D.M. Use of the indirect boundary elements method for isotropic plates...

    

   Решение задачи изгиба многосвязной пластины рассматривается как сумма основного (частного) и компенсирующего решений [6]:

   M

    

   W (t) = Wr (t) + ∑[∫

   i=1 Γi

   G(t, ζ)qi(ζ) −

   i

    

   ∂G(t, ζ)

   

   m (ζ)ds(ζ)], (4)

   ∂n1

   где Wr (t) = ∫

   Ω+

   n

   G(t, ζ)pz (ζ)dΩ(ζ) + ∑ G(t, ζi)Pi(ζi).

   i=1

   Здесь t(x, y) — точка внутренней области Ω+; ζ(ξ, η)) — точка контура Γ; G(t, ζ) — фундаментальные

   решения задачи изгиба пластины, лежащей на упругом основании Винклера или Пастернака — Власова; pz (ζ) – интенсивность нормального давления, действующего на пластину; Pi — модуль сосредоточенной силы, приложенной в точке ζi(i = 1, 2, ..., n).

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

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

   ∂W

   −

    

   W = 0, W,n =

   

   = 0 жесткая заделка, (5)

   ∂n

   W = 0, Mn = 0 − шарнирное опирание, (6)

   Vn = k1W, Mn = k2W,n −упругое закрепление. (7)

   Здесь Mn, Vn — изгибающий момент и обобщенная поперечная сила на контуре Γ; k1, k2 - постоянные, определяющие упругие свойства закрепления (при k1 = k2 = 0 получим граничные условия свободного края).

   Усилия и моменты на контуре пластины в глобальной декартовой системе координат X, Y

   определяются по формулам [4]:

    

   Mnτ = −

   Mn = Mx cos2 α + My sin2 α + Mxy sin 2α;

   MT = Mx sin2 α + My cos2 α − Mxy sin 2α;

   Mx − My sin 2α + M cos 2α; Q = Q cos α + Q sin α;

   2 xy n x y

   ∂Mnτ

   Vn = Qn + ∂s ,

   где для задачи изгиба изотропной пластины соотношения примут вид:

   Mx = −D(W, xx +νW, yy ); My = −D(W, yy +νW, xx ); Mxy = −D(1 − ν)W, xy ;

   ∂ ∂

   

   

   Qx = −D ∂x ∆W = −D(W, xxx +W, xyy ); Qy = −D ∂y ∆W = −D(W, yyy +W, xxy ).

   Для изотропной пластины, лежащей на упругом основании, моменты и усилия на контуре пластины можно также определить в локальной системе координат n, τ , построенной в точке t(x, y) контура, где записываются граничные условия. При этом ось n — внешняя нормаль к контуру в этой точке; а ось τ — касательная к контуру в этой точке.

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

   Mn = −D(∆W − (1 − ν)W, ττ ); Mτ = −D(∆W − (1 − ν)W, nn );

   ∂M nτ

   Mnτ = −D(1 − ν))W, nτ ; Qn = −D∆W,n Vn = Qn + ∂s ,

   2 2

   

   ∂n2

    

   +

    

   где ∆ = ∇2 = ∂

   ∂

   

   ∂τ 2

   • оператор Лапласа в локальной системе координат.

     

     ∂s

      

     При определении производной по дуговой координате контура ∂M nτ

     следует учитывать, что α =

     

     

     ∂s

      

     = ±

      

     ρ′

      

     = α(s); k = ∂α 1

   • кривизна контура в рассматриваемой точке (алгебраическая величина), ρ′ –

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

   ∂ ∂3W

   

   

   Vn = −D( ∂n ∆W + (1 − ν)[ ∂n∂τ 2 − k(2W,xy sin 2α + [W,xx −W,xy ] cos 2α)]) =

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

   Vestnik of Samara University. Natural Science Series. 2021, vol. 27, no. 2, pp. 33–47 41

    

   ∂ ∂3W

   ∂2W

   ∂2W

   

   

   

   = −D( ∂n ∆W + (1 − ν)[ ∂n∂τ 2 − k( ∂n2 −

   

   ∂τ 2 )]).

   Компенсирующие нагрузки определяются из граничных условий (5)–(7) на контуре пластины, которые записываются на контуре Γ для области Ω+.

   Известно, что выражение (4) в точности удовлетворяет дифференциальным уравнениям задачи изгиба изотропной пластины на упругом основании Винклера и Пастернака — Власова. Чтобы

   (4) было решением краевой задачи, необходимо из системы сингулярных интегральных уравнений, которая получается при подстановке (4) в краевые условия (5)–(7) на контуре пластины, определить компенсирующие нагрузки q(ζ), m(ζ). При этом нужно совершить предельный переход точки t(x, y) из внутренней области Ω+ на границу Γ. Тогда ядра интегральных уравнений будут определены во всех точках контура, за исключением точки, где t(x, y) = ζ(ξ, η). В этой точке ядра потенциалов будут иметь особенности.

   В области Ω+ изгибающие и крутящий моменты, а также поперечные силы для изотропной пластины, лежащей на упругом основании, по известным компенсирующим нагрузкам q(ζ), m(ζ), определяются по формулам:

   ∫ ∂2G

   ∂2G

   ∂3G

   ∂3G r

   Mx(t) = −D

   Γ

   

   

   

   −

    

   [( + ν )q(ζ) (

   ∂x2 ∂y2 ∂x2∂n1

   

   + ν

   ∂y2∂n1

   )m(ζ)]ds(ζ) + Mx (t);

   ∫ ∂2G

   ∂2G

   ∂3G

   ∂3G r

   My (t) = −D

   Γ

   

   

   

   −

    

   [( + ν )q(ζ) (

   ∂y2 ∂x2 ∂y2∂n1

   

   + ν

   ∂x2∂n1

   )m(ζ)]ds(ζ) + My (t);

   ∫ ∂2G

   ∂3G r

   Mxy (t) = −D(1 − ν)

   Γ

   

   

   −

    

   [ q(ζ)

   ∂x∂y ∂x∂y∂n1

   m(ζ)]ds(ζ) + Mxy (t);

   ∫ ∂∆G

   ∂2∆G r

   Qx(t) = −D [

   Γ

   

   −

    

   q(ζ)

   ∂x ∂x∂n1

   m(ζ)]ds(ζ) + Qx(t);

   ∫ ∂∆G

   ∂2∆G r

    

   где

   Qy (t) = −D [

   Γ

   

   −

    

   q(ζ)

   ∂y ∂y∂n1

   m(ζ)]ds(ζ) + Qy (t),

   2 r 2 r

   2 r 2 r 2 r

   ∂ W (t)

   Mr

   ∂ W (t) r

   ∂ W (t)

   ∂ W (t) r

   ∂ W (t)

   x (t) = −D(

   ∂x2 + ν

   ∂y2 ); My (t) = −D(

   ∂y2 + ν

   ∂x2 ); Mxy (t) = −D(1 − ν)

   ;

   ∂x∂y

   3 r 3 r

   3 r 3 r

   ∂ W (t)

   Qr

   ∂ W (t) r

   ∂ W (t)

   ∂ W (t)

   x(t) = −D(

   ∂x3 +

   ∂x∂y2 ); Qy (t) = −D(

   ∂x2∂y +

   ∂y3 ).

    

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

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

   С учетом предельных значений потенциалов по методике [7], а также результатов дифференцирования фундаментального решения задачи изгиба пластины на упругом основании запишем сингулярные интегральные уравнения, из решения которых определяются компенсирующие нагрузки q(ζ), m(ζ) для различных граничных условий. При этом предельные значения потенциалов берутся в области Ω+.

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

      ∫ ∫

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

      Γ Γ

      ∂G(t, ζ)

      ∂n1

      m(ζ)ds(ζ) + W

      r (t) = 0,

      ∫ ∂G(t, ζ)

       

      

      ∂n

      Γ

      ∫

      q(ζ)ds(ζ) −

      Γ

      ∂2G(t, ζ)

       

      

      ∂n∂n1

       

      m(ζ)ds(ζ) +

      ∂Wr (t)

       

      

      ∂n

       

      = 0.

      Великанов П.Г., Куканов Н.И., Халитова Д.М. Использование непрямого метода граничных элементов...

      42 Velikanov P.G., Kukanov N.I., Khalitova D.M. Use of the indirect boundary elements method for isotropic plates...

       

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

      ∫

       

      Γ

       

      ∫

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

      Γ

       

      ∂G(t, ζ)

      ∂n1

       

      m(ζ)ds(ζ) + W

       

      r (t) = 0,

      1 ∫

      

      2 m(t) − D

      Γ

      [∆G(t, ζ) − (1 − ν)

      ∂2G(t, ζ)

       

      

      ∂τ 2

      ]q(ζ)ds(ζ)+

       

   3. Свободный край

   ∫

   +D [

   Γ

   ∂∆G(t, ζ)

    

   

   ∂n1

   − (1 − ν)

   ∂3G(t, ζ)

    

   

   ∂τ 2∂n1

   n

    

   ]m(ζ)ds(ζ) + Mr (t) = 0,

   1 ∫

   

   2 m(t) − D

   Γ

   [∆G(t, ζ) − (1 − ν)

   ∂2G(t, ζ)

    

   

   ∂τ 2

   ∫

   ]q(ζ)ds(ζ) + D

   Γ

   ∂∆G(t, ζ)

   

   [

   ∂n1

   − (1 − ν)

   ∂3G(t, ζ)

    

   

   ∂τ 2∂n1

    

   ]m(ζ)ds(ζ)+

   1 ∫

   

   2 q(t) − D

   Γ

    

   ∂∆G(t, ζ)

   

   [

   ∂n

    

   + (1 − ν)[

    

   ∂3G(t, ζ)

    

   

   ∂n∂τ 2

   n

    

   +Mr (t) = 0,

   ∂2G(t, ζ)

   

   − k(t)( ∂n2 −

    

   ∂2G(t, ζ)

    

   

   ∂τ 2

   ∫

   )]]q(ζ)ds(ζ) + D

   Γ

    

   ∂2∆G(t, ζ)

   

   [ +

   ∂n∂n1

   ∂4G(t, ζ)

   

   1

    

   +(1 − ν)[ ∂n∂τ 2∂n

   ∂3G(t, ζ)

   

   1

    

   − k(t)( ∂n2∂n

   ∂3G(t, ζ)

   

   − ∂τ 2∂n1

   n

    

   )]]m(ζ)ds(ζ) + V r (t) = 0,

   где введены обозначения

   Mr r r

   n(t) = −D[W,nn (t) + νW,ττ (t)];

   3 r 2 r 2 r

   V r r

   ∂

    

   ∂ W (t)

   ∂ W (t)

   ∂ W (t)

   n (t) = −D( ∂n ∆W

   (t) + (1 − ν)[

   ∂n∂τ 2 − k(t)(

   ∂n2 −

   ∂τ 2 )]).

   Функции плотностей q(ζ) и m(ζ) удовлетворяют условию Гельдера:

   |q(t) − q(ζ)| � C|t − ζ|α; |m(t) − m(ζ)| � C|t − ζ|α,

   где C = const > 0, 0 < α � 1.

   Внеинтегральные члены q(t)/2, m(t)/2 в интегральных уравнениях получаются на основе анализа предельных значений потенциалов [7].

    

6. Численная реализация

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

   N

   ΓN =

   ∑ Γj и ΓN → Γ при N → ∞ или max |Γj | → 0, где N – количество элементов контура ΓN .

   j=1 j

   В пределах элемента компенсирующие нагрузки q(ζ) и m(ζ) считаются постоянными. Точки коллокации берутся в узлах, которые располагаются в серединах элементов Γj .

   Для решения систем сингулярных интегральных уравнений применяется метод механических квадратур [4]. При удовлетворении граничным условиям в узлах, после вычисления интегралов, решение задачи сводится к системе линейных алгебраических уравнений вида:

   N N

   

   ∑ q(ζj )Gij − ∑ m(ζj )Hij + Ai = 0, 1, N,

   j=1

   N

   j=1

   N

   

   ∑ q(ζj )Rij − ∑ m(ζj )Mij + Bi = 0, 1, N,

   j=1

   j=1

   где коэффициенты Gij, Hij, Rij, Mij в зависимости от вида граничных условий принимают соответствующие выражения и представляют собой интегралы по элементу Γj ; Ai, Bi — правые части системы уравнений, принимающие соответствующие выражения в зависимости от граничных условий; q(ζj ) = qi, m(ζj ) = mj — искомые значения компенсирующих нагрузок, принимающие постоянные значения в пределах каждого элемента Γj . Система состоит из 2N уравнений с 2N неизвестными, т. е. является замкнутой и имеет единственное решение.

   Интегралы Gij, Hij, Rij, Mij при i ̸= j не имеют особенностей и вычисляются по восьмиузловой

   квадратурной формуле Гаусса. При i = j ядра интегралов имеют особенности типа ln r, 1/r, 1/r2

   при r → 0 и, кроме того, необходимо учесть внеинтегральные члены уравнений. Такие сингулярные

   интегралы по Γj вычисляются аналитически.

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

   Vestnik of Samara University. Natural Science Series. 2021, vol. 27, no. 2, pp. 33–47 43

    

7. Вычисление сингулярных интегралов по элементам контура

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

   особенности при x → ξ, y → η. Считаем, что компенсирующие нагрузки в пределах этих элементов

   постоянны, поэтому они входят в интегральные уравнения в виде множителей перед интегралами от

   фундаментального решения и его производных.

    

   

    

   Рис. 2. Прямолинейный граничный элемент Fig. 2. Straight boundary element

    

   Пусть Γj (j = 1, 2, ..N ) — один из элементов контура Γ, представляющий собой отрезок прямой (рис. 2), соединяющей точки (xA, yA) и (xB, yB ), а t(x, y) — середина отрезка Γj :

   x = (xA + xB )/2; y = (yA + yB )/2.

   Введем координату t0, определяющую положение точки ζ(ξ, η) на Γj :

    

   где t0 ∈ [−1; +1].

   ξ = (xA + xB )

   2

   + (xB − xA)

   2

    

   t0; η =

   (yA + yB ) +

   2

   (yB − yA)

   2

    

   t0,

   Расстояние между точками t и ζ равно

   

   

   

   2

    

   B

    

   r = √(x − ξ)2 + (y − η)2 = |t0| √(x

    

   • xA

      

     )2 + (yB

      

   • yA

   

   )2 = |t0| δ,

   2

   где δ — длина отрезка Γj .

   Фундаментальное решение и его производные зависят от аргументов x − ξ, y − η, которые

   определяются выражениями

   x − ξ = −

   0 − −

    

   (xB − xA) t ; y η = 2

   (yB − yA) t .

   2 0

   Выразим элемент длины ds через координату t0

   

   ds = √ ξ′ )2 + (η′ )2dt0 =

   

   (

    

   √ (xB − xA)

    

   )2 + (

   (yB − yA)

   δ

   )2dt0 =

    

   dt0.

   ( t0 t0

   2 2 2

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

    

8. Результаты решения некоторых задач

Пример 1

Кольцевая пластина с внешним радиусом R = 0.1 м и внутренним радиусом r = 0.05 м лежит на упругом основании Пастернака — Власова с коэффициентами постели kz = 15400 MH/м3, kt =

= 10000 МН/м3 (рис. 3). Внутренний контур шарнирно оперт, а внешний — жестко защемлен.

На пластину действует постоянная равномерно распределенная поперечная нагрузка pz = 1 МПа. Физические параметры: h = 0.01 м, E = 2 · 105 МПа, ν = 0.3. При численном решении каждый контур

k

 

разбивался на 40 элементов, а частное решение принято в виде Wr (t) = pz . На рис. 4, а, б приведены

z

результаты решения поставленной задачи: графики распределения прогибов W , углов поворота Wr ,

радиального Mr и тангенциального Mφ моментов по радиусу пластины. Точки соответствуют численному решению НМГЭ, а непрерывные линии — точному решению задачи, которое ищется в виде:

r r r r r

W (r) = C1 ber( l , φ) + C2 bei( l , φ) + C3 ker( l , φ) + C4 kei( l , φ) + W

(r),

Великанов П.Г., Куканов Н.И., Халитова Д.М. Использование непрямого метода граничных элементов...

44 Velikanov P.G., Kukanov N.I., Khalitova D.M. Use of the indirect boundary elements method for isotropic plates...

 

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

 

 

Рис. 3. Кольцевая пластина для расчета Fig. 3. Ring plate for calculation

 

 

a б

 

Рис. 4. Результаты расчета Fig. 4. Calculation results

 

Пример 2

 

 

 

Рис. 5. Пологая сферическая оболочка Fig. 5. Shallow spherical shell

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

Vestnik of Samara University. Natural Science Series. 2021, vol. 27, no. 2, pp. 33–47 45

 

Квадратная в плане пологая изотропная сферическая оболочка радиуса R с угловым вырезом размером a с граничными условиями шарнирного закрепления (пунктирная линия), жесткой заделки (наклонная линия) и свободного края (нет линии) (рис. 5) лежит на упругом основании Винклера. Оболочка находится под действием постоянного температурного поля, которое по толщине пластины изменяется по линейному закону, и сосредоточенной силы P , приложенной в начале координат. Перепад

температуры по толщине пластины ∆T = 100 град. Коэффициент температурного расширения материала примем следующим: αt = 0.125 · 10−4 град−1. Геометрические и механические параметры пластины: a = 0.1 м, h = 0.01 м, R = 0.5 м, Kz = 15400 МН/м3, P = 4000 Н, E = 2 · 105 МПа, ν = 0.3.

Контур разбит на 60 элементов, по 10 на каждую сторону.

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

2 2 kth2 2 1 Eh

pz (x, y)

αt 2

∇ ∇ W −

4D ∇ W + D (kz + R2 )W =

D − (1 + ν) h ∇ ∆T,

конкретизируя его на случай основания Винклера (kt = 0).

На рис. 6 представлены результаты решения поставленной задачи: на рис. 6, а прогиб w, а на рис. 6, б напряжения (как нормальные, так и касательные) σx, σy, τxy при z = h/2 в сечении пластины x = 0. На графиках напряжений сплошная линия соответствует σx; пунктирные линии (линия с большей длиной пунктира и линия с меньшей длиной пунктира) – соответственно σy и τxy .

 

 

a б

 

Рис. 6. Результаты расчета Fig. 6. Calculation results

 

Пример 3

 

 

a б

 

Рис. 7. Двусвязная пластина сложной формы Fig. 7. Doubly connected plate of complex shape

Великанов П.Г., Куканов Н.И., Халитова Д.М. Использование непрямого метода граничных элементов...

46 Velikanov P.G., Kukanov N.I., Khalitova D.M. Use of the indirect boundary elements method for isotropic plates...

 

Двусвязная пластина сложной формы с различными граничными условиями на отдельных участках лежит на упругом основании и находится под действием термомеханического нагружения (сосредоточенной силы P , приложенной в начале координат и линейного по толщине температурного

поля) (рис. 7, а). Параметры пластины и поля: h = 0.01 м, a = 0.1 м, E = 2 · 105 МПа, ν = 0.3, kz =

= 15400 МН/м3, P = 5000 Н, ∆T = 100 град, αt = 0.125 · 10−4 град−1.

На рис. 7, б приведена поверхность распределения прогиба W с линиями уровня.

 

Выводы

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

About the authors

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

N. I. Kukanov

Ulyanovsk State Technical University

Email: kukanov_n_i@mail.ru
ORCID iD: 0000-0003-0880-4591

Candidate of Physical and Mathematical Sciences, assistant professor of the Department
of Industrial and Civil Engineering

Russian Federation, 32, Severny Venets Street, Ulyanovsk, 432027, 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

Supplementary files