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


Cite item

Abstract

The calculation tasks combining lightness, economy, high strength and reliability of thin-walled structures on an elastic base are relevant for modern mechanical engineering. In this regard, the use of isotropic materials on an elastic base seems justified, therefore their calculation is considered in this article. 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 (the nonlinearity of differential equations, the 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 compensating loads method for solving problems of the theory of isotropic plates on an elastic base of Winkler and Pasternak — Vlasov, based on the application of exact fundamental
solutions, is relevant.

Full Text

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

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

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

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

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

     

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

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

    αt

    D22W = pz (x, y) p(W ) (1 + ν)

    T,

    h

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

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

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

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

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

    2W

     

    image

    ∂x2

    g2

    2W

    image

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

    image

    image

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

    ∂x ∂y

    image

    image

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

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

    0

     

    где L1

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

      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

      image

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

      k=1

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

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

      image

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

       

      ++

      1

      image

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

      2π

      image

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

      −∞ −∞

      image

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

      (n+m)f (x, y)

      image

      image

      F [ ∂xn∂ym

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

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

      ++

      =

       

      1 ∫ ∫

      image

      2π

      image

      2π

       

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

      −∞ −∞

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

      image

      − −

       

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

      2π

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

      N

      image

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

      k=1

      image

      1

      image

      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...

       

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

      ++

      image

      2π

       

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

      f (ξ, η)ei(ξx+ηy)dξdη.

      −∞ −∞

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

      ++

      (n+m)f (x, y)

      image

      1 ∫ ∫ n m

      i(ξx+ηy)

      image

      image

      ∂xn∂ym = 2π

       

      −∞ −∞

      ()

      ()

      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

      image

      Res Q(a) =

       

      lim

      image

      dm1

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

      (m 1)! za dzm1

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

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

      ++

      1

      image

      I1(x, y) = 2π

      1

      image

      (ξ2 + η2) e

      i(ξx+ηy)

      dξdη,

      −∞ −∞

       

      ++

      1

      image

      I2(x, y) = 2π

      1

      image

      (ξ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

      image

      image

      1

       

      1 1 ηx

      (ξ2 + η2)= πi

      Res = πi

      (z2 + η2) 2z

      = e.

      2 η

      0 Imz�0

      1 1z1 =

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

      Θ(η) = ln η+ =

      { 0, η 0,

      ln η, η > 0.

      image

      η

       

      Производная этой функции есть обобщенная функция Θ(η) = 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η =

      image

      = 2 ln(x

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

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

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

      ++

      1

      image

      I2(x, y) = 2π

      1

      image

      (ξ2 + η2 + λ2) e

      image

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

      π

      cos ηydη

      cos ξxdξ

      image

      (ξ2 + η2 + λ2) .

      −∞ −∞ 0 0

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

      cos ξx π 1

      отсюда

      image

      xη2 +λ2

      image

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

      0

       

      I2(x, y) =

      image

       

      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) = ),

      image

      ∂x4

      image

      ∂x2∂y2

      image

      image

      ∂y4

      image

      D 0 D

      image

      image

      3

       

      где χ4 = K1 ; D = Eh

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

    D 12(1ν2 )

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

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

    4 4 4 4

    2 2

    2 2 2 2

    image

    image

    image

    image

    image

    image

    image

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

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

     

    L0G1(x, y) =

    δ(x, y)

    image

    .

    D

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

    2

    image

    ( x2 +

    2

     

    image

    ∂y2

    2) G˜1(x, y) =

    δ(x, y)

    image

    .

    D

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

    image

    1

     

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

    image

    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

    image

    image

    rK0(

    0

     

    (r0)H(2)(

    (r r0))dr0,

    0

     

    где H(2)

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

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

    4 4 4 4

    2 2

    2 2 2 2

    image

    image

    image

    image

    image

    image

    image

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

    2

    image

    ( ∂x2

    2

    image

    + ∂y2

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

    i

    δ(x, y)

    image

    ,

    D

    image

    G˜2(x, y) =

    image

    0

     

    H(2)(iχr),

    4D

    i G2(x, y) = 2D

    0

    image

    0

     

    rH(2)(

    (r0)K0(

    image

    (r r0))dr0.

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

    l2 r

     

    image

    где χ4 = 1

    image

    image

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

    image

    image

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

    l4 D l

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

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

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

    1

    image

    G(ξ, η) =

    1 i 1

    image

    image

    = (

    1

    image

    ).

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

    4πDχ2

    ξ2 + η2 + 2

    ξ2 + η2 2

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

    G(x, y) =

    i

    image

    4πDχ2

    image

    K0(iχr)

    1

    image

    82

    image

    0

     

    H(2)(iχr) =

    l2

    image

    2πD

    image

    r kei0( l );

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

    1

    image

    G(ξ, η) =

    1 1 1

    image

    image

    = Im .

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

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

    2πDχ2

    ξ2 + η2 + 2

    1

    image

    l2 r

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

    image

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

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

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

    Φ

    image

    ∂xi

    Φ

    image

    = (

    ∂xi

    )G + (1)i1δ(x1, . . . , xn)

    Γ

     

    Φdx1 . . . dxi

     

    1dx

     

    i+1

     

    . . . dxn,

    image

    ∂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

    image

    2p2 = g

    kth

    image

    =

    image

    image

    image

    , χ4 = K1 = kz = 1 .

    D 4D

    D D l4

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

    image

    4πs2 D

     

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

    image

    [K0(p2 + s2r)

    image

    image

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

    2

    image

    image

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

    image

    image

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

    2πβ2 D

    image

    4πDχ2

     

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

    image

    4πp2 D

     

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

     

    0

     

    kei (r/l, 0);

    image

    [ln r + K0(2pr)];

    2 p2

    image

    2

     

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

    kei0(r/l);

    2πD

    image

    8πD

     

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

     

    r2ln r.

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

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

    1

    image

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

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

     

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

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

     

    image

     

    Рис. 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, ζ)

    image

    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 =

    image

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

    ∂n

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

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

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

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

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

     

    M =

    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

    ∂M

    Vn = Qn + s ,

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

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

    ∂ ∂

    image

    image

    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

    M = D(1 ν))W, ; Qn = DW,n Vn = Qn + s ,

    2 2

    image

    ∂n2

     

    +

     

    где ∆ = 2 =

    image

    ∂τ 2

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

      image

      ∂s

       

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

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

      image

      image

      ∂s

       

      = ±

       

      ρ

       

      = α(s); k = α 1

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

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

    ∂ ∂3W

    image

    image

    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

    image

    image

    image

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

    image

    ∂τ 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

    Γ

    image

    image

    image

     

    [( + ν )q(ζ) (

    ∂x2 ∂y2 ∂x2∂n1

    image

    + ν

    ∂y2∂n1

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

    2G

    2G

    3G

    3G r

    My (t) = D

    Γ

    image

    image

    image

     

    [( + ν )q(ζ) (

    ∂y2 ∂x2 ∂y2∂n1

    image

    + ν

    ∂x2∂n1

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

    2G

    3G r

    Mxy (t) = D(1 ν)

    Γ

    image

    image

     

    [ q(ζ)

    ∂x∂y ∂x∂y∂n1

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

    G

    2G r

    Qx(t) = D [

    Γ

    image

     

    q(ζ)

    ∂x ∂x∂n1

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

    G

    2G r

     

    где

    Qy (t) = D [

    Γ

    image

     

    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, ζ)

       

      image

      ∂n

      Γ

      q(ζ)ds(ζ)

      Γ

      2G(t, ζ)

       

      image

      ∂n∂n1

       

      m(ζ)ds(ζ) +

      ∂Wr (t)

       

      image

      ∂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

      image

      2 m(t) D

      Γ

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

      2G(t, ζ)

       

      image

      ∂τ 2

      ]q(ζ)ds(ζ)+

       

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

    +D [

    Γ

    G(t, ζ)

     

    image

    ∂n1

    (1 ν)

    3G(t, ζ)

     

    image

    ∂τ 2∂n1

    n

     

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

    1

    image

    2 m(t) D

    Γ

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

    2G(t, ζ)

     

    image

    ∂τ 2

    ]q(ζ)ds(ζ) + D

    Γ

    G(t, ζ)

    image

    [

    ∂n1

    (1 ν)

    3G(t, ζ)

     

    image

    ∂τ 2∂n1

     

    ]m(ζ)ds(ζ)+

    1

    image

    2 q(t) D

    Γ

     

    G(t, ζ)

    image

    [

    ∂n

     

    + (1 ν)[

     

    3G(t, ζ)

     

    image

    ∂n∂τ 2

    n

     

    +Mr (t) = 0,

    2G(t, ζ)

    image

    k(t)( n2

     

    2G(t, ζ)

     

    image

    ∂τ 2

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

    Γ

     

    2G(t, ζ)

    image

    [ +

    ∂n∂n1

    4G(t, ζ)

    image

    1

     

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

    3G(t, ζ)

    image

    1

     

    k(t)( n2n

    3G(t, ζ)

    image

    ∂τ 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

    image

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

    j=1

    N

    j=1

    N

    image

    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 η. Считаем, что компенсирующие нагрузки в пределах этих элементов

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

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

     

    image

     

    Рис. 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 и ζ равно

    image

    image

    image

    2

     

    B

     

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

     

    • xA

       

      )2 + (yB

       

    • yA

    image

    )2 = |t0| δ,

    2

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

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

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

    x ξ =

    0 − −

     

    (xB xA) t ; y η = 2

    (yB yA) t .

    2 0

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

    image

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

    image

    (

     

    (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. При численном решении каждый контур

image

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 определяются из граничных условий на контурах пластины. Как видно, в пределах точности графиков результаты решений совпадают.

 

image

 

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

 

image

 

a б

 

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

 

Пример 2

 

 

image

 

Рис. 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 · 104 град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

image

∇ ∇ 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 .

 

image

 

a б

 

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

 

Пример 3

 

image

 

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 · 104 град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

  1. Velikanov P.G. Investigation of the isotropic plates bending lying on the complex two-parameter elastic foundation by boundary element method. Izvestiya of Saratov University. Mathematics. Mechanics. Informatics, 2008, vol. 8, issue 1, pp. 36–42. DOI: http://doi.org/10.18500/1816-9791-2008-8-1-36-42. (In Russ.)
  2. Kukanov N.I., Velikanov P.G. Temperature bending of single-connected and multi-connected isotropic plates by the indirect boundary elements method. In: Actual problems of modern science: proceedings of the 2nd International forum. Natural sciences. Parts 1–3. Samara: Izd-vo SGTU, 2006, pp. 178–183. (In Russ.)
  3. Velikanov P.G., Kukanov N.I., Khalitova D.M. Nonlinear deformation of a cylindrical panel of step-variable stiffness on an elastic base by the method of boundary elements. In: Actual problems of continuum mechanics-2020. Izhevsk: OOO "Print 2020, pp. 111–115. Available at: https://www.elibrary.ru/item.asp?id=44838517; https://kpfu.ru/portal/docs/F600474390/apcm_2020.pdf. (In Russ.)
  4. Artyukhin Yu.P., Gribov A.P. Solving of nonlinear deformation problems of plates and shallow shells by the boundary elements method. Kazan: Fen, 2002, 199 p. Available at: https://booksee.org/book/438702. (In Russ.)
  5. Shevchenko V.P. Integral transformations in the theory of plates and shells. Donetsk: Donetskii gosudarstvennyi universitet, 1977, 115 p. (In Russ.)
  6. Korenev B.G. Some problems of the theory of elasticity and thermal conductivity solved in Bessel functions. Moscow: Fizmatgiz, 1960, 458 p. (In Russ.)
  7. Panich O.I. Potentials for polyharmonic equations of the fourth order. Matematicheskii Sbornik. Novaya Seriya, 1960, Vol. 50 (92), Issue 3, pp. 335–354. Available at: http://mi.mathnet.ru/msb4795. (In Russ.)

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

Creative Commons License
This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.

This website uses cookies

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

About Cookies