автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Численное моделирование динамических процессов в твердых телах на основе схем повышенной точности

ГОСУДАРСТВЕННЫЙ КОМИТЕТ РОССИЙСКОЙ ФЕДЕРАЦИИ ПО ВЫСШЕМУ ОБРАЗОВАНИЮ КРАСНОЯРСКИЙ ГОСУДАРСТВЕННЫЙ ТЕХНИЧЕСКИЙ УНИВЕРСИТЕТ

•« р}. А1". -

""* 3 /, ]; | I

На правах рукописи

Вогульский Игорь Олегович

УДК 539.3

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ДИНАМИЧЕСКИХ ПРОЦЕССОВ В ТВЕРДЫХ ТЕЛАХ НА ОСНОВЕ СХЕМ ПОВЫШЕННОЙ ТОЧНОСТИ

05.13.18 - теоретические оснрвы математического моделирования, численные методы и комплексы программ

АВТОРЕФЕРАТ диссертации на соискание ученой степени доктора физико - математических наук

Красноярск - 1998 г.

Работа выполнена в Институте вычислительного моделированы СО РАН, Красноярск

Официальные оппоненты: доктор физико - математических наук, профессор доктор физико - математических наук, профессор доктор физико - математических наук, профессор

Р. П. Федоренко Б. Д. Аннин Е. А. Новиков

Ведущая организация:

Институт вычислительной математики и математической геофи зики СО РАН , Новосибирск

. Защита диссертации состоится " " 199-^ г. :

часов на заседании диссертационного совета Д 064.54.01 по за щите диссертаций на соискание ученой степени доктора наук в Красно ярском государственном техническом университете по адресу: 660074 г. Красноярск, ул. Киренского 26.

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

Автореферат разослан "___199_£. г.

Ученый секретарь диссертационного совета Д 064.54.01 д. т. н., профессор л /

Ловчиков А. Н.

Общая характеристика работы

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

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

Большей части указанных недостатков лишена, например, хорошо известная схема С. К. Годунова, имеющая, однако, только первый порядок аппроксимации. Для ряда задач использование этой схемы в вычислительном эксперименте вполне приемлемо и достаточно. Однако существуют задачи (к примеру, задачи сейсмики), решение ко-

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

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

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

Научная новизна. Новыми в диссертации являются:

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

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

3. Явные монотонные схемы решения двумерных (плоских и осесим-метричных) задач динамической теории упругости с минимальным (по

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

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

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

Практическая ценность.

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

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

Результаты диссертации вошли в лекционные курсы кафедры "Математическое моделирование в механике" Красноярского госуниверситета и кафедры "Диагностика и безопасность технических систем и инженерных сооружений" Красноярского государственного техниче-

ского университета.

Апробация работы. Результаты, вошедшие в диссертацию, докладывались на X, XI, XIII Всесоюзных и XV Международной конференциях "Численные методы решения задач теории упругости и пластичности" (Красноярск, 1987 г., Волгоград, 1989 г., Новосибирск, 1993 и 1997 гг.), Всесоюзной школе молодых ученых "Вычислительные методы и математическое моделирование" (пос. Шушенское, 1986 г.), УШ Всесоюзном симпозиуме по распространению упругих и упругопласти-ческихволн (Новосибирск, 1986 г.), Всесоюзной школе молодых ученых "Численные методы механики сплошной среды" (Абакан, 1989 г.), I и П Сибирских школах по современным проблемам механики деформируемого твердого тела (Новосибирск, 1988 г.; Якутск, 1990 г.), Ш и IV Республиканских школах - семинарах молодых ученых и специалистов по теоретической и прикладной гидромеханике (Алушта, 1988 и 1990 гг.), III Всесоюзной Школе молодых ученых по "Численным методам механики сплошной среды (пос. Дюрсо, 1991 г.), УШ Всесибирской школе по вычислительной математике (пос. Шушенское, 1993 г.), Международной конференции "Проблемы обеспечения качества изделий в машиностроении" (Красноярск, 1994 г.), Международном симпозиуме "Аспекты использования углеводородных и минеральных ресурсов" (Красноярск, 1995 г.), Международной конференции "Математические модели и численные методы механики сплошной среды" (Новосибирск, 1996 г.), II Сибирском конгрессе по прикладной и индустриальной математике (Новосибирск, 1996 г.) и Международной конференции "Математические модели и методы их исследования (Красноярск, 1997 г.), а также на семинаре "Проблемы математического и численного моделирования" ВЦ СО РАН в г. Красноярске (руководитель - акад. Ю. И. Шокин, 1986 г.; член-корр. РАН В. В. Шайдуров, 1991, 1994 гг.), расширенном семинаре отдела "Нелинейных задач механики" ИВМ СО РАН (руководитель - профессор В. К. Андреев, 1997 г.), семинаре по механике деформируемого твердого тела в Институте гидродинамики им. М. А. Лаврентьева СО РАН (руководитель - профессор О. В. Соснин, 1986, 1997 гг.), семинаре "Большие задачи математической физики" ВЦ СО РАН в г. Новосибирске (руководитель - член - корреспондент РАН А. Н. Коновалов, 1986 г.), на семинаре по динамике сплошной среды

Института проблем механики РАН (руководитель - профессор В. Н. Кукуджанов, 1986, 1991, 1998 г.г.), семинаре кафедры волновой и газовой динамики МГУ (руководитель - акад. РАН Е. И. Шемякин, 1998 г.), объединенном семинаре Института вычислительных технологий СО РАН и Института вычислительной математики и математической геофизики СО РАН (руководители - акад. Ю. И. Шокин, проф. В. М. Ковеня и член-корр. РАН А. Н. Коновалов, 1997 г.), и на семинаре Института прикладной математики им. М. В. Келдыша РАН (руководитель - проф. А. В. Забродин, 1998 г.).

Публикации. Основные результаты диссертации опубликованы в работах [1] - [27], в том числе в монографии [21].

Объем работы. Диссертация состоит из введения, пяти глав, заключения и списка литературы, изложенных на 282 страницах, включает 79 рисунков. Список цитируемой литературы на 14 страницах содержит 156 наименований.

Краткое содержание работы

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

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

нии задач гидро- и газодинамики, отраженные в работах А. А. Самарского, Н. Н. Яненко, Ю. И. Шокина, О. М. Белоцерковского, С. К. Годунова, А. В. Забродина и др. Численному решению задач динамической теории упругости и пластичности посвящено большое число работ. Их подробный обзор и анализ различных подходов наиболее полно представлен в работах В. Н. Кукуджанова и В. И. Кондаурова. Существующие методы достаточно условно можно представить в виде трех направлений: - методы конечных элементов; - характеристические и сеточно-характеристические методы; - конечно-разностные методы. Эти три подхода не противопоставляются и их взаимное проникновение все более заметно в последнее время.

Метод конечных элементов получил большее распространение при решении задач статики. Для решения динамических задач он активно развивается за рубежом (К. Bathe, J. Hallquist и др.). Среди отечественных работ этого направления наиболее известны работы Н. Г. Бураго, В. Г. Баженова. Характеристические и сеточно-характеристические подходы основаны на записи исходной системы в характеристической форме с последующей их конечно-разностной аппроксимацией и наиболее полно отражены в работах А. С. Холодова, К. М. Магомедова, В. Н. Кукуджанова, И. Б. Петрова.

Среди сеточных методов различают явные и неявные схемы. Среди практически используемых неявных схем наибольшее применение получили схемы расщепления и схемы, основанные на методе дробных шагов, развитие которых связано в первую очередь с работами А. Н. Коновалова, Н. М. Горского. Значительное число явных схем имеет в своей основе метод "распадаразрыва", предложенный С. К. Годуновым. Большое развитие для решения динамических задач теории упругости этот метод получил в работах В. Г. Чабана, В. К. Римского. П. О. Сабодаша и др.

Достаточно внимания уделено вопросу использования схем повышенного порядка точности и, тесно связанной с этим, проблеме монотонности решения. Среди способов борьбы с паразитными осцилляци-ями отмечена идея введения искусственной вязкости (П. Лаке, Б. Вея-дрофф), подходы на основе коррекции потоков (Д. Бук, Дж. Борис, А. И. Жмакин, А. А. Фурсенко и др.), различные способы монотониза-

ции (В. П. Колган, Б. вал Лир, А. В. Родионов). Среди используемых при решении задач динамической теории упругости методов высокого порядка точности известен подход Р. Клифтона и гибридные схемы, исследованные в работах А. А. Гурьянова.

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

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

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

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

и = ио + ЩТ), а = а0 + а-11], и'— и'0 + и'^, а'~а'0 + а\^ (1) которые удовлетворяют системе уравнений

ва".

ди да' да 20и'

Р~Ш = ~дх'' Л = рс'дх

дх'

(2)

Величины и, а удовлетворяют начальным, а и , а' - граничным условиям; г] и £ - локальные координаты в w.

Формулы пересчета на следующий слой по времени в силу (1),(2) принимают вид

= uj+1 + - = "Щ + - (3)

а необходимые для определения величин Uj, Oj дополнительные уравнения формулируются на основе записи разностного аналога закона сохранения механической энергии

г,рди2 ди\ г _ , / д(и'а') ,

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

. т да1, да' , , 2 гди'.ди'

О = К - - + (-о ~ ^ - рс _)_.

В качестве дополнительных принимаются уравнения

, г да' 71 Зет' <Эи' , т ди' 74 ди' да' . .

и°~и}+1~2~дх = 2 (4)

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

Из множества полученных схем выделяется двухпараметрическое семейство явных схем

1 71 + Д

«Я-1 = £ К1 +72Н+§ + ~ + 1 ~

1 74 4- Д

4+1 = ¿К1 ~ + (1 +72)а,ц)+рс 2 ~ «я-*)»

722 = 1 -(71 + Д)(74 + Д),

включающее как частный случай и схему С. К. Годунова (7! = 74 = 1 —Д) и схему, имеющую нулевую мощность искусственной диссипации

(71 = 74 = 0), что в данном случае соответствует схеме второго порядка точности на гладких решениях. Удается выделить область изменения параметров, в которой расчеты показывают монотонный профиль решения, в то время как диссипация существенно меньше, чем по схеме "распада разрыва". Показано, что при условии диссипативности имеет место сходимость построенного приближенного решения к точному в энергетической норме. Алгоритм иллюстрируется тестовыми расчетами.

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

" „ , "+1 « = е(«2+«Ь)А(О, ..., «' = гима ...,

*=0 ¿=0

где Рк{0 ~ полиномы Лежандра степени к.

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

В четвертом параграфе предложенный выше подход обобщается на случай решения одномерных систем гиперболических уравнений общего вида

Л1н+в1Гх+Си = 0>

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

форме

дУ дУ1

+ + — 0. (5)

Принципиально новым элементом при построении схем для (5) является процедура независимой аппроксимации младших (недифференциальных) членов уравнений

У = Ц + IV/, V' = У0' + V".

В итоге формулы пересчета решения на верхний по времени слой принимают вид

+г) = - 1мут - V,] -

а величины с целочисленными индексами и константы V" к определяются из сформулированных на основе разностного аналога энергетического тождества уравнений

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

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

12

устойчивой в энергетической норме. Схема I (С. К. Годунова)

\¥1 = -т1, Г2 = г(2т - 1)7, = г(27 - 1)/

Схема II

И? + Г2 = 0 Г1 = 7убу, ^=0, Г2 = О, 1У2 = 0, 7кк>0. Схема III

Г1=7У«У. Г2 — т1, \¥2 = 0.

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

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

В первом параграфе главы рассматриваются безразмерные уравнения нестационарного деформирования пластины с постоянными по толщине деформациями сдвига (модель С. П. Тимошенко).

дЧ/ дЯ _ дУУ

Ж + _х дм _ ды т ~ вх 6 т ~ дх'

г>0, 0<*<1. £г = /г2/12 ь2.

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

13

применение даже устойчивых схем I и II может давать быстрое накопление ошибок округлений.

Показано, что схема III с независимой аппроксимацией младших членов, построенная в §4 предыдущей главы, обладает свойством сильной устойчивости (хорошей обусловленности) при решении задач для уравнений С. П. Тимошенко. На примере решения задачи об импульсном деформировании пластины показано, что в схеме III влияние искусственной диссипации приводит к тому, что при неограниченном росте времени значения скорости прогиба стремиться к нулю, в то время как структура диссипации в схемах I и II такова, что с течением времени напряжения (усилия и моменты) стремятся к постоянным значениям, близким к статическому решению, а скорости - к постоянным, но не равным нулю значениям, приводящим ж неограниченному росту прогиба пластины.

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

ди _ д{г<тг) ^ даг _ du ^^и dav _ ^ди и dt дг dt дг г' dt дг г'

Проблемы, возникающие при создании алгоритмов для таких задач, связаны с тем, что коэффициенты младших членов уравнений этого класса имеют особенность в точках оси симметрии. Численные расчеты показывают, что применение в этом случае схем, лишенных недостатков в плоских задачах (схемы I и II §1.4), может приводить к решениям, обладающим дополнительными (нефизичными) локальными экстремумами в окрестности оси симметрии. При решении двумерной осесимметричной задачи это, например, приводит к тому, что фронт плоской волны, распространяющейся вдоль оси симметрии, деформируется вблизи этой оси. Оказывается, что независимая аппроксимация младших членов (схема III) позволяет сформулировать схемы, свободные от указанных недостатков.

С целью тестирования алгоритма проведено сравнение численных результатов, полученных на основе построенной схемы и аналитического (асимптотического) решения задачи о схлопывании цилиндри-

ческой упругой волны.

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

В математическом плане необходимо получить численное решение системы уравнений Максвела

тоЬЕ = --~, го1= В = у.Н, б = ёЁ,

с от с at

где ё и ¡1 - тензоры диэлектрической и магнитной проницаемости. В

данном случае проблема сводится к решению двух самостоятельных

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

ЭЕХ ^д Ну ЭЕу цдНх

дг с с^ дг с с^ '

дНу а{0) дЕх дНх е± дНх

дг с д1 ' дг с '

где

а{9) =

€± зт20 + £ц соэ^ 6' Существенная и сильно меняющаяся неоднородность среды приводит к необходимости включения в алгоритм автоматической дискретизации расчетной области, согласованной со свойствами конкретной слоистой структуры. Контакт рассматриваемой сверхрешетки с вакуумом моделируется путем включения слоев вакуума в расчетную область и формулировки на концах области "неотражающих" граничных условий.

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

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

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

ди да2 dv drzx

P~dt=lh' P~di~~~d7

да* _ лди i Bdv дт" _ вди i cdv dt dz dz' dt dz dz'

A = a cos4 a + 26 sin2 a cos2 a + с sin4 a + /i2 sin2 2a, С = - 2b + c) sin2 2a + /z2 cos2 2a,

В = — cos2 a — b cos 2a — с sin2 a) sin 2a + Ц2 cos 2a sin 2a,

_ (1 ~ vi)E*l ^ U2EiE2 й Е2(1-рх)-2Е1иГ E2{\ - vi) - 1ExvV _ E,{E-2 — v%Ei)_

Здесь a - угол наклона кристаллографической оси к оси 2, /i2, щ. v2, Ei, E'¿ - упругие константы материала. Алгоритм численного решения таких уравнений описан в §4 первой главы и основан на приведении системы к каноническому виду. Однако, в случае системы большой размерности, такая процедура сопряжена со значительными техническими трудностями. Поэтому в работе, с одной стороны, строится и иллюстрируется примерами решения ряда задач метод, сформулированный для систем общего вида в §1.4, а с другой стороны предлагается итерационная процедура решения полученной гиперболической

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

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

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

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

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

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

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

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

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

В элементарном параллелишшеде ш в качестве приближенного решения принимаются функции

и=и° -f- и1!], v = v° + v1т], ах = + air),... u' = u0 + uif, v' = v'0 + v[£, dx = cr'xQ +

п // и н п // // и //

и ^Uq + щС, V ^Uo-Ъг^с, ory = <T!/0 + fffflC,... которые удовлетворяют уравнениям

ди^да^ Sir;, dv _ дт'ху да"у

P8t ~ дх ду ' Pdt дх ду'

дах 1/ди дь" . дау п/л^м ду'\

дтху 2 ди дь ~дГ = Р ~8у ~дх Функции без штрихов удовлетворяют начальным условиям, функции с одним штрихом - краевым условиям на гранях, перпендикулярных оси х, с двумя штрихами - на гранях, перпендикулярных оси у.

Для формулировки необходимых дополнительных уравнений выписывается разностный аналог закона сохранения энергии.

Алгоритм основывается па расщеплении полной задачи на четыре независимые одномерные задачи вида (4), причем каждая из этих одномерных задач содержит максимальный набор констант диссипации, так что в целом в алгоритм заложен достаточно большой произвол.

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

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

по сравнению с ограничением в схеме "распада разрыва"

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

- с аппроксимацией как усилий, так и напряжений в соответствующих элементарных объемах;

- и с аппроксимацией младших, недифференциальных членов уравнений.

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

19

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

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

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

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

^ 2 „ dudv

д = ... + а(_) _2г__ + /3(_) +

для неотрицательности которых необходимо выполнение неравенств а > 0, (3 > 0, а/3 > г2,

что исключает возможность выбирать ais. ¡3 нулевыми.

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

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

20

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

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

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

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

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

Для иллюстрации работы алгоритма решен ряд задач об удар-

21

Эи. 2 ду> '

(6)

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

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

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

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

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

тры диссипации а и /? для квадратных ячеек принимают значения

а = р = }1/ср — т и а = Р — Н/2св — т

для продольных и для поперечных волн соответственно. Это означает, что для реальных материалов, у которых значение коэффициента Пуассона близко к 1/3 (т.е. отношение с8/ср близко к 1/2), на одной и той лее квадратной сетке, выбором к/ср = т удается получить численное решение, не обладающее искусственной диссипацией на продольных и с практически нулевой диссипацией на поперечных волнах.

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

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

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

Работоспособность предложенного алгоритма проиллюстрирована примерами решения ряда модельных задач.

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

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

23

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

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

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

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

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

теризующую отражение возмущения от каждой из границ раздела и в итоге - временной разрез рассматриваемой среды.

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

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

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

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

оо

Д = Jl<p(t) - v*(t))2dt, (7)

о

<p(t) = 1 * /(f) = j f(T)dr, cp*(t) = 1 * Г it) = J f*(r)dT.

о 0

Процедура решения состоит в следующем. Заранее в достаточ-

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

Заключение

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

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

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

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

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

- с помощью введения независимой аппроксимации младших членов построены схемы решения динамических задач с осевой и сфериче-

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

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

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

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

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

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

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

телом существенно пространственной формы;

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

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

Основные результаты диссертации опубликованы в следующих работах:

1. Вогульский И. О. Построение монотонной схемы решения задач для гиперболических уравнений - Красноярск, 1982. - 18 с. - (Препринт / ВЦ СО АН СССР, N 26)

2. Вогульский И. О. Построение монотонных схем решения одномерных задач для гиперболических уравнений па основе аппроксимации полиномами Лежандра. - В сб.: Труды I Всесоюзной школы-семинара по многомерным задачам МСС. - Деп. в ВИНИТИ, 1983 - N 4623/83. - С. 59 - 65.

3. Вогульский И. О. и др. Математическое моделирование процессов трехмерного проникания. - Красноярск-Новосибирск, 1983 (отчет N1.1.16.1, per. N81017684).

4. Вогульский И. О. Об одном семействе явных схем решения задач динамики упругих тел на основе аппроксимации линейными полиномами. - Деп. в ВИНИТИ, 1985. - N 687 - 85 Деп. - 35 с.

5. Вогульский И. О. Повышение точности решения плоских динамических задач упругости в рамках аппроксимации линейными полиномами. - Деп. в ВИНИТИ, 1986. - N 65 -В86. - 52 с.

6. Вогульский И. О. Схемы решения двумерных задач динамики упругих тел // В сб.: Доклады VIII Всесоюз. конф. по распространению упругих и упругопластических волн - Новосибирск, 1986. - С.7-12

7. Вогульский И. О. Схемы повышенной точности решения задач динамики упругих тел. - Новосибирск, 1986. - 24 с. (Автореферат канд. диссерт.)

8. Вогульский И. О. Об одном семействе явных монотонных схем

решения задач динамики упругих тел. - В сб.: Численный анализ и пакеты прикладных программ. - Красноярск, 1986. -С. 42 - 54.

9. Вогульский И. О. Мопотонная схема второго порядка решения задач динамики упругих тел. - Деп. в ВИНИТИ, 1986.- N 64 - 86 Деп.-56 с.

10. Вогульский И. О. Асимптотическое поведение цилиндрических ударных волн вблизи оси симметрии - Красноярск, 1988. - С. 16 - 20. (Препринт / ВЦ СО АН СССР, N 2)

11. Вогульский И. О. Об одном алгоритме решения двумерной динамической задачи механики деформируемого твердого тела. - В сб.: Численные методы механики сплошной среды. Ч. 2. - Красноярск, 1989. - С. 42 - 45.

12. Вогульский И. О. Об одной схеме расщепления решения двумерной задачи динамики твердого тела. - Деп. в ВИНИТИ, 1989. - N 1816 -В89. - 40 с.

13. Вогульский И. О. Численное моделирование распределенного ударного воздействия цилиндрических ударников на упругую плиту // В сб.: Доклады II Всесибирской. школы по современным проблемам МДТТ - Новосибирск, 1990. - С.7-9

14. Вогульский И. О. Об одном алгоритме решения двумерной осе-симметричной задачи динамики твердого тела.- Деп. в ВИНИТИ, 1990, N 5234 - В90. 35 с.

15. Вогульский И. О. Численное моделирование распределенного ударного воздействия на упругую плиту // Деформирование и разрушение современных материалов и конструкций (Динамика сплошной среды), 1991. - N 103.- С.ЗО - 35.

16. Bogulskii I. О. On numerical method for solving of two-dimensional problems of deformable solids dynamics //Modelling & Control, В., AMSE Press,France,1993,- Vol. 47, N 4 - P.21 - 40.

17. Bogulskii I. O. A monotonicity schemes of second-order accuracy for solving of problems of deformable solids dynamics.// Modelling & Conrtol, В., AMSE Press, France, 1994. - Vol. 53, N 2, - P. 19 - 28.

18. Bogulskii I. O. The use of second-order polynomials for local approximation of solution of two-dimensional problem of dynamic elastic theory //Modelling & Control,B, AMSE Press,1994. - Vol.53, N 2 - P.29-38.

19. Анисимов С. А., Вогульский И. О. Алгоритм независимой аппроксимации недифференциальных членов при численном решении краевых задач для систем гиперболических уравнений // Динамика сплошной среды.- Новосибирск, 1994. - Вып. 109.- С. 34 - 48.

20. Bogulskii I. О. The testing the scheme for solving two-dimensional problems //Modelling & Control, В., AMSE Press,France,1995.- Vol. 60, N 2 - P.21 - 28.

21. Анисимов С. А., Вогульский И. О. Численное решение задач динамики упругих тел. - Новосибирск: Изд-во НГУ, 1995. - 154 с.

22. Вогульский И. О. и др. Разработка компонентов системы сей-смоакустического мониторинга бурящейся скважины. - Красноярск, 1995 (отчет НИФТИ КГУ, гос.рег. N 16-94-41/1).

23. Вогульский И. О. О точности численного решения двумерной задачи динамической теории упругости // В сб.: Научные исследования на математическом факультете Красноярского госуниверситета. -

24. Вогульский И. О. и др. Математическое моделирование, информационные подходы в газодинамике больших скоростей и механике деформируемого твердого тела // В сб.: Актуальные проблемы информатики, прикладной математики и механики, Ч. II. Мат. моделирование. - Новосибирск-Красноярск, 1996. - С.157-197

25. Вогульский И. О. Алгоритмы высокой точности решения многомерных задач динамики твердых тел // В сб.: Математические модели и численные методы механики сплошных сред - Новосибирск, 1996. -С.158-159.

26. Вогульский И. О. Об одном численном алгоритме решения задач распространения сейсмических волн в вертикально-неоднородной среде // Геология и геофизика, 1997. - Т. 38, N 9. - С.1549-1560.

27. Вогульский И. О., Ветров С. Я., Шабанов А. В. Электромагнитные волны в неограниченных и конечных сверхрешетках // Оптика и спектроскопия, 1998. - Т. 84, вып. 5.

Деп. в ВИНИТИ, 1995. N1072. - С.72-83.

Вогульский Игорь Олегович

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

Подписано к печати Формат 60x84/16

Уч. изд. л. 2.0 Тираж 100 экз. Заказ N

Участок оперативной полиграфии ИВМ СО РАН 660036 Красноярск, Академгородок

й> iL ьгм

: . ' Y"ie ' y*ÎO С . i. ' Ц'-

В г - : 1льник упр^вл-:- Í:Í? БА*. - ,

.....Jdtukmï^^

Л / / . >

' \ / . /-./у

/ - ч

РОССИЙСКАЯ АКАДЕМИЯ НАУК Сибирское отделение Институт Вычислительного Моделирования

На правах рукописи

Вогульский Игорь Олегович

УДК 539.3

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ДИНАМИЧЕСКИХ ПРОЦЕССОВ В ТВЕРДЫХ ТЕЛАХ НА ОСНОВЕ СХЕМ ПОВЫШЕННОЙ ТОЧНОСТИ

05.13.18 - Теоретические основы математического моделирования, численные методы и комплексы программ

Диссертация на соискание ученой степени доктора физико-математических наук

Красноярск, 1998 г.

ОГЛАВЛЕНИЕ

Введение ............................................................- 4

Глава I. Схемы повышенной точности решения одномерных гиперболических систем..........< —.....26

§1.1 Одномерная гиперболическая система первого порядка.

Метод С. К. Годунова.........................................26

§1.2 Численное решение на основе нескольких аппроксимаций

неизвестных функций. Точность и монотонность .............33

§1.3 Монотонная схема второго порядка точности .................. 48

§1.4 Численное решение краевых задач для одномерных

систем гиперболических уравнений ...........................58

Глава II. Численное моделирование одномерных динамических процессов ...................................70

§2.1 Нестационарное деформирование пластинки с постоянными по

толщине деформациями сдвига (уравнения П.С. Тимошенко) 70 §2.2 Одномерные упругие задачи с осевой и сферической

симметрией....................................................78

§2.3 Моделирование процессов распространения электромагнитных

волн в слоистых диэлектриках ................................ 94

§2.4 Моделирование распространения плоских ударных волн

в анизотропной упругой среде ................................ 106

Глава III. Схемы рбшения двумерных задач на основе нескольких аппроксимаций

полиномами ................................................120

§3.1 Плоская и осесимметричная задача двумерной

динамической теории упругости .............................121

§3.2 Схемы решения плоской задачи на основе нескольких

аппроксимаций линейными полиномами .....................126

§3.3 Схема решения двумерной осесимметричной

задачи ........................................................140

§3.4 Алгоритм построения численных решений в областях,

состоящих из произвольных четырехугольников ............ .147

Глава IV. Итерационная процедура решения двумерйых задач ........................................ .163

§4.1 Корректировка решения путем последовательного

приближения .................................................164

§4.2 Симметричный вариант расщепления.........................170

§4.3 Двухэтапная процедура решения

осесимметричной задачи......................................182

§4.4 О точности решения двумерных задач ........................188

§4.5 Локальная аппроксимация решения полиномами

порядка выше первого........................................197

§4.6 Структура вычислительного алгоритма для

неоднородной области ........................................ 205

Глава V. Численное моделирование многомерных процессов распространения волн

в неоднородных упругих телах ......................... . 213

§5.1 Моделирование множественного ударного воздействия

жестких ударников на упругую плиту........................213

§5.2 Численное решение задачи распространения сейсмических

волн в вертикально-неоднородной среде......................241

§5.3 Определение некоторых механических и геометрических

характеристик слоисто-неоднородной упругой среды ........252

Заключение ........................................................267

Литература ........................................................269

ВВЕДЕНИЕ

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

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

Противоречие между попыткой повысить порядок аппроксимации

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

Однако существует ряд задач, для которых требования к точности решения оказываются достаточно жесткими. Одной из таких задач, к примеру, является прямая задача сейсмики, характерной особенностью которой можно считать достаточно малую зону действия источника возмущений по сравнению с масштабами всей расчетной области. Кроме того, сильно осциллирующий характер источника накладывает ограничения на шаг по времени, что, в свою очередь, приводит к необходимости достаточно мелкой дискретизации пространства. Для ее решения традиционно применяется полуаналитический метод, основанный на сведении исходной многомерной задачи к семейству одномерных краевых задач с последующим восстановлением решения [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], при расчете волновых процессов с большими градиентами все равно возникают ограничения на шаг по времени, вызванные соображениями не устойчивости, но точности, которые не позволяют выбирать его достаточно крупным. В то же время, использование неявных схем с шагами интегрирования, сравнимыми с допускаемыми явными схемеми, менее экономично ввиду дополнительных сложностей реализации первых.

Среди практически используемых неявных разностных схем наибольшее применение получили схемы расщепления [5