Численный алгоритм исследования напряженно-деформированного состояния цилиндрических оболочек железнодорожных цистерн Текст научной статьи по специальности «Математика»

УДК 539.3 Вестник СПбГУ. Математика. Механика. Астрономия. 2019. Т. 6 (64). Вып. 2

МБС 74К25

Численный алгоритм исследования напряженно-деформированного состояния цилиндрических оболочек железнодорожных цистерн

П. В. Герасименко, В. А. Ходаковский

Петербургский государственный университет путей сообщения Императора Александра I, Российская Федерация, 190031, Санкт-Петербург, Московский пр., 9

Для цитирования: Герасименко П. В., Ходаковский В. А. Численный алгоритм исследования напряженно-деформированного состояния цилиндрических оболочек железнодорожных цистерн // Вестник Санкт-Петербургского университета. Математика. Механика. Астрономия. 2019. Т. 6(64). Вып. 2. С. 308-317. https://doi.org/10.21638/11701/spbu01.2019.213

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

Ключевые слова: метод сеток, ортогональные развертки, краевая задача, напряженно-деформированное состояние оболочки вращения, метод развертки, векторы ортогона-лизации по Граму — Шмидту.

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

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

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

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

(¡5 Санкт-Петербургский государственный университет, 2019

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

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

s=0 в

B0y (0, в) = 0; Bny (1,в)= bN. (1)

Здесь а, в — координаты срединной поверхности соответственно в меридиональном и окружном направлениях; A(s) (а) —квадратные матрицы 8-го порядка (выражения для их элементов не приводятся из-за громоздкости); yT (а, в) = [Na, Sa, Qa, Ma, u, v, w, — транспонированный вектор столбец усилий Na, Sa, Qa, момента Ma, перемещений u,v,w и угла "ва поворота; f (а, в) — вектор нагрузки; l — длина оболочки; Bk (k = 1,N) —прямоугольные матрицы 4 х 8, характеризующие условия закрепления оболочки на краях и имеющие по восемь ненулевых элементов

= Y(k); &2+s = 1 - Y(k), s = 1, 2, 3, 4; k =1,M,

(k)

где Ys — числа, принимающие значения, равные единице, если s-я компонента вектора y (а, в) задана, и равна нулю, если не задана.

Компоненты вектора нагрузки имеют вид

f 1(а, в) = 2nRpMla; f 3(а, в) = 2пЕ2аржа(1 + sin в); f 2(а,в) = f 4(а,в) = ... = f 8(а,в)=0,

где рж — плотность перевозимой жидкости, рм — плотность металла цилиндрической оболочки, R — радиус цилиндрической оболочки, a — ускорение, возникающее при интенсивном торможении цистерны, l — длина цилиндрической оболочки, вектор (bN)T = [Nнаг 0 0 0 0 0 0 0], а его ненулевая компонента равна Nнаг = 2nR2lpyKa.

Решение поставленной задачи возможно только приближенными методами, из числа которых широкое применение нашли численные методы. Среди численных методов приемлемым методом для рассматриваемых задач является метод сеток. Этот метод достаточно хорошо разработан теоретически. Одним из направлений применения метода является сочетание его с методом прогонки [2]. Вместе с тем, при практической реализации его появлялся ряд проблем, которые создавали немалые трудности. Эти трудности были вызваны потребностью большого времени решения задачи, необходимостью устранять накопление вычислительных погрешностей в силу большого числа вычислительных операций и возможной неустойчивостью избранной разностной схемы. В современных условиях время решения сокращено за счет применения компьютеров. Накопление вычислительных погрешностей, возникающее из-за необходимости «прогонять» большое число векторов размерностью

8М, число которых составляет 4М + 1, удается избежать за счет применения ортогональной прогонки [2].

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

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

2

Уi+1,j =

в=-2

Аг,3+з уг,з + з + ^г,

0, 2, — 1; з = 1, 2,..., М,

(2)

где г, N и з, М — номер и число узлов соответственно в меридиональном и окруж-

ном направлениях, у/. = ,

и>г,з, ] — транспонирован-

ный вектор-столбец компонент напряженно-деформированного состояния в узле (г,3), = [//., /¿2з-,..., /8.] — транспонированный вектор-столбец произвольных силовых, моментных и температурных нагрузок, действующих на оболочку, Аг, (в = -2, —1,0,1, 2) — квадратные матрицы восьмого порядка, зависящие в общем случае от номера узла сетки как в меридиональном, так и в окружном направлениях, N¿,3, ¿г., Qilj, Мг. —усилия и момент в узле (г,з), иг., , и>г., — перемещения и угол поворота в узле (г, з).

Элементы матриц (в = —2, —1,0,1, 2) не приводятся ввиду их громозд-

кости.

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

Во,. уо,з = 0 Воуо,з = 0, В^,. ум,з = ъ^-, Вдт= ь^,..

(3)

Здесь Во,., — прямоугольные матрицы порядка 4x8, характеризующие условия закрепления оболочки на краях; 0 — нулевой вектор-столбец;

В0,.

10000000 01000000 00100000 00010000

10000000 01000000 00100000 00010000

Во2з

Во ,з

0 0 0 0 1 0 0 0 00000100 00000010 00000001

10000000 00000100 00000010 00000001

(Ъ^з)Т = [Nнаг 0 0 0 0 0 0 0], з = 1, 2,

N.

Значения компонент вектора fi}j определяются по следующим зависимостям: flj = 2nRpMla,i; fj = 2пВ2а.ржа.(1 + ът вj); fj = fj = ... = ff¿ =0. Систему уравнений (2), (3) можно записать в более компактной форме:

yi+i = Aiyi + fi, г = 0, 1, 2,..., N - 1, (4)

где Ai —в основном блочная пятидиагональная матрица в каждом г-м узле:

Ai =

Ai,i Ai, 2 Ai,3 0 0 0 0 0 0 0 Ai ,M-i Ai,M

Ai, i Ai, 2 Ai,3 Ai, 4 0 0 1 0 0 0 0 0 Ai,M

Ai,i Ai, 2 Ai,3 Ai, 4 Ai, 5 0 1 0 0 0 0 0 0

Ai,2 Ai,3 Ai,4 Ai, 5 Ai,e 1 0 0 0 0 0 0

0 0 0 0 0 0 . Ai MS Ai ,M-4 Ai ,M-3 Ai , M-2 Ai ,M-i 0

0 0 0 0 0 0 0 Ai ,M-4 Ai ,M-3 Ai , M-2 Ai ,M-i Ai,M

Ai, i 0 0 0 0 0 1 0 0 Ai ,M-3 Ai , M-2 Ai ,M-i Ai,M

Ai, i Ai,2 0 0 0 0 1 0 0 0 Ai , M-2 Ai ,M-i Ai,M

УТ = [Ум> УТ21 уТ,Ы ] и ^ = , ^2 1 ] - вектор-столбцы.

Матрица А —блочная, почти пятидиагональная. Пятидиагональность нарушают элементы первых двух и последних двух строк. Учитывая осевую симметрию оболочки, в матрице появляются элементы: А^м-1, А,м в первой строке и А.1,м во второй строке. В последних двух строках введены элементы: А^1 в предпоследней строке и А^1, А^2 в последней строке.

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

Во,¿Усу = 0, В2^ук = 0, 2 = М/2 + 1, М; (5)

B

0,j

" Bl,1 o ;

0 ,2 i0, B

0 0 1

0 o :

0 0 1

0 0

RÍ---

B0,M/2

О

~ "О" ~

0

I

0

_ _0 _ _ 1 _ о о Т~о

R2 1

_°j.M/2+i j_

0

2

B0,M

BN ,j yN,j

bi , bN,j,

j = 1, M/2;

BN2 ,jyN,j

b2 , bN,j,

j = 1, M/2;

(6)

BN,j =

BNi ,i 0 1 0 0 0

0 Bi TV,2 1 0 ■"Rl 1 TV,M/2 0 0

и 0 0 0

0 0 : o B2 N ,M [2 + \. 0

0 0 ! 0 0

Для решения системы (4)—(6) в работе [4] предложен метод ортогональной прогонки, который позволяет избежать накопления вычислительных погрешностей. Предложенный метод базируется на методе С.К.Годунова [5]. Он использует ортогональную прогонку для решения численным методом обыкновенных дифференциальных уравнений. Метод широко используется при решении многих задач теории

0

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

уг

М 4

*(0)+ЕЕ в

3=1 к=1

(3,^(3,^)

0, 1, ..., N,

(7)

где z(з'fc) (к = 1, 2, 3, 4; з = 1, 2,..., М)—система ортонормированных векторов, z(0) — вектор, ортогональный к системе векторов z(з'fc) (к = 1, 2, 3, 4; з = 1, 2,..., М), вг(3'к) (к = 1, 2, 3,4; з = 1, 2,..., М) — произвольные постоянные коэффициенты в узле г,

( уОд \

0

0 0

„(3,*0

< ¿8 ^

(к) у0,.

0

0

(М,к)

0 0

0 0

\ у(к) /

у0,М

(8)

Векторы z03'fc) (к = 1, 2, 3,4; з = 1, 2,..., М), задаются таким образом, чтобы они были ортонормированы, и каждый из них удовлетворял граничному условию (5). В качестве начальных векторов можно выбрать

(0)

5о,з

/1 — 71 0)\ 0 0 0

71 0 0 0

( )

,(2) 5о,з

/

0

1

72

( )

\

72 0 0

( )

,(з) =

50,3 =

1

7з

( )

Остальные векторы z,

V

(3,к)

0 0 0

7з 1 )

(0)

находятся по следующим формулам:

у

(3,к)

Аг-1Аг_ 1 ,

,(4) = 50,3 =

1

74

(1)

0 0 0

741)

(г = 1, 2,. .., N; з

/

1, 2,..

(0)

,М; к = 1, 2, 3,4)

Аг-1 zf\ +

0

0

0

0

0

0

= Ы3 — ЕЕ ^Ч-Ч

V р=1в=1

,(р,3)

,,(р,3) у(3,к) ^ , уг

, в = 1, 2, ...,п; р = 1, 2, ...,з,

м(г-3) =

у(3,к) у(3,к)

уг , уг

—

Р=1 в=1

,(р,к)

, к = 1, 2, 3,4; з = 1, 2,..., М,

М4

,(°)_„(°) ^\^,,(3,к)„(3,к)

у(0) —££ , 3 = у(0), *) , г = 0,1,..., N

3=1 к=1

где п =

к — 1, з = 1,

к, з > 1.

При решении системы разностных уравнений (4)—(6) методом, описанном в пособии [3], на каждом шаге ортогональной прогонки необходимо сохранять в памяти

компьютера систему векторов z ортонормирования

(0)

zг(3'fc) (к

1, 2, 3, 4; з = 1, 2,..., М) и матрицу

«г =

«1,1 0 0

0 0Т

«1,2 «1,3 «2,2 «2,3 0 «3,3

0 0Т

0 0Т

«1,М «2,М «3,М

«

М,М

,(1) ,(2) ,(3)

(9)

Здесь 0 — нулевая матрица четвертого порядка;

«г,3 =

«р

,(г,3) , ,(г,3) , ,(г,3) , ,(г,3)

1,1

0 0 0

1,1

2,1

[^0 ]

1,2

1,3

,(г,3) , ,(г,3)

2,2 0 0

2,3

,(г,3) , ,(г,3)

1,4

(г,3)

2,4

3,3 0

3,4

(г,3)

4,4

1,2

1,3

1,4

2,2

2,3

2,4

3,2

3,3

3,4

4,2

4,3

(3) ,,(3) ,,(3)

0,2 , 0,3 , 0,4

с^А

4,4 з

,(Р,3) . ,(Р,3) . ,(Р,3) . ,(Р,3)

,(Р,3) . ,(Р,3) . ,(Р,3) . ,(Р,3)

,(Р,3) . ,(Р,3) . ,(Р,3)

,(Р,3) , ,(Р,3) , ,(Р,3) , ,(Р,3) 3, 1 , ,(р,3) ^4,1

- Г, ,(3)

- И0,1,

з = 1, 2,

1, 2,

2, 3,

М.

М;

М;

1, 2,

з — 1;

После вычисления системы векторов zг

э(3,к)

(°) „(3,к) /,. _

(к = 1, 2, 3, 4; з = 1, 2,..., М) произвольные постоянные вги'К'' (к = 1, 2, 3,..., 4М) в узле г = N определяются из следующей системы алгебраических уравнений:

М4

\ \ в(3,к) В „(3,^) = В „(°)

3=1 к=1

2

к

0

0

1

з

р

Рис. 1. Схема нагружения цилиндрической оболочки.

В остальных узлах сетки г = N — 1, N — 2,..., 0 эти постоянные находятся по следующей рекуррентной формуле:

= [а]-1, г = N — 1^ — '2,..., 0,

где [ДР]Т = |в1(1'1),в1(1'2),в1(1'3),в1(1'4),в1(2'1),..., 1].

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

В качестве примера в работе выполнен расчет напряженно-деформированного состояния круговой цилиндрической оболочки с постоянной толщиной (рис. 1).

Рис. 2. Зависимость относительного прогиба оболочки от окружной координаты.

Рис. 3. Изменение относительного прогиба оболочки по меридиональной координате.

Оболочка была жестко защемленная на всей окружности, где и = v = w = v = 0 в сечении а = 0, и шарнирно опертая, где Na = Ma = v = w = 0 на части окружности ß = ±0.8п (144° ^ ß ^ 216°) в сечении а = l. На оставшейся части окружности она имела свободное крепление (Na = S = Q = Ma = 0).

Оболочка находилась под действием распределенных поверхностных нормальных qn = qo cos ß и касательных qT = qo sin ß сил, которые возникают вследствие инерционной нагрузки. Необходимо отметить, что такая расчетная схема может выступать при оценивании напряженно-деформированного состояния оболочки котла железнодорожной цистерны при условии движения состава по железнодорожному пути с большими зазорами на стыке рельсов.

Расчеты выполнялись для модельной оболочки при следующих параметрах: l/R = 1, R/h = 50. На рис.2 показано изменение относительного прогиба w = w/(qo/E) в зависимости от окружной координаты ß в сечениях а = 1/2 и а = l.

На рис. 3 представлена зависимость относительного прогиба от меридиональной координаты а в сечениях ß = 0 и ß = п.

При решении данной задачи учитывалась симметрия напряженно-деформированного состояния относительно плоскости ß = 0, п, что позволило использовать равномерную прямоугольную разностную сетку, у которой N = 20, а M = 10. Для аппроксимации частных производных применялась двумерная разностная схема второго порядка точности в обоих направлениях.

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

Литература

1. Григоренко Я. М., Василенко А Т., Беспалова Е. И. и др. Численное решение задач статики ортотропных оболочек с переменными параметрами. Киев: Наукова думка, 1975.

2. Герасименко П. В. Численный алгоритм исследования устойчивости дискретно подкрепленных шпангоутами ортотропных оболочек вращения // Прикладная механика. 1975. Т. 11, №11. С. 50-57.

3. Годунов С. К., Рябенький В. С. Разностные схемы (введение в теорию): учебное пособие. М.: Наука, 1977.

4. Герасименко П. В., Павлов Г. А. Об одном численном методе решения задач статики оболочек вращения // Прикладная механика. 1980. Т. 16, №5. С. 62-66.

5. Годунов С. К. Метод ортогональной прогонки для решения систем разностных уравнений // Журн. вычисл. математики и мат. физики. 1962. Т. 2, №6. С. 972-980.

6. Виноградов А. Ю. Численные методы решения жестких и нежестких краевых задач. М.: National Research, 2017.

Статья поступила в редакцию 14 ноября 2018 г.;

после доработки 18 декабря 2018 г.; рекомендована в печать 20 декабря 2018 г.

Контактная информация:

Герасименко Петр Васильевич —д-р техн. наук, проф.; pv39@mail.ru Ходаковский Валентин Аветикович —д-р техн. наук, проф.; hva1104@mail.ru

Numerical algorithm of investigation of stress-strain state of cylindrical shells of railway tank

P. V. Gerasimenko, V. A. Khodakovskiy

Emperor Alexander I St. Petersburg State Transport University, Moskovsky pr., 9, St. Petersburg, 190031, Russian Federation

For citation: Gerasimenko P. V., Khodakovskiy V. A. Numerical algorithm of investigation of stress-strain state of cylindrical shells of railway tank. Vestnik of Saint Petersburg University. Mathematics. Mechanics. Astronomy, 2019, vol. 6(64), issue 2, pp. 308-317. https://doi.org/10.21638/11701/spbu01.2019.213 (In Russian)

On the basis of synthesis of methods of grids and orthogonal sweeps of grids system of algebraic equations of big dimension the algorithm of the decision of a boundary value problem in partial derivatives describing the stress-strain state of a shell of rotation of a copper of the railway tank is offered. The shell is cylindrical, loaded arbitrarily by inertial forces and pressure, has combined fastening conditions at the ends. The method of grids according to the explicit scheme allowed the system of equations of the theory of shells reduced to eight differential equations of the first order with respect to the meridional coordinate and the fourth order with respect to the circumferential, to lead to a system of algebraic equations having a five-diagonal matrix, non-zero elements of which are eighth order matrices. To solve a system of algebraic equations with a rare matrix of non-zero elements, the method of run is applied, in which to eliminate the accumulation of computational errors, the or-thogonalization of "run" vectors by Gram — Schmidt is applied, which allows to exclude the formation of a singular matrix from the "run" vectors when calculating the coefficients of solving the boundary value problem. An example is considered in which the stress-strain state of the shell of the boiler of a railway tank tests the inertial forces variable along the meridional and circumferential coordinates and has a rigid fastening at each end on the part of the circle and a free state on the remaining part.

Keywords: method of grids, orthogonal sweeps, boundary value problem, the stress-strain

state of the shell of revolution, sweep method, orthogonalization vectors by Gram —

Schmidt.

References

1. Grigorenko Y. M., Vasilenko A.T., Bespalova E.I. et al., Numerical solution of statics problem of orthotropic shells with variable 'parameters (Naukova Dumka Publ., Kiev, 1975). (In Russian)

2. Gerasimenko P. V., "Numerical algorithm for investigation of the stability of discretely reinforced frames of orthotropic shells of revolution", Soviet Applied Mechanics 11 (11), 1175—1180 (1975). https: //doi.org/10.1007/BF00883277

3. Godunov S. K., Ryaben'kii V. S., Difference schemes (introduction to theory) (Nauka Publ., Moscow, 1977). (In Russian)

4. Gerasimenko P. V., Pavlov G.A., "Numerical method of solving statics problems for shells of revolution", Soviet Applied Mechanics 16 (5), 409-413 (1980). https://doi.org/10.1007/BF00884172

5. Godunov S. K., "A method of orthogonal successive substitution for the solution of systems of difference equations", USSR Computational Mathematics and Mathematical Physics 2 (6), 1151-1165 (1963). https://doi.org/10.1016/0041-5553(63)90339-2

6. Vinogradov A. Y., Numerical methods of .solving stiff and non-rigid boundary value problems (National Research Publ, Moscow, 2017). (In Russian)

Received: November 14, 2018 Revised: December 18, 2018 Accepted: December 20, 2018

A u t h o r's i n f o r m a t i o n:

Petr V. Gerasimenko — pv39@mail.ru Valentin A. Khodakovskiy — hva1104@mail.ru