Интегрирование динамических задач механики деформируемого твердого тела обобщенными методами Рунге - Кутты

Yu.V. Nemirovskii and A.P. Yankovskii

Khristianovich Institute of Theoretical and Applied Mechanics SB RAS, Novosibirsk, 630090, Russia

The original numerical procedure of integration of initial-boundary problems of mechanics of deformable solids is developed based on a generalization of Runge – Kutta methods. The efficiency of the proposed numerical method is confirmed by numerous calculations of the elastic and non-elastic dynamics of homogeneous and composite thin-walled structures.

Для анализа динамических процессов в задачах механики деформируемого твердого тела применяются различные методы интегрирования разрешающих систем уравнений. Выбор наиболее эффективного метода тесно связан с моделями поведения рассматриваемых элементов конструкций, идеализацией реальных свойств материала и типа внешних воздействий. Известные методы решения задач динамического деформирования твердых тел разделяются на точные и приближенные (аналитические и численные) методы. Аналитические методы интегрирования разработаны в основном для решения линейно-упругих [1, 2] или жесткопластических задач [3, 4], для тел простой конфигурации при некоторых ограничениях на вид внешних воздействий. Однако решение нелинейных задач динамики точными или приближенными (вариационными [5]) аналитическими методами наталкиваются на существенные трудности, поэтому широкое применение получили численные методы интегрирования.

Среди существующих в настоящее время численных методов интегрирования систем гиперболических уравнений, описывающих динамическое поведение, например, упругих тел, наиболее распространенными являются метод характеристик, метод конечных разностей (МКР), метод конечных элементов (МКЭ) и вариационноразностные методы.

Метод характеристик [6] хорошо разработан для решения линейных и нелинейных одномерных задач о переходных волновых процессах в стержнях, цилиндрических и

конических оболочках [7, 8]. Обобщение его на задачи динамики составных конструкций и многомерные задачи наталкивается на значительные трудности математического и алгоритмического характера [9].

Более универсальным методом решения задач динамического деформирования является МКР. Так как область определения в нестационарных задачах, кроме пространства, включает время, то методы численного решения на основе МКР предполагают дискретизацию определяющей системы уравнений и по этой переменной. Возможны два способа дискретизации: одновременное и последовательное пространственно-временное деление области интегрирования. При одновременном разбиении области определения [10] задача сводится к системе алгебраических уравнений, в которую входят неизвестные на всех временных слоях. Такой способ приводит к определенным трудностям (в частности, увеличивается объем хранимой информации, повышается вероятность зацикливания при решении жесткопластических задач методами линейного программирования и др.), которые ограничивают его применение на практике. Используемые до настоящего времени методы последовательной дискретизации базируются на разбиении области определения в первую очередь по пространственным переменным , сводя исходную задачу к системе обыкновенных дифференциальных уравнений по времени, для интегрирования которых применяются как явные, так и неявные разностные схемы. При решении прикладных задач динамического деформирования тонкостенных конструкций типа оболочек и пластин обычно используются явные трехслойные схемы типа «крест» второго порядка точности относительно шага по времени [9, 11]. Однако эффективность этого метода существенно снижается его условной устойчивостью, что приводит к необходимости использования шагов по времени на 3 … 4 порядка меньших характерного размера сетки по пространственным переменным [9]. Кроме того, недостатки схемы «крест» связаны с трудностями аппроксимации граничных условий и требованием постоянства шага по времени на всем временном промежутке интегрирования [9], что не позволяет без потери порядка точности гибко управлять этим шагом во временном интервале приложения, например, кратковременного высокоинтенсивного нагружения и после него — при движении конструкции по инерции.

Целью настоящего исследования является разработка на основе обобщенных методов Рунге — Кутты [12, 13] принципиально новых численных процедур решения начально-краевых динамических задач механики деформируемого твердого тела, которые базируются на последовательной дискретизации области определения сначала по времени, а затем по пространственным переменным, и лишены указанных недостатков схемы «крест».

В общем случае система N разрешающих уравнений движения деформируемого твердого тела в матричной форме имеет вид:

Rv,t + L(v,t ) = Q (a, t; v, u) + D(v, u), u,t = v(a, t),(1)

где

TTT u = iul’ u 2, ••', uN }  ; v = {vi’ v 2, •••, vN}  ; Q = {Q1, Q 2, •••, QN}

TTT a = {ai, a2, аз} , L ={ L, L,, •••, Ln } ; D ={ Di, D 2, •••, dn }

Здесь также приняты обозначения: R (a) — матрица N X N приведенных объемных плотностей; u, v — обобщенные перемещения и их скорости (в качестве обобщенных перемещений могут выступать и электрические потенциалы в задачах электроупругой динамики и т.п.); α, t — вектор пространственных координат и время соответственно; Q - матричный оператор (в частности, векторная функция), характеризующий обобщенные внешние распределенные массовые нагрузки, действующие на конструкцию, которые в общем случае могут зависеть от v, u (например, при магнитно-импульсном нагружении [11]), сюда же отнесены реакции вязкого, упругого или вязко-упругого основания (может быть использована модель типа Винклера, двух- или многопараметрическая модель основания [14] и др.), а также слагаемые, порожденные геометрической нелинейностью задачи (важной особенностью операторов Qi (1 < i < N) является то, что они не содержат производных по времени от неизвестных векторных функций u, v и не порождают высших производных от этих функций в правой части уравнения (1) по пространственным переменным α); D — матричный дифференциальный оператор эллиптического типа, не содержащий производных по времени и содержащий высшие производные от u, v по переменным α (может быть нелинейным, если рассматривается физически нелинейное поведение материала конструкции); L — линейный матричный дифференциальный оператор по переменным α, перестановочный с оператором дифференцирования по времени t (второе слагаемое в левой части (1) возникает, например, при учете инерции вращения поперечного сечения балок [15], подчиняющихся гипотезам Бернулли, а также кирхгофовских пластин и оболочек или при решении динамических задач электроупругости пьезоэлектрических оболочек [16]).

Индекс после запятой означает частное дифференцирование по времени t .

Для однозначного интегрирования системы (1) необходимы начальные

v ( a , t о ) = v o ( а ) , u ( a t о ) = u o ( ^а )

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

Для численного интегрирования по времени t начально-краевой задачи, соответствующей системе (1), используем обобщенные методы Рунге — Кутты, разработанные авторами [12, 13]. Не будем приводить схемы первого порядка по τ (τ — шаг по времени), так как устойчивыми являются лишь неявные схемы с опережением [12, 13, 17, 18], а перейдем сразу к построению схем второго порядка, которые, как будет показано ниже, требуют примерно тех же вычислительных затрат, что и схемы с опережением.

Используя двустадийный обобщенный метод Лобатто IIIA — метода трапеций (здесь и далее будем использовать терминологию, принятую в работе [19]), получим следующую схему интегрирования системы (1) с точностью второго порядка по τ:

Rv n ^{+ 1} ( а ) + L ( v n ^{+ 1} ( a ) ) = Rv n ( a ) + L ( v n ( a ) ) +

+ 0,5 т [ q ( a , t n + 1 ; v n ^{+ 1} ( a ) , u n ^{+ 1} ( a ) ) + Q ( a , t n ; v ⁿ ( a ) , u n ( a ) ) + + D ( v n ^{+ 1} ( a ) , u n ^{+ 1} ( a ) ) + D ( v n ( a ) , u n ( a ) ) ] ;

un+1 (a) = un (a)+0,5т(vn+1 (a) + vn (a)), n = 0,1, 2,..., где (см. (2)) un (a) = u(а, tn), vn (a) = v(a, tn), tn+1 = tn + t, n = 0,1,2,... (t0 = 0).

Шаг по времени t >  0 может быть переменным ( t_{n + 1} = t_n + t _{n + 1}).

Уравнение (5) запишем в виде

-t D ( v ^{n + 1} ( a ) , u ^{n + 1} ( a ) ) -T Q ( a , t n +p v ^{n + 1} ( a ) , u ^{n + 1} ( a ) ) +

+ 2 Rv ^{n + 1} ( a ) + 2 L ( v ^{n + 1} ( a ) ) = 2 Rv ⁿ ( a ) + 2 L ( v ⁿ ( a ) ) +                               (7)

+t Q ( a , t_n ; v ⁿ ( a ) , u ⁿ ( a ) ) +T D ( v ⁿ ( a ) , u ⁿ ( a ) ) .

Если на n -ом по времени слое векторные функции v ⁿ ( a ) , u ⁿ ( a ) известны, то система уравнений (6), (7) определяет решение на следующем ( n + 1)-ом слое. Недостатком уравнения (7) является то, что для вычисления его правой части необходимо к известным функциям v ⁿ ( a ) , u ⁿ ( a ) применять дифференциальные операторы D и Q . Чтобы исключить из правой части уравнения (7) эти операторы, введем в рассмотрение векторные функции:

P ⁿ ( a ) = -T D ( v ⁿ ( a ) , u ⁿ ( a ) ) -T Q ( a , t_n ; v ⁿ ( a ) , u ⁿ ( a ) ) + 2 Rv ⁿ ( a ) + 2 L ( v ⁿ ( a ) ) ,

_T                                                                    (8)

P ⁿ = { P 1 ⁿ , P n ,..., p n } , n = 0,1,2...

Тогда уравнение (7) примет вид

-t D ( v ^{n + 1} ( a ) , u ^{n + 1} ( a ) ) -T Q ( a , t_{n + 1}; v ^{n + 1} ( a ) , u ^{n + 1} ( a ) ) + 2 Rv ^{n + 1} ( a ) +

+2L(vn+1 (a)) = Pn+1 (a), n = 0,1, 2,..., где правая часть известна и определяется формулой

P ^{n + 1} ( a ) = - P ⁿ ( a ) + 4 Rv ⁿ ( a ) + 4 L ( v ⁿ ( a ) ) ,                                        (10)

полученной в результате сопоставления выражения (8) и правой части (7), которая не содержит операторов D и Q . (Если в разрешающих уравнениях (1) выполняется условие L ( i ) = 0, то правая часть в (9) определяется по рекуррентной формуле (10), вообще не содержащей дифференциальных операторов.)

Запишем выражение для векторных функций в начальный момент времени t ₀ . Учитывая выражения (3), (4), (8), получим

P ⁰ ( a ) = -T D ( v ₀ ( a ) , u ₀ ( a ) ) -T Q ( a , t g ; v ₀ ( a ) , u ₀ ( a ) ) + 2 Rv ₀ ( a ) + 2 L ( v ₀ ( a ) ) . (11)

При нулевых начальных условиях v0 (a) = v0 (a) = 0, u0 (a) = u0 (a) = 0                                             (12)

из(11)следует, что

P ⁰ ( a ) = -T Q ( a , t о ; 0 , 0 ) ,

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

Дальнейшее преобразование уравнения (9) зависит от выбранной модели поведения материала конструкции. Будем различать два случая поведения: вязкопластическое и упруго-вязко-пластическое.

Случай 1 . При вязкопластическом деформировании материала конструкции (или материалов всех фаз, если конструкция композитная) определяющие уравнения строятся на основе известной зависимости напряжения о при одноосном напряженном состоянии от скорости деформации Ё : о = о ( ё ) . При этом важной особенностью уравнения (9) является то, что его левая часть содержит высшие производные по переменным а лишь от скоростей смещений v n ^{+ 1} ( a ) (эти производные порождаются оператором D ). Поэтому в рассматриваемом случае целесообразно исключить из уравнения (9) векторную функцию u n ^{+ 1} ( a ) , используя для этого равенство (6):

-t D n ^{+ 1} ( a ; v n ^{+ 1} ( a ) ) -T Q n ^{+ 1} ( a ; v n ^{+ 1} ( a ) ) + 2 Rv n ^{+ 1} ( a ) +

+2L (vn+1 (a)) = Pn+1 (a), n > 0.

Здесь

D ^{n + 1} ( a ; v ^{n + 1} ( a ) ) = D [ v ^{n + 1} ( a ) , ( u ⁿ ( a ) + T v ⁿ ( a ) /2 ) + t v ^{n + 1} ( a ) /2 ] , Q ^{n + 1} ( a ; v ^{n + 1} ( a ) ) = Q [ a , t n _{+ 1}; v ^{n + 1} ( a ) , ( u ⁿ ( a ) + T v ⁿ ( a ) /2 ) + t v ^{n + 1} ( a ) /2 ] .

Таким образом, в случае вязкопластического деформирования конструкции для определения скоростей смещений ее точек на ( n + 1 )-ом по времени слое необходимо проинтегрировать систему уравнений (14) с известной правой частью (10)-(13) при соответствующих граничных условиях, которые получаются из граничных условий для системы уравнений (1) формальной заменой v на v ^{n + 1} и u на u ^{n + 1} с использованием равенства (6). Если скорости смещений v известны на n -ом и ( n + 1)-ом слоях по времени, то из системы (6) с учетом (4) можно определить перемещения u в момент времени t n _{+ 1}.

Система N уравнений (14) является системой квазилинейных уравнений эллиптического типа и ее можно интерпретировать как систему уравнений равновесия установившейся ползучести рассматриваемой конструкции на вязком (в общем случае, нелинейно-вязком) основании. В силу известного из монографии [20] формального сходства определяющих уравнений установившейся ползучести (в рамках теории течения) и уравнений теории упругопластических деформаций равенство (14) формально совпадает с системой N уравнений равновесия конструкции на упругом (в общем случае, нелинейно-упругом) основании при ее упругопластическом деформировании, если под v ^{n + 1} понимать перемещения точек тела. Поэтому для интегрирования граничной задачи, соответствующей системе уравнений (14), можно использовать известные методы статики или установившейся ползучести.

Для линеаризации операторов (14) используем известные методы: для оператора D n ^{+ 1} ( a ; i) — метод секущего модуля, предложенный в работе [20] для решения задач установившейся ползучести и качественно аналогичный методу переменных параметров упругости, широко используемому при решении упругопластических задач статики [21]; для Q n ^{+ 1} ( a ; •) — идею метода Ньютона [22]. После элементарных преобразований получим разрешающее уравнение на ( к + 1 )-вой итерации

-tDn+1 (a; v^k+^D (a))+ [2R -TQn+1 (a; v"+ (a))] v^k+^D (a) + 2L(vn^ (a)) = = Pn+1 (a)+ t[Qn+1 (a; v"+ (a))- tq. (a; v . (a)) v    (a)], n > ° где Qi'1 — линейный матричный оператор вида

Q " ^{+ 1} ^||d Q ^{n + 1}/ d v ^{n + 1}||,     ( Q ) i   d Q i ^{+ 1}/ d v " ⁺¹,     i , j = 1, 2,..., N .

После такой линеаризации на каждой итерации систему уравнений (15) можно рассматривать как линейную эллиптическую систему уравнений статического деформирования некоторой фиктивной неоднородной конструкции на линейно-упругом винклеровском или многопараметрическом (если L ( • ) ^ 0 ) основании и использовать для ее интегрирования хорошо разработанные для таких задач численные (МКР, МКЭ), вариационные и другие методы [14].

В качестве начального приближения vnO)1 (a) для вектор-функции vn+1 (a) выберем решение на предыдущем по времени слое vi+W v" (a), а в случае L (•) = 0 — векторную функцию vioK a ) = 3v" (a)- R-1Pn (a), получающуюся по формуле Тейлора   vn+1 (a) = vn (a)+Tv,t (a tn)+O (T2)   в предположении, что на предыдущем n-ом по времени слое решение задачи уже известно. При этом учитываем выражения для производной v,t в разрешающих уравнениях (1) и для оператора D (vn (a), un (a)) + Q (a, tn; vn (a), un (a)) (L (• ) = 0) в формуле (8). Здесь R-1 — матрица, обратная матрице R.

В случае вязкопластического деформирования конструкции система разрешающих уравнений динамики (1) является системой квазилинейных уравнений составного типа [23]. То есть, если операторы D и Q в (1) не зависят от u (в случае геометрической линейности), то эта система распадается на две подсистемы, интегрировать которые можно последовательно, причем первая подсистема (1) в этом случае является системой квазилинейных уравнений параболического типа относительно скоростей смещений v. Из работы [24] известно, что для квазилинейных дифференциальных уравнений и систем общая теория устойчивости и сходимости конечно-разностных схем разработана недостаточно, поэтому основным критерием доверия той или иной конечно-разностной схеме служат приближенные решения для тестовых (модельных) задач, аналитические решения которых известны.

Авторам пока не удалось доказать устойчивость численной схемы (6), (9) в общем случае, когда операторы D , Q нелинейны. Но в пользу устойчивости этой схемы говорят физическая корректность (непротиворечивость) результатов многочисленных расчетов, проведенных ими для различных тонкостенных конструкций, и удовлетворительное совпадение результатов с известными аналитическими и апробированными численными решениями [15, 25-30]. В модельном случае линейной вязкости ( о = ие ), а также в предположении о линейности оператора Q в (1) и его независимости от и доказать спектральную устойчивость схемы (7), (9) при L ( • ) = 0 можно, повторив все рассуждения в [12, 13], касающиеся доказательства устойчивости обобщенных методов Рунге — Кутты при решении задачи нестационарной теплопроводности, которая описывается параболическим уравнением, содержащим производную по времени t только первого порядка (подобно первому уравнению (1) при d D / d u = 0 , d Q / d u = 0 ). Из спектральной устойчивости следует устойчивость по начальным данным, из которой, в свою очередь, для двухслойных схем, к которым относится и построенная выше схема, следует устойчивость по правой части [24]. При известных же v n , v n ^{+ 1} схема (6) устойчива [19].

Случай 2 . При линейно-упругом ( о = E е ), нелинейно-упругом или упругопластическом ( а = о ( е ) ), а также при упруго-вязко-пластическом ( о = а ( е , е ) ) деформировании материала конструкции (или материалов всех фаз, если конструкция композитная) левая часть уравнения (9) содержит высшие производные по переменным а лишь от перемещений и n ^{+ 1} ( а ) . Следовательно, в этом случае, используя (6), целесообразно исключить из (9) векторную функцию v n ^{+ 1} ( а ) . Тогда после элементарных преобразований получим:

-t2Dn+1 (а; un+1 (а))-t2Qn+1 (а; un+1 (а)) + 4Run+1 (а) + 4L(иn+1 (а)) = = tPn+1 (а) + 2R[2un (а) + tvn (а)] + 2L[2un (а) + tvn (а)], vn+1 (а) = 2t-1 (un+1 (а)-un (а))- vn (а), n = 0,1,2,...,                              (17)

где

D ^{n + 1} ( а ; u ^{n + 1} ( а ) ) = D [ 2 t^{- 1} u ^{n + 1} ( а ) - ( 2 t^{- 1} u ⁿ ( а ) + v ⁿ ( а ) ) , u ^{n + 1} ( а ) ] ,

Q ^{n + 1} ( а ; u ^{n + 1} ( а ) ) = Q [ а , t n _{+ 1}; 2 t^{- 1} u ^{n + 1} ( а ) - ( 2 t^{- 1} u ⁿ ( а ) + v ⁿ ( а ) ) , u ^{n + 1} ( а ) ] .

Правая часть уравнения (16) с учетом формул (10)-(13) — известная функция.

Таким образом, в рассматриваемом случае деформирования конструкции (тела) для определения перемещений ее точек на ( n + 1 )-ом по времени слое необходимо проинтегрировать систему уравнений (16) с известной правой частью при соответствующих граничных условиях, которые получаются из граничных условий для системы уравнений (1) формальной заменой и на и ^{n + 1} и v на v ^{n + 1} с использованием равенства (17). Если перемещения и известны на n -ом и ( n + 1)-ом слоях по времени, то из системы (16), (17) можно определить скорости смещений v точек конструкции в момент времени t_{n + 1} .

Система N уравнений (16) является системой линейных или квазилинейных уравнений эллиптического типа и ее можно интерпретировать как систему уравнений равновесия упруго или упругопластически (в рамках деформационной теории) деформируемой конструкции на упругом (в общем случае нелинейно-упругом) основании.

Для линеаризации в (16) оператора D n ^{+ 1} ( a ; i ) используем метод переменных параметров упругости [21], а для линеаризации оператора Q n ^{+ 1} ( a ; i ) — идею метода Ньютона [22]. После элементарных преобразований получим разрешающее уравнение на ( к + 1)-й итерации:

2 т\п + 1       -,n + 1                      22-т n + 1       ..n + 1            n + 1

• -Т D ( ^a ; u ( к + 1) ^{( a )} ) ⁺ ^ ^{4 R — T Q} u ( ^a ; u ( к ) ^{( a )} ) J u ( к + 1) ^{( a ) +}

• + 4 L ( u nk + 1) ( a ) ) = T P ^{n + 1} ( a ) + 2 R [ 2 u ⁿ ( a ) + t v ⁿ ( a ) ] + 2 L [ 2 u ⁿ ( a ) + t v ⁿ ( a ) ] +      (18)

• ^+T 2 [ ^{Q n + 1} ( a ; u "n_k + ¹ ( ^a ) ) ^{- Q} n ^{+ 1} ( ^a ; u ⁿ n_fc + 1 ( ^a ) ) u nk + 1 ( ^a ) ] , n = 0,1, ², -^,

где Q U ^{+ 1} — линейный матричный оператор вида

Q U ^{+ 1} =|d Q ^{n + 1}/ d u ^{n + 1}|, ( Q n ^{+ 1} ) i =d Q n ^{+ 1}/ d u j ⁺¹,    i , j = 1,2,..., n .

После такой линеаризации на каждой итерации систему уравнений (18) можно рассматривать как линейную эллиптическую систему уравнений равновесия некоторой фиктивной неоднородной конструкции на винклеровском или многопараметрическом (если L ( • ) ^ 0 ) основании и использовать для ее интегрирования различные численные (МКР [17, 18, 31, 32], МКЭ [33]), вариационные и другие методы [14]. В некоторых линейных задачах упругого деформирования конструкций систему уравнений (16) на каждом шаге по времени можно проинтегрировать аналитически [17, 18].

В качестве начального приближения u n ^)¹ ( a ) для вектор-функции u ^{n + 1} ( a ) выберем функцию