Численное моделирование динамических процессов в твердых телах на основе схем повышенной точности

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

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

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

Однако существует ряд задач, для которых требования к точности решения оказываются достаточно жесткими. Одной из таких задач, к примеру, является прямая задача сейсмики, характерной особенностью которой можно считать достаточно малую зону действия источника возмущений по сравнению с масштабами всей расчетной области. Кроме того, сильно осциллирующий характер источника накладывает ограничения на шаг по времени, что, в свою очередь, приводит к необходимости достаточно мелкой дискретизации пространства. Для ее решения традиционно применяется полуаналитический метод, основанный на сведении исходной многомерной задачи к семейству одномерных краевых задач с последующим восстановлением решения [1, 2]. Попытки же решения этой задачи с использованием классических разностных схем наталкиваются на серьезные трудности.

Анализ задачи и тестовые расчеты показывают, что ее численное решение на основе разностной схемы принципиально возможно только в случае, если положительно решены следующие проблемы:

— алгоритм дает монотонное решение и, в частности, не искажает его в окрестности оси симметрии (задачу, как правило, приходится Ь/ рассматривать в осесимметричной постановке);

— метод допускает естественную формулировку физических краевых условий (отсутствуют сложности, связанные с введением фиктивных слоев ит. д.);

— схема не обладает искусственной диссипацией ни на «продольных», ни на «поперечных» волнах (опыт показывает, что наличие даже малой искусственной диссипации очень скоро приводит к полному затуханию возмущения);

— алгоритм допускает обобщение на случай существенно неоднородной среды;

— условия на границах выделенного для численного решения объема эффективно обеспечивают «неотражение» от этих границ.

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

Существенное влияние на методы численного решения задач механики твердого тела оказали подходы, разработанные при решении задач газовой динамики и механики вязкой жидкости. Исторически, вычислительные методы в этих областях стали применяться раньше, и в настоящее время здесь накоплен существенно больший опыт. Достаточно полный обзор существующих методов численного решения задач гидрои газодинамики, анализ и удобная классификация конечно-разностных схем по определенным признакам даны в работах [3 — 9]. В то же время, как отмечается в [10], существенная сложность определяющих уравнений твердого тела и специфика этих задач не позволяют непосредственно переносить результаты из области гидромеханики на задачи твердого тела.

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

Численному решению задач динамической теории упругости и пластичности посвящено большое число работ. Их подробный обзор и анализ различных подходов можно найти, например, в работах [10, 12−14]. Существующие методы решения задач динамики твердых тел достаточно условно можно представить в виде трех направлений [10]:

— методы конечных элементов;

— характеристические и сеточно-характеристические методы;

— сеточные или конечно-разностные методы.

Следует сказать, что все три подхода не противопоставляются, и их взаимное проникновение все более заметно в последнее время.

Под методами конечных элементов понимают подходы, основанные на дискретизации расчетной области и формировании конечных соотношений между искомыми величинами (действующими в узлах си-, / лами и их перемещениями) на основе законов механики в вариационной / форме, минуя стадию формулировки краевых задач для систем дифференциальных уравнений. Такой подход дает определенные преимущества при описании процесса деформирования тел со сложной геометрией. Метод надежно зарекомендовал себя для решения статических задач и сейчас интенсивно развивается в области исследования нестационарных процессов в деформируемых твердых телах. Среди отечественных работ этого направления необходимо отметить, например, работы [15−26].

Как сочетание и обобщение методов конечных элементов и вариационно-разностных, можно упомянуть дискретно-вариационный метод [27, 28], разработанный для исследования нестационарной динамики в слоистых и композиционных средах. Конечно-элементный подход активно развивается за рубежом. В качестве примера здесь можно указать работы [29 — 34].

Детальному изложению, подробному обзору и анализу характеристических и сеточно-характеристических методов посвящена монография [9]. Эти подходы основаны на записи системы дифференциаль- / ных уравнений в характеристической форме с последующей их конечно-разностной аппроксимацией. Различают прямой и обратный характеристический метод [36, 37].

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

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

В обратном характеристическом методе решение для всей области вычисляется в точках, отвечающих одному и тому же шагу по времени. При этом используется конечно-разностная аппроксимация соотношений на характеристических плоскостях, касательных к характеристическим конусам, выходящим из этой точки на предыдущий / («нижний») слой по времени. Метод требует интерполирования на предыдущем слое, при этом расширяется область зависимости разностной задачи. В некоторых работах [9, 38] метод называется сеточно-характеристическим.

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

Среди работ, посвященных применению сеточно-характеристи-ческих методов для решения динамических задач деформирования упругих и упругопластических тел, можно указать работы [46 — 50].

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

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

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

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

Среди практически используемых неявных разностных схем наибольшее применение получили схемы расщепления [51 — 56] и схемы, основанные на методе дробных шагов [57 — 66]. Среди решенных с использованием этих схем задач можно упомянуть [67 — 69]. Достаточно эффективными оказываются подходы, допускающие использование комбинированных схем, в частности, схем явных в одном направлении и неявных в другом.

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

Метод выделения разрыва, позволяющий рассчитывать разрывные решения без размазывания скачков, был предложен С. К. Годуновым [ТО, 71] и основан на использовании подвижных сеток. В расчетной области, с помощью известного соотношения на скачке, выделяется поверхность разрыва. Течение за фронтом является гладким, и расчет его по явным или неявным схемам не вызывает больших проблем. Метод широко и эффективно используется при расчете газодинамических течений [71], однако, применение его чрезвычайно затруднительно в задачах, для которых характерно присутствие различных поверхностей разрыва, положение и конфигурация которых неизвестны.

При выборе численной схемы сквозного счета для исследования распространения ударных волн, их взаимодействия и т. д., следовало бы отдать предпочтение схемам повышенного порядка точности, позволяющим более тонко описывать картину течения, экономить время решения задач на ЭВМ. С помощью широко известных схем второго и выше порядка точности [71 — 75] либо их модификаций решается значительное число задач газовой динамики. Среди примеров решения неодномерных задач механики твердого тела с помощью методов повышенного порядка аппроксимации следует отметить схему [50] и работы [76, 77].

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

Одним из способов подавления нефизичных эффектов является процедура введения в сами уравнения дополнительных членов, которые принято называть искусственной вязкостью [71, 8]. Иногда к полученному решению применяют некоторое сглаживание [71], которое может быть каким-то образом связано с характером решения [79], либо применяться вообще безотносительно к уравнениям. Анализ диссипативных и дисперсионных свойств таких схем можно найти в [80].

Эффективный способ построения монотонных схем, имеющих на гладких решениях второй порядок аппроксимации, предложен Борисом и Буком [81, 82], и развит в последующих работах, например в [83 -85]. Метод названный ими ЕСТ-метод (метод коррекции потоков) представляет собой нелинейное консервативное сглаживание, имеющее два характерных этапа — введение в схему диффузионного оператора и последующее исключение диффузии. Метод применялся в основном для решения задач газодинамики. В [76, 77], на основе идеи коррекции решения, разработан гибридный конечно-разностный алгоритм для расчета динамических процессов деформирования оболочечных конструкций.

Еще одним из способов борьбы с немонотонностью схем высокого порядка аппроксимации является процедура монотонизации, представляющая собой подстройку численного алгоритма в зависимости от характера решения на предыдущем временном слое. В результате строится нелинейная разностная схема, сохраняющая высокий порядок точности. Целое семейство таких схем с переменным шаблоном получено на основе принципа минимальных значений производной, предложенного в [86 — 88]. Их обзор можно найти в [89 — 92].

На основе аналогичного подхода построены монотонные схемы второго и выше порядка аппроксимации в работах [93, 94]. К этому же семейству методов можно отнести опубликованные в [95 — 100 ] и изложенные в настоящей диссертационной работе схемы второго порядка точности, строгая монотонность которых обеспечивается специальным выбором входящих в схему параметров — констант диссипации, в зависимости от характера решения на нижнем слое по времени.

Недостатков, связанных с немонотонностью решения, лишен метод, предложенный С. К. Годуновым для расчета одномерных [70] и многомерных [11] задач газовой динамики. Метод представляет собой двухшаговую схему типа предиктор-корректор. На каждом слое решение рассматривается как кусочно-постоянное, а для вычисления некоторых вспомогательных («больших») величин на промежуточном этапе используются формулы распада произвольного разрыва. Схема обладает свойством монотонности, но имеет только первый порядок точности. В [11] выписана и схема решения плоской динамической задачи теории упругости в декартовых координатах. На основе метода С. К. Годунова и его модификаций получено решение ряда задач динамической теории упругости как в плоской геометрии, так и в криволинейных системах координат [101 — 104].

Для численного решения двумерных задач динамической теории упругости Г. В. Ивановым была предложена идея использования нескольких локальных аппроксимаций неизвестных функций полиномами [105]. В отличие от классических конечно-элементных подходов процедура нескольких аппроксимаций для каждой из искомых функций дает возможность сформировать достаточную для обеспечения монотонности численного решения искусственную диссипацию с одновременным расщеплением многомерной задачи на ряд независимых одномерных задач по пространственным направлениям. Содержащиеся в одномерных задачах параметры — константы диссипации — позволяют регулировать в получаемых схемах величину искусственной вязкости. При этом процедура решения каждой из этих одномерных задач является независимой и может быть как явной, так и неявной.

При частном выборе констант диссипации, для случая регулярных сеток, получающаяся явная схема полностью совпадает со схемой С. К. Годунова. В то же время, одно из преимуществ подхода состоит в возможности построения схем, обладающих лучшими диссипативными свойствами по сравнению с методом «распада разрыва». Существенно, что при этом не происходит увеличение числа арифметических опера-: ций, т. е. помимо улучшения диссипативных характеристик возрастает экономичность схемы. Метод получил дальнейшее развитие [106 -134]. На его основе исследован ряд задач упругого и упругопластиче-ского динамического деформирования твердых тел.

Настоящая работа посвящена построению эффективных численных алгоритмов повышенной точности интегрирования одномерных и многомерных задач динамической теории упругости и моделированию на их основе динамических процессов в твердых телах. В ней обобщены результаты исследований, проведенных за период с 1977 г. по 1997 г. и опубликованных в работах [95 — 100, 114 — 134]. Исследования проводились в рамках госбюджетных тем Вычислительного центра СО РАН в г. Красноярске (per. N 01.9.20 15 428 и per. N 01.96 0.0 5 372), хоздоговоров (НИИ ПМ N 9206) и были поддержаны грантами Краевого фонда науки (1Е0133 1994;95 гг. и 5F0061 1996 г.) и РФФИ (N 97−01.

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

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

Основные результаты работы состоят в следующем.

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

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

2. Получено численное решение ряда одномерных задач динамики:

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

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

— на основе численного моделирования процессов распространения электромагнитных волн в слоистых анизотропных диэлектриках установлена зависимость между спектрами и структурой волн и ориентацией оптической оси жидкого кристалла;

— предложена эффективная схема и проведено численное моделирование процессов распространения плоских волн в слоистых трансверса-льно-изотропных упругих материалах.

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

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

5. На основе построенного численного алгоритма решен ряд многомерных задач динамического деформирования неоднородных упругих сред. А именно:

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

— создан комплекс программ для моделирования процессов распространения сейсмических волн в слоистой вертикально-неоднородной среде;

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

ЗАКЛЮЧЕНИЕ

.