Содержание к диссертации
Введение
Глава 1. Алгоритм декомпозиции областей с использованием конечных интегральных преобразований Фурье-Бесселя 12
1.1. Постановка задачи для полярной системы координат 12
1.2. Теоретические аспекты построения решения 13
1.3. Аналитическое решение задачи для однородной среды 15
1.4. Численно-аналитическое решение для неоднородной среды. 16
1.5. Построение общей системы для декомпозиции областей и определение решения исходной задачи 18
1.6. Некоторые аспекты численных расчетов 20
1.7. Результаты численної о моделирования 21
1.8. Основные выводы 22
Глава 2. Алгоритмы моделирования волновых полей в вязкоунругих средах для пространственных осесимметричпых задач первого и второго порядка 28
2.1. Постановка задачи для системы уравнений первого порядка
в скоростях смещений и напряжениях 29
2.2. Алгоритм решения задачи для функций последействия вида 30
2.3. Решение полученной системы уравнений с функциями последействия
для обобщенной модели стандартного линейного твёрдого тела 36
2.4, Некоторые вычислительные аспекты реализации алгоритма 41
2.5. Результаты численною моделирования для задачи первого порядка 44
2.6. Постановка задачи для системы уравнений второю порядка
в смещениях 51
2.7. Алюритм решения задачи 52
2.8. Некоторые вычислительные аспекты реализации 59
2.9. Результаты численного моделирования для задачи второї о порядка 60
2.10. Основные выводы „ 65
Глава 3. Моделирование сейсмических волновых полей для неоднородных сред , 66
3.1. Постановка задачи 67
3.2. Алгоритм решения задачи 68
3.3. Некоторые вычислительные аспекты реализации алгоритма 75
3.4. Результаты численного моделирования 77
3.5. Основные выводы , 83
Заключение
- Теоретические аспекты построения решения
- Построение общей системы для декомпозиции областей и определение решения исходной задачи
- Алгоритм решения задачи для функций последействия вида
- Некоторые вычислительные аспекты реализации алгоритма
Введение к работе
Обзор вычислительных методов при решении задач моделирования сейсмических волновых полей. Актуальность и новизна предлагаемых алгоритмов.
В наешящее время существует большое количество различных методов решения прямых динамических задач сеіісмики Одним иї наиболее популярных, является л}чсвой метод, предложенный и работах [2], [3|, |7| и в последствие развитый в |Ш], J41], |50|. Лучевой метод с учетом нулевого члена лучевою ряда треб>еі шачигельно меньше вычислительных затрат, чем дручие методи. Кроме юю, с помощью лучевої о метода можно леї ко j честь вклад тех или иных сейсмических волн в формировании сложною волновою поля. Однако, с развитием высокоточной широкополосной сейсмологической аппаратуры, помнились факіьі рстиггра-щіи "нел>чевых" сейсмических волн, которые не описываются нулевым членом л} новою ряда, н для их вычислении необходимо учесть последующие члены ряда. В работах [G] и (51| такие "нелучевые" волны были обнаружены с помощью численно-аналитических методов моделирования сейсмических полей.
В начале GO-x годов, с появлением высокопроизводительных ЭВМ, в іеофи піке ста. і и использоваться численные методы расчета волновых полей, например, методы конечных разностей и конечных элементов. Возникло новое направление -вычислительная іеофизика. Первые результаты применения численных меюдов в сейсмике показали, что для получения приемлемых во точное і и результатов при численной дискретизации необходимо 15-20 точек па минимальною длину волны для расчет волновых полей на расстяние 80-100 длин волн ([34], [35], [37|, [53], ]34|). Такие требования оказались неприемлемыми, ддя существовавших в то время ЭВМ, по вычислительным затратам и оперативной памяти даже при решении дв)хмсрпых задач сейсмодопш и сейсморазведки.
В начале 70-х і одо в в Вычислительном Центре СО РАН (в последствие имени) і Вычислительной Математики и Математической Геофишки) стали развиваться численно-аналитические методы решения задач сейсмики, основанные на расщеплении дв) мерных и трехмерных задач на серию независимых одномерных с помощью интетрадьных преобразований по юризонтальным координатам и с последующим решением их конечно-разностным методом (см. [1], [5|, [19], [22], [33]). Такой подход не требовал большой оперативной памяти, обладал хорошей точностью и был реализован на отечественных ЭВМ типа БЭСМ-б .для раиичных моделей сред: для ради ал ьно-неоднородных - в работе [5], для анизотропных - в [18], [G8], ддя иязко)нр)піх сред - в [30]-[32]. Обобщение вышесказанною подхода
Лін лщ мерных и трехмерных моделей сред дано и работах [2flJ, [21 J, [23], [70)-(72].
Ни ікай точность и большие вычислительные затраты при расчете сейсмических волновых нолей на большие расстоянии стандартными численными методами дали толчок к развитию за рубежом конечно-разностных методов с высоким порядком аппроксимации (см. [3G], |42], |Ш], [63]). Кроме тою, нол)чилп рл шипі р пгеидоснеьтральныс меюдьі решения задач сейсмики (см. [СО] - [62|). Если в консчпо-ралюстных методах прои шодные аппроксимируются конечшьралюстным отношением в дискретных точках фишческою просгрансгва, то в псендоспек-іральньїх они определяются на основе разложения в ряды по базисным (функциям, которые бесконечное число раї дифференцируемы но всей расчетной обласні. Обычно, і: качестве таких базисных функций выбираются тршономеїрическис функции, либо полиномы Чебышева [18|, так как для них существуют 3(Jk|k>kihis-ні,те методы быстрою преобразования (БПФ). Важным достоинспюм нсевдоспек-тральных методов является высокая скоросіь сходимости (теоретически зкспонсн-пиальпая), если решение обладаеі достаточной степенью гладкости.
Аналої ичной точностью обладаюі и спектральные методы, предложенные ранее для решения задач іейсмики в работах [20], (21|, [23|, [71[, [72[. И отличие оі шевдосиектральных методов в спектральном подходе все вычисления проводятся в спектральном пространстве беї возвращения в каждый момент времени в физическое пространство. Основное время расчета в этом случае идеї на вычислении с} мм типа свертки, кошрые вычисляются с помощью БПФ или на спецпроцессорах. Сочетание спекіральноіо подхода но юризонтальным координаїам и высокоточными разностными схемами с переменным шаіом по вертикальной коордииаіе и явной разностной схемой но времени оказалось эффективным для решения МНОІ их динамических задач сейсмики (см. [23], [74]).
Все выше}казанные меюды расчета сейсмических волновых полей основаны пл использовании я иной разностной схемы шорого порядка по времени и условно нальшаюіся методами расчета во временной области. Дальнейшее развшие -этих методов связано с попытками исполыовать более высокий порядок аппроксимации производных по времени. Дейеинпелыю, спектральный и исевдоспектраль-ныи меюды дают очень высокий порядок аппроксимации нространсівешилх про-п шодных, а накопление ошибки при расчею волновых нолей на больших расстояниях и временах рас пространен и я происходит, преимущественно, за счет второю порядка аппроксимации временных прои зводных. Попытки использован. четве]ь іьій порядок аппроксимации но времени наїалкиваются на большие оіраниченнн на шаі дискретизации для тою, чтобы схема была устойчивой. В работе [811 предложено использовать спектральный подход на основе полиномов Чебышева для аппроксимации временною оператора. Метод обеспечивает высокою ючность, но
треб) er больших вычисли іельньїх іатрат
В работах [13], [55], |5б] предложен спектрально-разноегиый алтритм сведения задачи распространения сейсмических ноли в неоднородных средах к задаче Ко-иш для сисіемьі обыкновенных дифференциальных уравнений. Іі(Міікіе> s>^і меюд матричной декомпозиции, cut іема расщепляется на jV независимых обыкновенных дифференциальных уравнений, которые решаются аналшически. Поскольку, декомпошция (диаюнализацня) не ынисит от пространственною расположения неючника, который задается в правой час і и системы, то эта трудоемкая процедура проводится только один раз для любою положения исючника. Таким образом, отпадает необходимость решаїь задачу каждый раз при изменении положения исючника. Такой подход является весьма эффективным при моделировании сейсмических волновых полей в сейсморазведке методом общей іл>бинной точки (ОГТ), коїда необходимо решать задачу для фиксированною строения среды, но при различных (несколько сотен) положениях источника. Как разниіие іакою подхода, в данной диссертации предложен алюриїм сшивки аналитических и численно-авалиіическоїо решения для волновою уравнения н полярной системы координаї, рсз>льіаіьі моделирования, на основе данною алгоритма, были ои)бликованы в работах [10), \Щ
Вторая большая rpjnna методов расчеіа сейсмических волновых нолей основана на различных подходах в частотной обласні. Посіє применения преобразования Ф^рье по времени мы получаем уравнения или систему )равнений Гельмюльца. Эти уравнения необходимо решим» для набора временных частот, количеемю ко-юрых зависит от ширины спектра сишала в источнике, шем эти решения просуммировать, чтобы полечить импульсную сейсмоірамму. Наиболее известными меюдами решения и члсюміои области является метод Томсона-Хаскелла [17], [84], в дальнейшем развитый в рабоїе [25], а также "рефлективный" метод [10], [77]. В работе [32] предложен пол^аналигический алтритм, основанный на использовании преобразования Ф>рье-Бесселя по і ори юнталыюй координате и введением новых функций, которые сводят решение редуцированной задачи к решению сисісмьі уравнений Рикапи. В каждом однородном стос решение предсіавлясі-ся аналитически в виде комбинации экспонент только с отрицательным шаком. Просіые формулы пересчета волнового поля с одною стоя на другой ньіводяїся на основе іраничних ;словий.
Большое чисн) работ посвящено численным методам решения уравнений іина Гельмюльца, (см., например, [С5]-[б7]). Замеїим, что численные методы расчет сейсмических волновых полей и частотной области имеют свои особенности. По сравнению с методами расчета во временной области здесь пег оіранпчепип на шаі но времени, в то время как во временной области он определяеня млксп-
малыши перепадом скорої in распространения сейсмических ноли в выбранной модели среди. Кроме тою, при решении прямих динамических мдач но временной области для ису пру і их сред с* нос'.іедейсгниеч возникают проблочы, СЕЇЯ ІИШІЬІС с наличием и уравнениях иитеїралов і шіа свертки. Эти трудности иреододеваюкя iij it'M аппроксимации иитеїралов конечной суммой и введением дополнителЕшой переменной [38|, [79], что приводит к существенному увеличениіо вычислительных . В то время, как в чаепшюй области эти проблемы леї ко преодолеваются введением комплексных у пруї их нарамоіров.
Остановимся на особенностях численною решения прямых динамических задач it частотной области. После применения преобразования Фурье но времени и аппроксимации пространен пенных прои іводеіьіх, например, конечно-ралносшыч методом, либо методом конечных элементов задача сводиіся к системе аліебраи-ческих уравнений большой размерности Так как, матрица эюй системы зависит от временной частоты, то при решении системы матрицу нужно іраисфорчнро-ваіь для каждой частоты, что нредсіаіияеі чрезвычайно трудоемкую процедуру. Чтобы уменьшить вычислительные штраты, используют конечно-ра шостую аппроксимацию пространственных прои людных высокою порядка точности. И] >и ном, размерность матрицы уменьшается, но растет количество ненулешлх диаю-налой в самой матрице В последние юдЕЛ проіресс в этой обласіи сними с применением высокоточных компактных рашостных схем, с исполыовапиеч коюрьіх уменьшается размерность матриці>і, но не происходит увеличение ненулевых диа-юналей (см. [52], [80]).
Кардинальное решение вышесказанной проблемы связано с применением ин-нчралыюю преобразования Лаіерра но времени, имеет преобразования Фурі>е. Такой подход предложен в работах |12], [73]. В отличие от преобразования Ф)рі>е, применение иптеїральною нреобраюванин Лаіерра по времени с последующей дискретизацией пространственных переменных позволяет свести исходную іада-чу к решению аліебраической сие гем ел уравнений, в которой параметр разделении т приіуїсівуег только в правой част уравнений и имеет рекуррепшую мвиси-мосчь. В настоящее время вышеуказанный подход развит д,ія упру і их и юіропньїх [58], ]59] и анизотропных [09] неоднородных сред, а также для вязкоупруїих неоднородных сред |21|, [75[, [76].
В преддаїаечой диссергациоиной рабоїе рассматриваю!ся аліоригміл моделирования волновых нолей в вязкоупруїих неоднороднілх средах для сисіем уравнений первою и второю порядка но времени, основанные на применении иинчраль-ноео преобразования Лагерра.
Главная трудность в осуществлении численшлх методов в вязко)прутом моделировании - ирису іствис нніоірала сверіки в уравнении движения. В паеюящео
время с)іцеств)ют дна основных подхода при численном моделировании расіцкь (гранении сейсмическою волновою поля в иязкоупруюй среде. Первый подход основан на численном моделировании сейсмических полей во временной области. В этом случае, уравнение движения можеі биті, записано в дифференциальной
форме, ВВОДЯ ДОПОЛПИТСЛЬПЫи ПСрСМСННЫС В ИСХОДПЫе форМ)ЛЫ, КОЮрые ПОЗВО-
ляюі исключить ш уравнений интеїральньїе свертки [39], [43], [15]. Полученные, іаким образом, динамические jравнения нязкоунрутости моїут бып. решены, используя конечно-разностный мегод аппроксимации производных (см., например, [44], ]85]), Также, использ>ется исевдоспектралышй метод, іде наряду с дискретизацией по пространству применяется полиномиальная интерполяция опюсннмь-ио времени [39], [82]. Псеидоспектрлльные методы обеспечивают более высок) ю іочносіь решения, потребуют больших вычислительных чаїраі.
Второй подход связан с численным моделированием сейсмических полей в вя t-коунруюй среде в частотной области. Эюг спектральный подход основан на применении преобразования Ф)рье но времени к динамическим уравнениям вязко-)ир)10сги. Чтобы преодолеть трудность численною нредсіавления иннчралон cue-pi кн, можно преобразовать )равнения в частотную область и моделирован, результирующие )равнения типа )равнения Гельмюльца. Это позволяет пред-(іавиїь часюгно-занисимьіе парамеїрьі поілощения в виде комплексных ) пру і их параметров [78] Использование конечно-разностноіо меюда и метода конечных цементов, для полученных таким образом уравнений распространения ynpjioii волны в частотной области, было впервые предложено в работе [С5] и позже раз-рабоыиы в работах [GG[, [67[. Пол)аналиіический меюд решения в часюїной области ,гдя вязко) пру і ой слоистой среды был предложен и работе [32|, іде решение исходной задачи сводится к решению сисіемьі уравнений Рикатти. После чет, в каждом ) пру і ом слое, решение задачи может быть представлено аналишческн в виде комбинации убывающих экспоненциальных функций. Общее решение сіро-иіся по с[)0рм)ла\1 пересчета волновою поля от нижней к верхней і ранние на основе заданных іраничньїх условий. В случае произвольно неоднородных сред, после тою, как полненные уравнения в частотной области дискрет шронаны по про(тран( і ценным координатам, общее решение задачи представляется и виде систем ) равнений чрезвычайно большою порядка для каждою значения частоты. Численное решение такой іадачи іребуеі больших вычислительных заіраі. Чю-бы облегчить решение задачи и )менынить обьём вычислений, в работе [80] предложена новая конечно-разное і пая схема решения и масло гной обласні, котрая (нпована на операторах вращения.
Чн> касается предлагаемых в диссертации алгоритмов моделирования волновых полей в вязко)пр)інх неоднородных средах, то их можно рассматривать как
ана-юг спектральною моделирования, іде вместо временной часини и мы имеем номер т - енчіень многочлєноеї Лаіерра. Главное отличии данною менш от спеьлралЕ>ноЕО моделирования заключается в том, что после применения шне-іральноіо преобразования Лаіерра к временным проичиодныч и для шпсі радон свертки но времени с последующей диекроїн зацией но пространственным коордн-наїам, мі>[ подjчаем систему аліебраических уравнений с маїриней, независимой ог параметра т. Так как, юдько правая часть системы имеет рск}ррепін}іо зависимость or параметра т. Таким образом, возможно исиольчоваїьбьн іріле методіл решения систем линейных аліебраических уравнений с большим числом правых сторон, например, на основе разложений по методу Холецкою В этом случае, численные операции по преобразованию матрицы осущсствляюіся только только одни раї и отличии оі сиокіралішою подхода, когда необходимо прообразовать матриц} для каждой частотной іарчоннки w, что требует больших пычиедшоль-ных затрат. Сравнивая пнтсчральные преобразовании Ф}рье и Лаіерра, можно сказан., чю применение преобразования Лаіерра, в этом случае, позволяет естественно сократить вычислительные іаграгьі. Дія осуществления такою подхода к решению динамических задач вязко}пр}іости потребо на. і ось доказать теорем} для преобразования Лаіерра иіпеїральной свертки. Эта теорема подобна іеореме преобразования Ф>рьс питої радон сверіки. Её использование даёт возчожносп. рассматривать саміле общие связи межд} тензорами напряжения и деформации, выраженными произвольными релаксационными функциями в интегральных соотношениях Больцмана При этом не требуется введение доиолншелыплх переменных, коюріле приводят к появлению дополнительных уравнений и и'м сачілч }величивают размерность решаемой сие і омы. Предлаїаомая методика может быть леї ко обьединена с рядом методов лея решения полученной, после преобразования Лаіерра, дифференциальной системы по пространственным координаїам, включая использование конечных рачносюй или спектральные меіоділ (Фурье, Ф}рьс-Бесселя или преобразование Лежандра). D предлаїаемой диссерьщии такой иод-ход к решению динамических задач вязкоупруюсти обсуждаеіея на примере решения просіранстненшлх осеснммеїричньїх задач в цилиндрических координаїах Лія вертикально-неоднородных сред и ;ыя 2.5D .моделей сред is Декартовых координатах с релаксационными ф}нкциями для стандартною линейною тердою і ела.
Содержание диссертации состоит из трех глав, содержащих 15 нараірафов, заключения и двух приложений,
В первой їлаве рассматривается алюритч решения двухмерно-неоднородной }И[)}Юй задачи с помощью декомпозиции областей на основе шшчральных преобразований но пространственным координатам. Дія просюты, иредллі аен я рас-
смоірегь описываемый алюри їм декомпозиции областей иа примере решения дія волновою уравнения в полярной снсіеме координат, используя разбиение просі рапсі пенной области на ірн ;часіка. В качестве иллюстрации воїможносін использования алюригмадля более сложно-построенных моделей сред, прпиодяїся рсззльтаты моделирования волновою поля для рад и ал ы і о-н сод пород 11 ой мололи Земли и случае построения решения при сшивки нескольких численных и аналитически решений.
Но [порой їлане рассматрнваюіся дна алюритма моделироиания волновых полей іі вязко>нр}іих средах для решения пространственных осесиммегричных за-дач первого порядка в скоросіях смешений и напряжениях и для ) раїшений второ-іо порядка в смещениях. Среда -задастся изотропной неріикально-неоднородной При эгом нолаїаетея, что механизм последействия задан и виде митральных соотношений Больпмана для проишольїшх функций последействия. Главная оі-лнчительная особенность предлаыемых алгоритмов заключается в использовании интсчралыюю преобразования Лаіерра по временной координате Построение аліориімоіз основывается на комплексировании интеїральньїх преобразова-ний и конечно-разностшлх методов дли сведения исходных интеїрально- дифференциальных систем к сисіемам линейных аліебраичоских j равнений, решения коюрых мої уг бы п. найдены наиболее эффективными известными численными мої одами (например, тина LU— разложения). Приводятся результаты численною моделирования для вязко} npj і их сред с ф) нкниями последействия для стандартною линейною твердої о тела с несколькими релаксационными механизмами. В заключи іельном пункте данной і лавы делаются выводы об эффективности применения предлагаемых алюритмов, Рассматриваются отличительные особенности и преимущества в сравнении с друїнмн и шести ими методами решении подобных задач.
В іреіі>ей і лаве рассматривается алгоритм моделирования волновых нолей дли 2.5D неоднородных вязкоупр>іих сред Построение предллілемоіо алюршма основывается на комплексиронании иніоіральною преобразования Лаіерра по временной координате и численной конечно-разностной схемы аппроксимации npotn-водных но пространственным координатам для сведения исходной интеїрально-дифференциальной системы к системам линейных алісбраических ураиненнй. Псь добное сведение задачи к хорошо обусловленной системе аліебраических уравнений с множеством правых частей позволяет исполі>зовать быстрые алтрнтмы решения на основе итерационных меюды тпа сопряженных ірадиеніОЕЗ, сходя-ншеся к решению задачи всею за несколько итерации. На данном этапе алюри їм был эффективно распараллелен В частности, была реализована распарал-леленная версия меюда сопряженных іраднентов. На уровне входных данных,
при задании модели среды, это равносильно декомпозиции исходной обласні па множесшо нодобласіей, равных количеству процессоров. Это дасі возможность распределения памяїн, как при задании входных параметром модели, іак и при дальнейшей численной реализации адюрніма в подобластях. В заключительных н>нкіах ідави приводятся результаты численною моделировании волновых нолей и вя!ьодпр)1их средах с ладанными функциями последействия и делаю і ся выводы об -зффектиипосги использования предложенною алюритма.
В дв>х последующих приложениях приводятся некоторые известные теореш-ческис постулаты и формулы, используемые для построения предлагаемых в диссертации алюритмов. В приложении 1 описываются некоторые евойсіва и приводятся основные формулы для интегральною преобразования Лаіерра. Расемаї-риваюіся і рафики (фикций Лаіерра и зависимости от разных параметров и их влияние на спектр коэффициентов разложения но полиномам Лаіерра заданных ф)нкциґі. В приложении 2 показывается принцип построения и вывода основных фор\і}л іеории вяіко)П|>ііосіи для законов помещения и дисперсии ( заданной функцией последействия.
Теоретические аспекты построения решения
Колнчесіио требуемых корней эшх уравнений для конкреіною шачения п = 0,1,..., Лг, иптльз}емы\нри вычислении функций 4(0 6 (і), » , (/), wf t) пофіин м}лам (1.24), (1.25) и построения решений Sn(r,t),\Vn(r,t) задач (1.7), (1.9), определяется исходи in ширины спектра иолноиою поля в исючнике, задаваемою ф}нкцней f(t). На этом основании, можно уменьшить количество суммируемых
Главными факторами, влияющими на поірешносгь определения функций ( ?, ((), Qi(t)iQn(t) Qi(0 а слеДвате-Ы10 и иа конечное решение, являются точноеп аппроксимации прои шодных но г при решении задачи (1.19) и выбираема» ширина перекрытия Аг решений Stl(r,t),]V„(r,t),Pn{r,t). Как показа.! авали Ї численных расчетов, в ел} чае аппроксимации нрои зводных со вторым порядком точносіи, оп-інм.иьиая Лг составляет четверть нросірансгвенной длины волшд ,(ля волновою поли, сооівеіспз}юіцем} mtn{v\,Vi}.
Смотрим результаты тесюиых расчетов на примере модели среды изображенной на рисунке 1.1. Значения скорости упруїих колебаний vp{r, p) для данной модели среды, были заданы как: Vp{w)= vi , г п, 0 180 Vi , п г г2, 0 ip tpi U f?i f 180 , vz, г r2, 0 9 180 іде vi = 1.5 км/сек, у2 = 3 км/сек, гі = 20 км, Ті = 40 км, і = 60, 2 = 120. /Ідя решения задачи (1.1)-(1.3), следуя описанному алюритму, данная моделі» была разбита на 3 }частка но координате г: [0,ri] , \т\ - Дг, г2 + Дг] и [г А]. То-іда на интервалах [0, rt] и [г2, А] решения моїут быть записані,! и явном виде но формулам (1.14) и (1.18)соответственно. Л на интервале [г1-Дг, rj+Дг] использовалось численное решение в виде (1.21). В результате конечное решение сіроніся І\Л основе декомпозиции трех областей, иснольЗ}Я прельщаемый адюритм сшивки решений (Ы4), (1.18) н (1.21). Ширина области перекрытия решений Дг была іалама на интернале- 5-тн узлов рлзнсхтной аппроксимации протводных для численною решения (1.21). Граничное значение г = Л было выбрано іак, чтобы на рассматриваемом временном интервале волна не достигла данной іранинм. (2т/о( -(о))
Результат расчетов волновою ноля дли мчанной модели среды продсілвле-ны на рисунках 1.2 - 1.4. На представленных рисунках изображены мпювенные (нимки волновою ноля для трех фиксированных момен іон времени. Сплошной линией показана іраииіщ раздела среды. Местоположение исючника обозначено темным іреуюльнпком. Временной сні нал к источнике задавался и виде импульса Пу зырёиа /(О = ехр V bin{2wf0[t - tQ)), (1.27) іде , = 1, /о = 1 Гц, tQ=\,5 сек.
Анализ порченных реїулілаюв показы пае і, что ноірешиость расчеюв за счеі сшивки аналитических и численною решений не превышает поіреншосіи построения численною решении (1.21),
Описаними it данной їлане алюритм может быть испольїонлн и для более сложно-построенных сред. Как, например, для мною-слоистой среды, исполмуя сшивку аналитических и нескольких численных решений. Ма рисунках 1.5- 1.10 предсмвлены результаты расчетов ,ия радиально-нооднородной модели Земли. Гра(])ик задаваемых значений скоростей в зависимости от кабины изображен па рисунке L.5. На рисунках I.G - 1.10 представлен ЕЛ мпювенные снимки шиповою поля для 5-ти фиксированных моменшв кремени. Источник расположен на по-верхінхіи Земли при у = 90s. Временной сигнал в источнике задавшіся в киде импульса Пу зырё ва по формуле (1.27) для част ОТЕЛ /0 = 0.1 Гц и to = 15 секунд
В заключении следует отметить, чіа предлагаемый алюритм декомпозиции областей может быть применен и ЛЕЯ сред с более сложной структурой скороеI-иою разреза, если ичеюісн участки с посюянной скороспао. Данный алюритм позволяет существенно сокраппь время вычиенпелыюю процесса и уменьши і ь требуемый обьем машиной онераіивной памяти. Основной объем вычислений для данною алюритма приходится на построение численною решения \.п\ двухмерно-неоднородной среды, іде используется разложение матрицы еисіечм на собспзен-ные вектора и собственные числа. Как известно, при численных расчеіах увеличение порядка матрицы в Дг раз приводит к увеличению времени вычислений порядка Дг1 раз Использование аналитическою решения на ишерналах [0,п] , [rj,/?j позколяег, при необходимое пі, увеличивать ірапицьі зтих обласіей без су-щеспзениою увеличения объема н времени вычислений.
Данная посвящена новым подходам к численному решению динамических задач генсмикн для неоднородных нязко)нр}іих сред, Актуальность решения подобных мдач вызвана необходимое і і.ю моделирования сейсмических нолей eyie-точ законов ноі лощения и дисперсии в реальных средах. Такая модемі, может бы и. описана в рачках общей теории линейной вязкоупруюсти. Основная і рудное 11. при решении подобных ілдач связана с наличием интеїрала свертки по времени вправлениях последействия Больцмана, Эта проблема, водном ел;чае, обходится с помощью введения дополнительных переменных, позволяющих исключшь ипіеіральн к свертку за счёт сведения исходных уравнений к системе дифференциальных уравнений большей рлзмерносіи но временной координате [39], [f 1], (83. Полненные, таким образом, динамические уравнения внзко нруюсіи мої} г б[ іть решены, используя конечно-ра зносі ный метод аппроксимации производных [11], [85] или исиольч я шевдоспектральный метод, где применяется полиномиальная шперноляция относшелыю времени [39, [82.
Построение общей системы для декомпозиции областей и определение решения исходной задачи
В аналитических преобразованиях Ф рье-Бесселя (2.4)-(2 7) и Лаіерра (2.13) при определении значений функции по их спектр) используются формулы обращения и пиле с)мм с бесконечным пределом. При численной реализации необходимым условием являемся определение требуемою количества членов с ммпр;емою ряда дли построения решении с заданной точностью. Так количество іармоник /.-„ и формулах (2 4)-(2.7) зависит оі просірансгвенной длины полны и моделир)е-мой среде и значения апертуры носі ілнаїзлииаемою поля, задаваемою конечными пределами ипнчральною преобразования. Кроме тою скороеіь сходимости суммируемою ряда заиипп оі їладкости функций моделируемою волновою поля. Количество іармоннк т по Лаіерру необходимых для определения ф;нкцнн по формулам (2.13) зависит оі злдлваемою cm нала в источнике /(() и значения временною интервала восстанавливаемою волновою ноля. Определи ь их коли-чесшо можно следующем образом. В начале находим коэффициенты разложения для заданной функции /(/) но форм;ле ос fm = Jf(t)(htyb«m(ht)d(ht), (2.31) о іде / - функции Лаіерра Далее можно вычислить коэффициенты разложения по Лаіеррі дія функции f(t + TQ), используя аналитическую формул} решения iiKjfiit iccKOfi падачи для плоской волны (12]: k Jfc — I Д = /т -м(ВД - / -m-l№) (2 35) Здесь L\ - орюнормированные полиномы Лаіерра нулевою порядка.
Анллизируя сходимость, вычисленных іаким образом, коэффициентов Д, можно найти треб е\юе количество іармоник дія восстановления волновою ПОЛЯ В заданный момент времени То- Функциональный вид козффициенюв Д = f(k) для разных Т0 представлен в приложении 1. Из анализа которою видно, чю чем больше время ветвления сигнала, і ем более сдвиЕіут вправо ею спектр по Лаіерру. Следовательно, ею более старшие іармоиики отнечлюг за восстановление сиінала для больших времен ею вступлений. Поэтому найденное КОЛИЧССІВО іреб}ЄМЬ1\ ілрмонпк лія восстановления волны с большим временем вступления ІІОШОДЯЄІ определить все волновое поле до данною момента времени. Подробнее зависимость характера этою спектра от параметров а и к, а также ею сходимосіь, рассмоірен в рабше (12j.
ІІроаналгіир)СМ влияние нарамеїров преобразования Лаіерра а и к на точності) вычислений и построение конечної о решения. Дія ЗІОІ О, В начале, рассмотрим адюрнтм численного решения. Как было показано шлше, решение исходной задачи (2.1)-(2.3) сводится к решению системы линейных алюбраическнх jравнений (2 33), В настоящее время существ}ст достаточно мною численных адюрит-мон решения подобных задач. Наиболее подходящими дли решения данной задачи являются алюритмы, позволяющие находи І Ь решения системы аліебраическнх сравнений для фиксиронанной матрицы сисіемьі и большою набора правых частей с наименьшими вычислительными шраіачи. Поэтому дія решения данных пкчем испольювался меюд основанный на представлении матрицы системы А и виде произведения верхней и нижней іре}іольньіх матриц, т.е. A = L (У, где L - нижняя тре юльная матрица с единичной диаюналью. Тої да решение ешчемы линейных уравнений Ах = b дія жданных правых частей Ь сводится к последовл-іельному решению уравнений Uj = Ь, Ux = y.
Их решение осуществляется мсюдомобраиюй итерации. Данный аліоршм но иль ляет }читывать ленточную структур} маїрицьі А при ее разложении на нижнюю и верхнюю тре}[Ольные матрицы, которые н результате имеют также ленточную (Чріміру. При лом ношожно при одном таком разложении, находи 11 решение системы дія миоіих правых часіей. Эю существенно уменьшает время раечёюв, ьтк как основное количество нычислин льных операций при решении сисіемьі дія заданной правой части приходится, именно, на разложение матрицы системы на тре юльные.
Как ншесіно, нее численные алюритмы для работы с матрицами оеновыва-іоіся па достаточно хорошей обусловленности этих матриц. При этом дія увеличения точности вычислений требуется уменьшение степени обусловленности \нллиі реі)ліітаюи численною моделирования волновых нолей и неоднородных средах показывает, что большое нлиянио на сіепень обусловленности матриц систем, пол} чаемых при решении аналої ичных задач, оказываю і шачения фіпиче-(кпх харакіерисгик в моделируемой среде. Особенно для сред имеющих юнкие слои и слои с резко-контрастными іраницами раздела, связанньїе сбол[ шим перепадом шаченин физических характерисіик этих слоев. Например, ірлнніш типа воід}х-вода-земля. В нредлаїаемом алюригме степень обусловленное! и матрицы дія системы }равнений (2.33), в большей степени, зависит от шачения иарамеїра її преобразования Лаіерра. Это достигается в результате диаіопальною преобладания шаченин /і/2 располаїающихея на їлавной диаюнали матрицы системы. С помощью выбора которого, можно уменьшать обусловленность маїрнц. Как показываю! ре))льгаты численных расчеши даже для сильно-контрастных сред, сіепень обусловленности таких матриц не превышает одною порядка и при увеличении шачения на[)аметра h значение об)словленности стремится к единице. Выбор (ooi веч с пукнцею шачения h осуществляется на основе оптимизации схо-ДПМОСІИ с;ммнр емых рядов дія носі роения решения.
Алгоритм решения задачи для функций последействия вида
Как видно, и J вышеприведённых рассуждений, решение исходной ыдачи своди і-ся к решению сисіемьі линейных аліебраических уравнений. Анализ полеченной, іаким образом, системы (2.66) покашвлет, что матрица сисіеміл А имеет 7-диаі опальную лен гочп} ю cip Kijpj и является симметричной и положи і ел ьно-опредеденной. При этом на главной диаі опали матрицы раснолаїаются компоненты входящие в уравнения системы, как слаїасмью имеющие в качестве сомножп-іедя параметр h2/4 (парамеїр нреобраїования но Лаіерру). Выбор значения h позволяет влиять на обусловленность матрицы системы (2.66). Так как, матрица системы является симметричной, то в отличии от задачи полученной и предыдущем пункте, здесь можно для решения системы аліебраических уравнении вое-нодьюваться методом Холецкою, основанном па представлении матрицы системы .4 в виде произведения тре юлыюй матрицы на её транспонированною A = L LT. Ныиолнив один paj данное разложение, можно затем найти решение дли всех правых частей.
Таким образом, можно определим, решение ддн N независимых систем уравнений (2 60) для каждого кп корня уравнения (2.46). Количеспю суммируемых іармоппк ки в формулах (2.11) ілнисит от пространственной длины волны в моде, in pje мой среде и значения anepijpbi восстанавливаемого ноля, задаваемою ьо-нечнымн предсіами интеїральною преобразования (2.45) Кроме нно скорость сходимости суммирземою ряда зависит от 1.1 ад кости функций моделируемою вол-повою поля. Количество ілрмоиик т по Лаіерру необходимых для определения функций по формулам (2.51) записні от задаваемою сиі нала в неючнике f(t) и значения временного интервала восстанавливаемою волновою поля. Количество членов в сверточпых суммах для определения правых частей системы (2.66) зависит о г сходимости коэффициенюв Лаіерра функций последействия, определяемых по формуле (2 53). Анализ тестовых расчетов показывает, что сходимость коэффициентов для заданною типа фзнкций последействия зависит от требуемой частотною диапазона функции добротности Qp(u), Qb(v), и котором они задаются постоянными для заданной часюпл пинала в источнике. Для моноча-СЮІНОЮ ешнала и при выборе оптимальною значения параметра h преобраю-пання JIaieppa требуемое количество іармоник и свёрючных суммах сосіавлнеї порядка 20-30.
Лплли! точности шлчислений представленного алгоритма проводился на примере сравнения расчетов между численным и аналитическим решениями в однородной изотропной вязкоупруон среде. Физические параметры тоеюной модели среды были заданы следующие плотность р = 1 і/см3, значения скорое і и продольной волны CP(UJ = 0) = 1.5 км/сек и скорое і и поперечной волны cs(u = 0) = 1.0 км/сек для нулевой частоты. Функции последствия в среде были описаны двумя релаксационными механизмами (Lp = Ls = 2). Значения нремени релаксации на-пря/кения а и деформации г в функциях последействия $\(t) и Q2(t) (2.57) были їадаїпл следующими: 1) лій Р-водны те1 = 0.4107 сек, т2 = 0.06G61 сек, таХ = 0.401 сек, та2 = 0.00501 гек, 2) для 5-волны г,і = 0.4160 сек, те2 = 0.00529 сек, таХ - 0.1012 сок, тоі = 0.003 (ЄК Времена релаксаций и моделируемой среде подбирались так, чтобы задать постоянными значения добротности Qp = СО ;іля продольной волны и Q , = 40 дли поперечной полны в частоіном диапаюно задаваемою сиінала н иеіочнике. Временной сні нал в источнике задавался как: ып(27г/о( - о)), (2.G7) іде-7=1, /0 = 1Гц, t0 =1.5 сек Графики функций добротности Qp{u), Qsiyi) и амплитудный спектр пинала нредстав-іешл на рисунке 2.1. Значения фазовілх скоростей Ср(и ) и cs(u) изображены на рисунке 2.2. Их значения вычислялись но формулам (2 39) и (2.10).
Рисунки 2.G и 2.7 представляют сравнение между численным и аналитическим решениями ,ия заданной іссювой модели среды. Лнадшнческое решении для вязко}npjiой среды было подучено, применяя принцип соответствия к аналитическом) спектральному решению j пруг ой однородной задачи [39]. На данных рисунках изображены расчётные сейсмотрассы для всріикалішон компоненты смешения uz, которые рсчистрируются на расстояниях источник-приёмник лія J"i = V ( i - -о)2 + (п - го)2 = 10 км (рис. 2.G) и х2 = \J{z2 - z0)2 + (r2 - г0)2 = км (рис. 2.7), Пространственные координаты источника (г0,с()) и приёмников (Пі іїі ( 21 2) были заданы идентичными ,ия аналитическою и численного решений задачи. Дли данной задачи (г0,с0) = (0,0) и (п, ) = (С,8), (г2,г2) = (12,1С).
На лсиой диаірамме рис нка 2.G показано сраішение между лналишческнч и численным решениями и шчке (п. і)» іде численное решение было полечено н результате суммирования 100 іармоник преобразования Лаіерра На правой дпаїрачме данною рисунка нредсіанлена сейсмотрасса для численною решения нол;чепною суммированием 150 іармоник Лаіерра. На рисенке 2 7, ючно так АС, слева показано сраішение между аналитическим и численным решениями і; точке (rj, гД где численное решение было получено суммированием 150 іарчо-ник, а справа - суммированием 300 іармоник Лаіерра, Значения амплт д были пронормированы по максимальному значению амплитуды сейсмотрассы и ючке (г],Сі). Вычисления были сделаны для параметров функций Лаіерра о = 8 и h = 20 Количество членов в снёршчных суммах было задано равным 25. Как можно заметить HJ представленных рисунков, чем больше интервал времени на ко юром требеегся полечить решение с нужной точностью, тем больше требе с юн количество суммируемых іармоник в ряде Лаіерра. Анализ тесюиых вычислений показывает, что, увеличивая число іармоник до 300, можно полечить решение с ючпопыо до 10_t для t 25 секунд. Следовательно, задавай нужное количество іармоник Лаіерра д. і я носеїаноізлепин волновою поля для наибольшею времени, возможно полечить решение на всем заданном интервале нремени с іребеемой точностью.
Некоторые вычислительные аспекты реализации алгоритма
Численные расчеты для задачи (3.1)-(3.3) проводились для модели среды соеюя-шей іпупруюю слоя и ди}х тонких иоілоіиающих слоев, лежащих на иттропном jupyioM пол}пространстве и разделенных между собой структурой інпа ньісіун. Моделі) среды и плоскости XZ и юбражена на рисунке 3 1 (по оси У среда однородна) Фіпические характеристики для сооїиеістнующих слоев среды были аданы с, іе, о ющнмн. 1) Д.ІЯ верхнею унр}іою слоя - р = 1 і /см3, Ср = 1.5 км/сек, cs = 1.2 км/сек. 2) Для левою и правого вязко}пр}іих слоев р = 1.2 г/см3, ср(и = 0) = 1.7 ьм/сек, (а(ш = 0) = 1 4 км/сек Поглощение для продольных и поперечных ноли в -них слоях было іадано соошентвенно: Й ЛІЯ девою слоя Qp = 15 и Qs = К); 0) для праною слоя Qp = 70 и Qs = 50. 3) Для нижнею у пру юі о пол} просі ранета-/? =1.5 г/см3, Ср = 2 км/сек, rs = 1.6 км/сек.
Толщина верхнею слоя-0.75 км. Толщина ноілошаюших слоен- 0.25 км Функции последействия для поілоіцакнцих слоев описывались тремя релаксационными механн імами и сооївеїсгву юіцие значения времен релаксаций т и та подбирались так, чтобы в частотном диапазоне моделируемою сшнала в источнике задать нсь стоянии\ш значения добротности Qp и Qs для продольных и поперечных волн (ооівеїственпо. Форм}ЛЫ л їн вычисления частотной заиисимосш (функций доб-роїносіи Qp(u), Qs(u) и фазовых скоросіей rp(w), ca(w) для піюдольньїх и поперечных волн приведены в ііреділ,гуіцей главе. Волновое ноле моделировалось от ючечпою исючника типа центра давления с координатами х0 = 1.5 км, yQ = 1 км, ZQ = 0.25 км. Временной сипыл в исючнике задавался в виде импульса Пушрёна с основной частотой /0 = 30 Гц.
Расчеты волновою ноля для заданной модели среды проводились на основе предложенного алюритма распараллеливания па многопроцессорном комплексе МВС-1ШЮ/М. При проведении вычислений было чадано 120 ілрмоннк Фурье-преобразонания по координате Y и сетка размерностью 600 узлов по коордишие Лг и 100 }Тлов но Z с ніаюм Ах = Az = 50 метров. Время расчетов ,іля данной размерности чадачи сое ганило 32 часа машинной) времени на 20 процессорах типа DEC Alpha 21264 (833 МГц/1 Мб SLC). Результаты расчётов предетавлены на pucjHKax 3,2 - 3,8.
На рисунках 3.2-3 4 представлены множенные снимки шиповою ноля для вершкальной компоненты скорости смещений uz(x, у, z) дчя трех фиксированных моментов времени / = 0.2, t = 0.5 и t = 0.9 секунд соответственно Волновое поле изображено « виде цветной ідмми соответствующих значений uz компоненты на трех перпендикулярных плоское і ях X — х0 — 1.5 км, У = г/о = 1 км и Z = 1 км
На рисунках 3.5 и 36 представлены мпювенные снимки волновою ноля дія вертикальной компоненты скорости смещений иг в момент премсни г = 0.9 сек в плоскости у = у0 = 1 км (рис, 3 5) и в плоскости z = 1 км (рис, 3G). Границы слоев заданной модели среды на зіих рисунках показаны сплошной линией.
На рис;пках 3.7 и 38 представлены расчётные сейсмотрассы для г/.(() - компоненты (рис. 3.7) и ux(t) - компонешы (рис. 3.8) рсчистрир емые па приемниках расположенных на линии пересечения іррх перпендикулярных плоскостей изображенных на рисунках 3.2 - 3,4. В плоскости XZ, изображённой на рисунке 3.5, линия расположения приёмников проходит но нижней іранице ноілоіпаїоіпих слоев z = 1 км. Шаг по коордннаїе А между приёмниками равен 100 чеірлм, МЮ\У-динага первою приёмника х\ = 0. рассмотрения представленных рис нкон можно заметить различие дейп ния заданных функций последействия на характер распространения волновою ноля в слоях с разным поглощением. Данные различия в основном проявляются в виде раз ною амплитудпоіо чат уха и ия волн в поглощающих слоях. Так, например, явно заметно большее ослабление амплитуд персотражённых волн в левом вязко-}пр іом слое, іде значения поглощения дія продольных и поперечных ноли было задано большими, чем в нравом. Явление разности дисперсии скорости и помещающих слоях менее заметно. Это связано с небольшим раишчием фазовых скоростей дія моночастотного сшнала, которое проявляется в виде небольшою сдвим времен вступлений соогкектв}тощих типов воли в левом и нравом поілощающих
Эффекшвность использования преобразования Лаіерра ддя задач внзко}нр}-іосіи }жеоімечалась в предыдущей і лаве. Можно добавить, что такой подход может быть также иснолілован и при решении задач для трёхмерно-неоднородных сред При этом решение получаемой после преобразования Лаі ерра дш] ференцн-альной задачи но пространственным координатам может быть осущесівлено, как и в описанной в этой їлане способом комнлексирования аналитических и конечно-разностным методом, так и с использованием трёхмерных нросірлік і венных сеток ,ия разностной аппроксимации [44], [85. В этом случае, процесс раснараллелива-ния решения производится на «ЛИР нахождения решения системы лліебраических j равнение. Пспользземый в преддаїасмом алюри гме метод сопряжённых і радиен-тон но зволяег по леї ко сделать Такой подход распараллеливания был опробован на данном ллюршме при моделировании волновых полей для больших моделей (ред. Так как при решение іаких задач іребзегея большой объём оперативной намят. А при данном типе распараллеливания, на каждом процессоре ирои зводшея построение вектора решения только для ограниченною пространственною \часі-ка среды и следоваїельпо только для определённою }частка матрицы решаемой сисгемы. При этом требуется обмен межґгу процессорами найденными значениями вектора решения только в двух соседних пространственных злах разносіноіі сегки, так как от параметра Лаіерра матрица решаемой сие і ем ы не записи і, а правая часть }рашіений насчитывается независимо для решаемой части сисіемьі на данном процессоре
