1
РЕШЕНИЕ ЖЕСТКИХ КРАЕВЫХ ЗАДАЧ СТРОИТЕЛЬНОЙ МЕХАНИКИ
(РАСЧЕТ ОБОЛОЧЕК СОСТАВНЫХ И СО ШПАНГОУТАМИ) МЕТОДОМ
ВИНОГРАДОВЫХ (БЕЗ ОРТОНОРМИРОВАНИЯ)
Алексей Юрьевич Виноградов, кандидат физико-математических наук, freelance
Юрий Иванович Виноградов, доктор физико-математических наук,
профессор МГТУ им. Н.Э.Баумана, Москва, AlexeiVinogradov@yandex.ru
Предлагается простейший метод решения жёстких краевых задач строительной механики для
расчета тонкостенных составных оболочек и оболочек со шпангоутами. Не требуются процедуры
ортонормирования, что достигается за счёт разделения интервала интегрирования на сопрягаемые
участки.
Ключевые слова: жесткие краевые задачи, оболочки, шпангоуты, ортонормирование.
SOLUTION OF STIFF BOUNDARY VALUE PROBLEMS OF BUILDIND MECHANICS
(CALCULATION OF SHELLS COMPOUND AND WITH FRAMES) BY VINOGRADOVS’
METHOD (WITHOUT ORTHONORMALIZATION)
Alexei Yurievich Vinogradov, candidate of physical and mathematical sciences, freelance
Yuri Ivanovich Vinogradov, doctor of physical and mathematical sciences,
professor of MSTU n.a. Bauman, Moscow, AlexeiVinogradov@yandex.ru
The simplest method of solution of the stiff boundary value problems of building mechanics is offered for
calculation of compound shells and shells with frames. The method does not need orthonormalization
because of dividing an integration interval into matching sectors.
Keywords: stiff boundary value problems, shells, frames, orthonormalization.
1. Введение.
Введение
дается
на
примере
системы
дифференциальных
уравнений
цилиндрической оболочки ракеты – системы обыкновенных дифференциальных
уравнений 8-го порядка (после разделения частных производных методом Фурье).
Система линейных обыкновенных дифференциальных уравнений имеет вид:
Y ( x) AY ( x) F ( x) ,
2
где Y (x) – искомая вектор-функция задачи размерности 8х1, Y (x) – производная
искомой вектор-функции размерности 8х1, A – квадратная матрица коэффициентов
дифференциального уравнения размерности 8х8, F ( x) – вектор-функция внешнего
воздействия на систему размерности 8х1.
Краевые условия имеют вид: UY (0) u,VY (1) v ,
Y (0) – значение искомой вектор-функции на левом крае х=0 размерности 8х1, U –
где
прямоугольная горизонтальная матрица коэффициентов краевых условий левого края
размерности 4х8, u – вектор внешних воздействий на левый край размерности 4х1,
Y (1) – значение искомой вектор-функции на правом крае х=1 размерности 8х1, V –
прямоугольная горизонтальная матрица коэффициентов краевых условий правого края
размерности 4х8, v – вектор внешних воздействий на правый край размерности 4х1.
В случае, когда система дифференциальных уравнений имеет матрицу с
постоянными коэффициентами A =const, решение задачи Коши имеет вид [1]:
x
Y ( x) e A( x x 0 )Y ( x0 ) e Ax e At F (t )dt ,
x0
где e A( x x 0) E A( x x0 ) A 2 ( x x0 ) 2 / 2! A 3 ( x x0 ) 3 / 3!... , где E - единичная матрица.
Матричная экспонента ещё может называться матрицей Коши или матрициантом и
может обозначаться в виде K ( x x0 ) K ( x x0 ) e A( x x 0) . Тогда решение задачи Коши
может быть записано в виде Y ( x) K ( x x0 )Y ( x 0 ) Y ( x x0 ) ,
x
где Y ( x x 0 ) e
Ax
e
At
F (t )dt это вектор частного решения неоднородной системы
x0
обыкновенных дифференциальных уравнений.
Вместо формулы для вычисления вектора частного решения неоднородной
x
системы
дифференциальных
уравнений
в
виде
[1]
Y ( x x 0 ) e Ax e At F (t )dt
x0
3
предлагается использовать следующую формулу для каждого отдельного участка
xj
интервала интегрирования Y ( x j xi ) Y ( x j xi ) K ( x j xi ) K ( xi t ) F (t )dt .
xi
Вычисление
этого
вектора
частного
решения
неоднородной
системы
дифференциальных уравнений производиться при помощи представления матрицы Коши
под знаком интеграла в виде ряда и интегрирования этого ряда поэлементно. Это
справедливо для случая системы дифференциальных уравнений с постоянной матрицей
коэффициентов A =const. Для случая дифференциальных уравнений с переменными
коэффициентами
в
приведенной
использоваться
осредненная
выше
матрица
формуле
для
Ai A( xi )
каждого
участка
коэффициентов
может
системы
дифференциальных уравнений.
Вектор F (t ) может рассматриваться на участке ( x j xi ) приближенно в виде
постоянной величины F ( xi ) constant , что позволяет вынести его из-под знака интеграла,
что приводит к совсем простому ряду для вычислений на рассматриваемом участке.
Рассмотрим вариант, когда шаги интервала интегрирования выбираются достаточно
малыми, что позволяет рассматривать вектор F (t ) на участке ( x j xi ) приближенно в
виде постоянной величины F ( xi ) constant , что позволяет вынести этот вектор из-под
знаков интегралов. Тогда получаем ряд для вычисления вектора частного решения
неоднородной системы дифференциальных уравнений на малом участке ( x j xi ) :
Y ( x j xi ) K ( x j xi ) ( E A( xi x j ) / 2! A 2 ( xi x j ) 2 / 3!...) ( x j xi ) F ( xi ).
Если участок ( x j xi ) не мал, то его можно поделить на подучастки и тогда
можно предложить следующие рекуррентные (итерационные) формулы для вычисления
частного вектора: имеем Y ( x) K ( x x0 )Y ( x 0 ) Y ( x x0 ) , также имеем формулу для
4
xj
отдельного подучастка
Y ( x j xi ) Y ( x j xi ) K ( x j xi ) K ( xi t ) F (t )dt . Тогда:
xi
Y ( x1 ) K ( x1 x0 )Y ( x0 ) Y ( x1 x0 ) , Y ( x 2 ) K ( x 2 x1 )Y ( x1 ) Y ( x 2 x1 ) .
Подставим Y ( x1 ) в Y ( x 2 ) и получим:
Y ( x 2 ) K ( x 2 x1 ) K ( x1 x 0 )Y ( x0 ) K ( x 2 x1 )Y ( x1 x0 ) Y ( x 2 x1 ) .
Сравним это выражение с формулой Y ( x 2 ) K ( x 2 x0 )Y ( x0 ) Y ( x 2 x0 )
и получим, очевидно, что K ( x 2 x 0 ) K ( x 2 x1 ) K ( x1 x0 ) и для частного вектора
получаем формулу Y ( x 2 x0 ) K ( x 2 x1 )Y ( x1 x 0 ) Y ( x 2 x1 ) . То есть вектора
подучастков Y ( x1 x0 ), Y ( x 2 x1 ) не просто складываются друг с другом, а с
участием матрицы Коши подучастка.
Аналогично можем записать Y ( x3 ) K ( x3 x 2 )Y ( x 2 ) Y ( x3 x 2 ) и подставить
сюда формулу для Y ( x 2 ) и т.д.
Из теории матриц [1] известно свойство перемножаемости матричных экспонент
(матриц Коши): K ( xi x0 ) K ( xi xi 1 ) ... K ( x 2 x1 ) K ( x1 x0 ) .
В случае, когда система дифференциальных уравнений имеет матрицу с
переменными коэффициентами A A( x) , решение задачи Коши можно искать (как это
известно из теории матриц) при помощи свойства перемножаемости матриц Коши. То
есть интервал интегрирования разбивается на малые участки и на малых участках
матрицы Коши приближенно вычисляются по формуле для постоянной матрицы в
экспоненте. А затем матрицы Коши, вычисленные на малых участках, перемножаются:
K ( xi x0 ) K ( xi xi 1 ) ... K ( x 2 x1 ) K ( x1 x0 ) ,
где матрицы Коши приближенно вычисляются по формуле:
K ( xi1 xi ) e A( xi ) xi exp( A( xi ) xi ) , где xi xi 1 xi .
5
В этом случае использования перемножения матричных экспонент, то есть в случае
кусочно-константной
аппроксимации
матрицы
A( x)
системы
обыкновенных
дифференциальных уравнений (ОДУ), оценка точности в литературных источниках не
встречается.
На
практике
можно
удостоверяться
в
точности
расчетов
путем
сравнительного расчета при разных шагах, на которые делится интервал интегрирования,
или можно сразу разделять интервал интегрирования на очевидно малые шаги
интегрирования.
В тоже время можно производить вычисления и с заранее известной точностью.
Для этого следует вычислять матрицы Коши не как матричные экспоненты от
осредненных аргументов, а можно использовать для вычисления векторов, входящих в
матрицы Коши, методы типа методов Рунге-Кутты. Для вычисления матриц Коши
методами типа Рунге-Кутты используется стартовая (начальная) единичная матрица. А
для вычисления вектора частного решения неоднородной системы ОДУ берется
начальный нулевой вектор. В случае применения методов типа Рунге-Кутты оценки
погрешностей хорошо известны.
2. Метод решения жестких краевых задач без ортонормирования – метод сопряжения
участков, выражаемых матричными экспонентами (матрицами Коши).
Идея преодоления трудностей неустойчивого счета путем разделения интервала
интегрирования на сопрягаемые участки принадлежит д.ф.-м.н. Юрию Ивановичу
Виноградову [2]. А выражение идеи разделения и сопряжения через формулы теории
матриц, то есть через матричные экспоненты (матрицы Коши) принадлежит к.ф.-м.н.
Алексею Юрьевичу Виноградову.
Разделим интервал интегрирования краевой задачи, например, на 3 участка. Будем
иметь точки (узлы), включая края x0 , x1 , x 2 , x3 .
Имеем краевые условия в виде UY ( x 0 ) u, VY ( x3 ) v.
Можем записать матричные уравнения сопряжения участков (при переходах от
левых точек к правым точкам участков интервала интегрирования, то есть при переходах
типа от x0 к x1 , хотя вполне можно было бы записать и уравнения при переходах в
обратном направлении, то есть при переходах типа от x1 к x0 ):
6
Y ( x1 ) K ( x1 x0 )Y ( x0 ) Y ( x1 x0 ) ,
Y ( x 2 ) K ( x 2 x1 )Y ( x1 ) Y ( x 2 x1 ) ,
Y ( x3 ) K ( x3 x 2 )Y ( x 2 ) Y ( x3 x 2 ) .
Это мы можем переписать в виде, более удобном для нас далее:
EY ( x1 ) K ( x1 x 0 )Y ( x 0 ) Y ( x1 x 0 ) ,
EY ( x 2 ) K ( x 2 x1 )Y ( x1 ) Y ( x 2 x1 ) ,
(1)
EY ( x3 ) K ( x3 x 2 )Y ( x 2 ) Y ( x3 x 2 ) .
где E - единичная матрица.
В итоге получаем систему линейных алгебраических уравнений:
U
K ( x1 x0 )
0
0
0
0
E
u
0
0
Y ( x0 )
Y ( x1 x0 )
0
0
Y ( x1 )
Y ( x 2 x1 ) . (2)
K ( x 2 x1 )
E
0
Y ( x2 )
Y ( x3 x 2 )
K ( x3 x 2 ) E
0
Y ( x3 )
V
v
0
0
Эта система решается методом Гаусса с выделением главного элемента.
Оказывается, что применять ортонормирование в рамках предлагаемого метода не
нужно.
Это объясняется тем, что на каждом отдельном участке интервала интегрирования
вычисления матриц Коши (матричных экспонент) ведется каждый раз отдельно
(независимо от других участков) со стартом численного интегрирования (либо в виде
матричной экспоненты либо методами типа Рунге-Кутты) от начальной единичной
матрицы, то есть от ортонормированной (единичной) матрицы, что происходит
естественно без каких-либо дополнительных операций типа операций ортонормирования
векторов, которые применяются в методе С.К.Годунова, что существенно упрощает
программирование предлагаемого метода по сравнению с методом С.К.Годунова.
7
В точках вблизи узлов решение находится путем решения соответствующих задач
Коши с началом в i-ом узле Y ( x) K ( x xi )Y ( xi ) Y ( x xi ) .
3. Расчет составных оболочек вращения.
Пусть имеем, например, 3 участка интервала интегрирования составной оболочки
вращения, где каждый участок может выражаться своими дифференциальными
уравнениями и физические параметры могут выражаться по-разному – разными
формулами на разных участках:
Рис. 1
В общем случае (на примере участка 12) физические параметры участка (вектор
P12 ( x) )
выражаются
через
искомые
параметры
системы
обыкновенных
дифференциальных уравнений этого участка (через вектор Y12 ( x) ) следующим образом
P12 ( x) M 12Y12 ( x) , где матрица M 12 - квадратная невырожденная.
При переходе точки сопряжения можем записать в общем виде (но на примере
точки сопряжения x1 : P01 ( x1 ) P0112 L0112 P12 ( x1 ) , где P0112 - дискретное приращение
физических параметров (сил, моментов) при переходе с участка «01» на участок «12», а
матрица L0112 квадратная невырожденная диагональная и состоит из единиц и минус
единиц на главной диагонали для установления правильного соответствия принятых
положительных направлений сил, моментов, перемещений и углов при переходе с участка
«01» на участок «12», которые могут быть разными (в разных дифференциальных
уравнениях разных сопрягаемых участков) – в уравнениях слева от точки сопряжения и в
уравнениях справа от точки сопряжения.
Два последних уравнения образуют M 01Y01 ( x1 ) P0112 L0112 M 12Y12 ( x1 ) .
8
В точке x 2 аналогично получим M 12Y12 ( x 2 ) P12 23 L12 23 M 23Y23 ( x 2 ) .
Если бы оболочка состояла бы из одинаковых участков, то мы могли бы записать в
объединенном матричном виде систему линейных алгебраических уравнений в форме (2).
Но в нашем случае оболочка состоит из 3 участков, где средний участок можно считать,
например, шпангоутом, выражаемым через свои дифференциальные уравнения. Тогда
вместо векторов Y ( x0 ) , Y ( x1 ) , Y ( x 2 ) , Y ( x3 ) мы должны рассмотреть вектора:
Y01 ( x0 ), Y01 ( x1 ), Y12 ( x1 ), Y12 ( x 2 ), Y23 ( x 2 ), Y23 ( x3 ) (смотри рис.1).
Тогда матричные уравнения UY ( x0 ) u,VY ( x3 ) v и (1) примут вид:
UY01 ( x0 ) u, VY23 ( x3 ) v ,
EY01 ( x1 ) K 01 ( x1 x 0 )Y01 ( x 0 ) Y01* ( x1 x 0 ) ,
M 01Y01 ( x1 ) P0112 L0112 M 12Y12 ( x1 ) ,
EY12 ( x 2 ) K 12 ( x 2 x1 )Y12 ( x1 ) Y12* ( x 2 x1 ) ,
M 12Y12 ( x 2 ) P12 23 L12 23 M 23Y23 ( x 2 ) ,
EY23 ( x3 ) K 23 ( x3 x 2 )Y23 ( x 2 ) Y23* ( x3 x 2 ) .
Тогда получаем итоговую систему линейных алгебраических уравнений:
U
K 01 ( x1 x0 )
0
0
0
0
0
0
E
M 01
0
0
0
0
0
0
0
0
L0112 M 12
0
K 12 ( x 2 x1 ) E
M 12
0
0
0
0
0
0
0
0
0
u
0
Y01 ( x0 )
*
Y01 ( x1 x0 )
0
Y01 ( x1 )
P0112
0
Y12 ( x1 )
*
Y12 ( x 2 x1 )
0
Y12 ( x 2 )
L12 23 M 23
P12 23
0
Y23 ( x 2 )
*
Y23 ( x3 x 2 )
K 23 ( x3 x 2 ) E
Y23 ( x3 )
V
v
0
4. Расчет оболочек со шпангоутами, когда шпангоуты выражаются алгебраическими
уравнениями.
Рассмотрим случай, когда шпангоут (в точке
x1 ) выражается не через
дифференциальные уравнения, а через алгебраические уравнения.
9
Рис. 2
Выше мы записывали, что P 01 ( x 1 ) P 01 12 L 01 12 P12 ( x 1 ) .
Можем представить вектор P01 ( x1 ) силовых факторов и перемещений в виде:
P01 ( x1 )
R01 ( x1 )
S 01 ( x1 )
,
где R01 ( x1 ) - вектор перемещений, S 01 ( x1 ) - вектор сил и моментов.
Алгебраическое уравнение для шпангоута имеет вид GR S , где G – матрица
жесткости шпангоута, R – вектор перемещений шпангоута, S – вектор силовых
факторов, которые действуют на шпангоут.
В точке шпангоута имеем R 0, S GR , то есть нет разрыва в перемещениях
R 0 , (где 0 – нулевой вектор), но есть результирующий вектор силовых факторов
S GR , который складывается из сил и моментов слева плюс сил и моментов справа от
точки шпангоута. Тогда получаем:
P01 ( x1 )
0
R
L0112 P12 ( x1 ) , P01 ( x1 )
L0112 P12 ( x1 ) .
S
GR
R01 ( x1 )
M 11p
p
Запишем P01 ( x1 )
M 01Y01 ( x1 ) M Y01 ( x1 )
S 01 ( x1 )
M 21p
M 12p
Y01 ( x1 ) , где для
M 22p
удобства было введено переобозначение M 01 M p . Тогда можем записать:
R01 ( x1 ) M 11p
g*
0
GR01 ( x1 )
M 12p Y01 ( x1 ) и введем в рассмотрение вектор g * :
0
G M 11p
M 12p Y01 ( x1 )
0...0
G M 11p
M 12p
Y01 ( x1 )
Запишем матричные уравнения для этого случая (рис.2):
UY01 ( x0 ) u,VY12 ( x 2 ) v ,
0
G M 11p
M 12p
Y01 ( x1 )
10
EY01 ( x1 ) K 01 ( x1 x0 )Y01 ( x 0 ) Y01* ( x1 x 0 ) ,
M 01Y01 ( x1 ) g * L0112 M 12Y12 ( x1 ) ,
EY12 ( x 2 ) K 12 ( x 2 x1 )Y12 ( x1 ) Y12* ( x 2 x1 ) .
Распишем вектор g * , тогда: ( M 01
0
GM
p
11
M 12p
) Y01 ( x1 ) L0112 M 12Y12 ( x1 ) .
Для обеспечения негромоздкости введем обозначение ( M 01
0
GM
p
11
M 12p
) M*
и получим M *Y01 ( x1 ) L0112 M 12Y12 ( x1 ) .
Таким образом, получаем итоговую систему линейных алгебраических уравнений:
U
0
K 01 ( x1 x0 ) E
0
M*
0
0
0
0
0
0
u
0
Y01 ( x 0 )
*
Y01 ( x1 x0 )
0
Y01 ( x1 )
.
0
L0112 M 12
0
Y12 ( x1 )
*
K 12 ( x 2 x1 ) E
Y12 ( x 2 x1 )
Y12 ( x 2 )
0
V
v
Если к шпангоуту приложено внешнее силовое-моментное воздействие g p , то
g*
0
0
0
следует переписать в виде g *
, тогда:
p
GR01 ( x1 ) g p
GR
GR g
g*
GM
p
11
0
0...0
0
.
p
p
p Y01 ( x1 )
M Y01 ( x1 ) g
G M 11 M 12
gp
p
12
Тогда получаем M *Y01 ( x1 ) L0112 M 12Y12 ( x1 )
0
.
gp
Итоговая система линейных алгебраических уравнений примет вид:
U
K 01 ( x1 x0 )
0
E
0
0
M*
0
0
0
u
Y01 ( x0 )
Y ( x1 x0 )
0
0
0
Y01 ( x1 )
p
. (3)
L0112 M 12
0
g
Y12 ( x1 )
K 12 ( x 2 x1 ) E
Y12* ( x 2 x1 )
Y12 ( x 2 )
V
0
v
0
0
*
01
11
5. Расчет оболочек составных и со шпангоутами, когда их уравнения выражаются не
через абстрактные вектора, а через вектора, состоящие из конкретных физических
параметров.
Рассмотрим случай, когда части оболочечной конструкции и шпангоут выражаются
через вектора состояния (типа Y12 ( x) ), которые (в частном случае) совпадают с векторами
физических параметров (типа P12 ( x) - перемещения, угол, силы, момент). Тогда матрицы
типа M 12 будут единичными: M 12 E . И пусть положительные направления физических
параметров одинаковы для всех частей оболочки и шпангоута ( L0112 E ).
Тогда будем иметь уравнения P12 ( x) M 12Y12 ( x) , P01 ( x1 ) P0112 L0112 P12 ( x1 ) ,
M 01Y01 ( x1 ) P0112 L0112 M 12Y12 ( x1 ) в виде P12 ( x) EY12 ( x) , P01 ( x1 ) P0112 EP12 ( x1 ) ,
EY01 ( x1 ) P0112 EY12 ( x1 ) , где E – единичная матрица.
Тогда: EY01 ( x1 ) g * EY12 ( x1 ) , EY01 ( x1 )
R01 ( x1 )
E
EY01 ( x1 ) M pY01 ( x1 )
S 01 ( x1 )
0
P01 ( x1 )
g*
0
GR01 ( x1 )
0
EY12 ( x1 ) ,
GR01 ( x1 ) g p
0
Y01 ( x1 ) , R01 ( x1 ) E 0 Y01 ( x1 ) ,
E
0
0
0
0
0
,
Y
(
x
)
01
1
G E 0 Y01 ( x1 )
GE 0
gp
gp
gp
M *Y01 ( x1 ) EY12 ( x1 )
0
0
, где ( E
) M*
p
GE 0
g
Итоговая система линейных алгебраических уравнений (3) примет вид:
U
K 01 ( x1 x0 )
0
E
0
0
0
M*
0
0
u
0
0
*
Y01 ( x 0 )
Y01 ( x1 x0 )
0
0
0
Y01 ( x1 )
p
,
E
0
g
Y12 ( x1 )
K 12 ( x 2 x1 ) E
Y12* ( x 2 x1 )
Y12 ( x 2 )
V
0
v
12
где M * ( E
8 x8
0
4 x8
G E
4x4 4x4
) (E
0
8 x8
4x4
M *Y01 ( x1 ) EY12 ( x1 )
0
4 x8
G
4x4
0
)
4x4
4x4
4x4
4x4
G
E
0
4x4
принимает
вид
G
gp
4x4
E
4x4
G
4x4
E
0 R01 ( x1 )
4x4
E S 01 ( x1 )
0
4x4
4x4
4x4
0
E
4x4
E
0
4x4
E
. Это означает, что уравнение
Y01 ( x1 ) EY12 ( x1 )
4x4
0
или
gp
0 R12 ( x1 )
0
p .
E S12 ( x1 )
g
4x4
4x4
Это означает, что R01 ( x1 ) R12 ( x1 ) 0 , (нет скачка в перемещениях и угле) и
GR01 ( x1 ) S 01 ( x1 ) S12 ( x1 ) g p
(равновесие шпангоута), то есть
R01 ( x1 ) R12 ( x1 )
(перемещения и угол: нет разрыва) и S 01 ( x1 ) S g p S12 ( x1 ) , где S GR01 ( x1 ) (силы
и момент: равновесие), что и означает наличие шпангоута.
СПИСОК ЛИТЕРАТУРЫ
1.
Гантмахер Ф.Р. Теория матриц. – М.: Наука, 1988. – 548 с.
2. Виноградов Ю.И. Методы вычислений и построение алгоритмов решения краевых задач
строительной механики// Докл. АН СССР, 1988. -Т.298. -№2. -С.308-313.
Отзывы:
Авторизуйтесь, чтобы оставить отзыв