авторефераты диссертаций БЕСПЛАТНАЯ БИБЛИОТЕКА РОССИИ

КОНФЕРЕНЦИИ, КНИГИ, ПОСОБИЯ, НАУЧНЫЕ ИЗДАНИЯ

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


Pages:   || 2 | 3 | 4 | 5 |   ...   | 15 |
-- [ Страница 1 ] --

Посвящается пионерам освоения космоса

Предисловие

Фактически в настоящее время закладываются основы решения фундамен-

тальных проблем, связанных с

Землей как средой обитания человечества в XXI

веке. Они приобретают все б льшую актуальность и международное значение,

о

так как от их решения зависят и проблема выживания человечества и ответ

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

Еще в 1938 году В. И. Вернадский опубликовал очерк «Научная мысль как планетное явление», в котором пророчески писал:

«Исторический процесс на наших глазах коренным образом меняется.

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

Меняется человеком физически и химически воздушная оболочка суши, все ее природные воды... В будущем нам рисуются, как возможные, сказочные мечтания: человек стремится выйти за пределы своей планеты в космическое пространство. И, вероятно, выйдет... »

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

В теоретических прикладных исследованиях стал применяться термин «глобальная система» (ГС). Необходимы анализ и синтез знаний о развитии планетарной цивилизации в рамках теоретико-прикладных концепций об орга низации глобальной действительности. Материальная основа ГС базируется на двух фундаментальных уровнях: атомно-молекулярном, формы связи – – электромагнитные силы и поля;

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

«Своеобразным, единственным в своем роде, отличным и неповто ряемым в других небесных телах представляется нам лик Земли – ее – изображение в космосе, вырисовывающееся извне, со стороны, из дали бесконечных небесных пространств. В лике Земли выявляется поверхность нашей планеты, ее биосфера, ее наружная область, ограничивающая ее от космической среды. Лик Земли становится видным благодаря проникающим в него световым излучениям небесных светил, главным образом Солнца. Мы едва начинаем сознавать их разнообразие, понимать отрывочность и неполноту наших представлений об окружающем и проникающем нас в биосфере мире излучений, об их основном, с трудом постижимом уму, привыкшему к иным картинам мироздания, значении в окружающих нас процессах...» («Биосфера как область превращений космической энергии» из очерка В. И. Вернадского «Биосфера в космосе», 1926 г.) Солнечное излучение в коротковолновом диапазоне спектра длин волн – – один из неотъемлемых факторов жизнеобеспечения человека, животного и растительного мира на Земле. Радиационное поле Земли – одна из опреде – ляющих компонент земной экосистемы и биосферы, для поведения которых характерно взаимодействие отдельных компонент с проявлением синергизма (обратных связей, которые иногда приводят к взаимоусилению различных процессов). Поле солнечного излучения влияет на механизмы изменчивости (динамические процессы: циркуляция, конвекция, турбулентный перенос;

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

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

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

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

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

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

В задачах атмосферной оптики, климата, метеорологии и прогноза по годы впервые численные методы расчета радиационных характеристик были предложены профессором Е.С. Кузнецовым. Дальнейшее развитие матема тических методов в теории переноса осуществлено М. В. Масленниковым, Т. А. Гермогеновой, М. С. Малкевичем, Е. М. Фейгельсон, Л. М. Романовой, М. Г. Кузьминой, Т. А. Сушкевич и другими учениками Е. С. Кузнецова.

Научная школа Е. С. Кузнецова повлияли на автора, что отразилось на содержании книги.

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

Список литературы содержит только публикации автора совместно с учени ками и коллегами.

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

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

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

В книге изложен один из эффективных подходов к решению широкого класса многомерных плоских задач с горизонтальными неоднородностями на основе рядов теории возмущений и интегрального преобразования Фурье – – метод функций влияния (ФВ) и пространственно-частотных характеристик (ПЧХ), позволивший сконструировать полуаналитические решения в виде интегралов суперпозиции и фурье-интегралов. Исходя из трехмерной плоской общей краевой задачи теории переноса аналитически, на физическом уровне строгости, построены комплексные уравнения переноса для линейных и нелинейных ПЧХ, а также уравнения и представления через ПЧХ для фундаментальных скалярных и векторных решений линейных и нелинейных приближений теории возмущений по вариациям коэффициентов, источников и граничных условий в горизонтальных плоскостях. Метод ПЧХ позволяет редуцировать трехмерную по пространству, пятимерную по фазовому объему краевую задачу для уравнения переноса к параметрическому набору рекур рентных, одномерных по пространству, трехмерных по фазовому объему задач для комплексной, непрерывной, бесконечно дифференцируемой по параметру функции, т. е. параметрических одномерных комплексных уравнений переноса, различающихся источниками в зависимости от вида «возмущения».

В рамках концепции линейных систем построены оптические передаточные операторы (ОПО), связывающие явно горизонтальные вариации интенсивности (вектора Стокса) излучения с вариациями коэффициентов, граничных условий, источников. Ядра ОПО – универсальные характеристики системы переноса – – – ПЧХ и ФВ, инвариантные относительно горизонтальных вариаций параметров системы. Тем самым разделяются проблемы вычисления ядра оператора и численного расчета решения исходной задачи в виде функционалов при различных пространственных структурах входных данных.

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

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

8 Предисловие Развит ИМХ и проведен качественный анализ (разрешимость, сходимость итераций, асимптотика, непрерывность, дифференцируемость по параметру, симметрия, влияние разностной сети и т. п.) решения комплексного уравнения переноса. Смоделировано влияние оптических параметров среды на ПЧХ и ФВ, описывающие передаточные свойства среды.

Предложены алгоритмы корректного расчета прямого ОПО при различных структурах вариаций альбедо и коэффициентов рассеяния, в том числе для набора «тест-объектов», и обратного ОПО при произвольных вариациях интенсивности с привлечением аппарата интегральных преобразований обоб щенных функций и определением ПЧХ как амплитудно- и фазо-частотных характеристик (АЧХ и ФЧХ).

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

Замечание о нумерации формул: внутри главы нумерация формул сплош ная, при этом первая цифра – номер главы.

– Текст монографии в TEX подготовлен автором. Искренне благодарю С. А. Стрелкова за активное участие в работах по поляризации, высоко профессиональные ценные консультации и обеспечение условий на ПК для конвертирования текстов из других редакторов и языков, а также Е. В. Влади мирову за помощь в подготовке текстов в TEX. Особую благодарность выражаю С. В. Максаковой, при участии которой в сложные для российской науки 90-е годы XX века проводилась большая научная работа по проектам РАН, ГКНТ и Миннауки РФ, грантам РФФИ и были подготовлены публикации, вошедшие в монографию.

Основные результаты получены автором в Институте прикладной мате матики имени М. В. Келдыша Российской академии наук.

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

Т.А. Сушкевич ГЛАВА Одномерные плоские задачи.

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

§ 1.1. Скалярная плоская задача Рассмотрим задачу переноса монохроматического излучения в плоскопа раллельном вертикально-неоднородном плоском слое трехмерного евклидова пространства. На верхнюю границу слоя z = 0 под углом 0 = arccos 0 падает внешний параллельный поток мощности S в азимутальной плоскости 0, s0 = {0, 0 }. Внутри слоя и на его границах могут быть диффузные источники излучения. Предполагается, что рассеяние в акте происходит без изменения частоты. Подстилающая поверхность на границе слоя z = H (H – толщина– слоя) отражает падающее изнутри излучение по заданному закону. Требуется определить интенсивность многократно рассеянного излучения внутри слоя, а также излучения, отраженного и пропущенного слоем.

Направление распространения излучения s в точке слоя z [0, H] описываем сферическими координатами и, где [0, 2] – азимут (за – = 0 принимается азимут внешнего потока излучения 0 = 0), [0, ] – – 10 Глава 1. Одномерные плоские задачи зенитный угол, отсчитываемый от направления внутренней нормали к плос кости z = 0 с единичным вектором n0, совпадающей с осью z. Направления с 0 = + 90 соответствуют нисходящему, пропущенному излучению, а с 90 = 180 – восходящему, отраженному излучению. Вводим – обозначения:

+ = cos + (0, 1], = cos [1, 0) = cos [1, 1], и множества:

+ = {(, ) : [0, 1], [0, 2]} = [0, 1] [0, 2] – полусфера направлений нисходящего излучения;

– = {(, ) : [1, 0], [0, 2]} = [1, 0] [0, 2] – полусфера направлений восходящего излучения;

– = + = {(, ) : [1, 1], [0, 2]} = [1, 1] [0, 2] – единичная сфера направлений распространения излучения.

– Интенсивность излучения (z,, ) на уровне z [0, H] в направлении s = {, } находится как решение общей краевой задачи для интегродиффе ренциального уравнения переноса Dz = S + F (z,, ), (1.1) = S (s s0 ) + f0 (, ), = q RH + fH (, ).

0 H Решение задачи (1.1) в геометрии z,, определяется на множестве z = {(z,, ) : z [0, H], s } = [0, H] [1, 1] [0, 2].

Дифференциальный оператор переноса и интеграл столкновений Dz S s (z) (z,, ) (z, cos ) ds, + t (z), z где ds = d d, в качестве коэффициентов содержат физические харак теристики среды: коэффициент рассеяния s (z), коэффициент ослабления (экстинкции) t (z) = s (z)+a (z), в который входит коэффициент поглощения a (z), и индикатрису рассеяния (z, cos ) с нормировкой (cos ) ds = 1;

(s, s ) ds = (cos ) d cos =, где угол рассеяния из направления s = {, } в направление s = {, } находится из соотношения cos = + sin sin cos( ).

1.1. Скалярная плоская задача Если ввести оптическую толщину z (z) = t (z ) dz, то с учетом равенства дифференциалов d = t (z)dz в уравнении (1.1) можно перейти к пространственной переменной :

D S ( ) (,, ) (, cos ) ds ;

+ 1, (1.2) альбедо акта рассеяния s ( (z)) ( ) =.

t ( (z)) В геометрии,, решение определяется на множестве = {(,, ) : [0, H ], s } = [0, H ] [1, 1] [0, 2] с оптической толщиной слоя H (z = H).

Для средней по азимуту или нулевой азимутальной гармоники интенсив ности излучения 0 (z, ) = (z,, ) d интеграл столкновений S0 = s (z) 0 (z, ) 0 (z,, ) d вычисляется через нулевую гармонику индикатрисы рассеяния 0 (z,, ) = (z, cos ) d.

Граничные условия в задаче (1.1) определяются на множествах s + } {(z,, ): z = 0, в геометрии z,,, {(,, ): + } s = 0, в геометрии,,, 0 = {(z, ): (0, 1]} в геометрии z,, z = 0, {(, ): (0, 1]} в геометрии, ;

= 0, s } {(z,, ): в геометрии z,,, z=H, {(,, ): } s в геометрии,,, = H, H = {(z, ): [1, 0)} в геометрии z,, z=H, {(, ): [1, 0)} в геометрии,.

= H, 12 Глава 1. Одномерные плоские задачи Если на уровне z = h, 0 h H, располагается граница раздела двух слоев, то для задания граничных условий вводятся множества:

s + } в геометрии z,,, {(z,, ): z = h, h + = s + } в геометрии,, ;

{(,, ): = (h), s } в геометрии z,,, {(z,, ): z = h, h = s } в геометрии,,.

{(,, ): = (h), Закон отражения излучения от подложки слоя описывается оператором || (H,, ) (,, ) d d.

RH = + Для ламбертовой поверхности (,, ) = 2||, q – альбедо, – при френелевском отражении (,, ) = 2( + )( )Z( ) ;

RH = Z(||)(H,, ).

Для удобства и большей наглядности величины, определенные для на правлений с = + (0, 1], помечаются меткой «+», а для = [1, 0) – – меткой «», например:

I + I(, +, ), I I(,, ).

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

Для учета пространственной неоднородности слоя по высоте в виде плоскопараллельных слоев вводятся зоны двух типов, а именно:

1) зоны с равномерным шагом пространственной разностной сети;

их границы обозначаются через zl1, l1 = 1 L1, zl1 =1 = 0, zl1 zl1 +1, zl1 =L1 = H;

шаги разностной сети zl1 = (zl1 +1 zl1 )/Nl1, где Nl1 – количество равных – интервалов по z в l1 -зоне;

размерность пространственной сети zk равна L1 Nl1 + 1;

K= l1 = 2) зоны с постоянными по z индикатрисами рассеяния (z, cos ) (в частности, элементами матриц рассеяния (z, cos ));

их границы обозначаются через zl2, l2 = 1 L2, zl2 =1 = 0, zl2 zl2 +1, zl2 =L2 = H;

для направлений с 0 полагаем (z, cos ) = l2 (cos ), если zl2 z zl2 +1 ;

для направлений с 0 полагаем (z, cos ) = l2 (cos ), если zl2 z zl2 +1.

Не нарушая общности физических постановок задачи, считаем, что границы l2 -зон входят в число границ l1 -зон, т. е. L1 L2.

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

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

МОДЕЛЬ 1: s (z) и t (z) постоянны в l1 -зонах;

полагаем t (z) = tl1, s (z) = sl1 для zl1 z zl1 +1, если 0, и zl1 z zl1 +1, если 0.

МОДЕЛЬ 2: на каждом интервале разностной сети [zk, zk+1 ] коэффициенты считаются постоянными: s (z) = sk, t (z) = tk, k = 1 (K 1).

МОДЕЛЬ 3: значения t (z) и s (z) на интервале [zk, zk+1 ] задаются табли цей (или формулой), т. е. на одном шаге zk = zk+1 zk задается несколько значений t (z) и s (z);

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

+ Для всех точек разностной сети zk формируются массивы sk для + направлений + и sk для направлений ;

на границах слоев индекс sk определяется из условия zk z zk+1, а sk – zk z zk+1. Эта процедура – вводится для учета разрывов коэффициентов s (z) и t (z) на границах зон в вычислениях интегралов столкновений: при интегрировании по s + в точке zk берется s (z) слоя, прилегающего сверху к zk, а по s – слоя, – прилегающего снизу к zk.

В случае задачи, приведенной к оптической толщине, условно делается замена z на, t на 1, s ( ) на ( ), H на H.

Анизотропия рассеяния характеризуется такими величинами как (cos = 1) ± (1) = (cos = 1) – отношение вероятностей рассеяния в направлениях вперед и назад;

– 1 ± = (cos ) d cos (cos ) d cos – отношение вероятностей рассеяния в переднюю и заднюю полусферы;

– 1 g1 = (cos ) cos d cos (cos ) d cos 1 – средний интегральный косинус угла рассеяния, причем 3g1 = 1, где 1 – – – коэффициент разложения индикатрисы рассеяния по полиномам Лежандра N P0 (1) = 1, (1.3) (cos ) = P (cos ), P1 (1) = cos.

= 14 Глава 1. Одномерные плоские задачи Суммарная индикатриса рассеяния может быть задана молекулярной R и аэрозольной a индикатрисами раздельно с соответствующими молекулярной R и аэрозольной a,s оптическими толщинами (или весами) и вычисляется по формуле (послойно) a,s a () + R R (), s a,s + R.

() = s s Если значения коэффициента s (z) постоянны по слоям, то в пределах слоев () = a,s a () + R R (), s a,s + R.

s s Каждая из индикатрис рассеяния входит в оптическую модель либо в виде таблицы значений для набора углов рассеяния [0, 180 ], либо в виде коэффициентов разложения по полиномам Лежандра (1.3), либо в виде аналитической формулы с параметрами, например, индикатриса Хеньи Гринстейна 1 g (1.4) HG () = (1 + g 2 g cos ) 2 3/ или рэлеевская (1 + cos2 ). (1.5) R () = 1.1.1. Алгоритм интегрирования по характеристике. Разностная сеть.

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

(1.6) Dz = F (z, ), = f0 (), = fH ().

0 H Сначала решается однородное уравнение Dz = 0 методом разделения переменных (для удобства полагаем (z) t (z)):

d (z) (z) = dln|| = dz, dz с учетом физического условия 0. Из равенства дифференциалов следует, что неопределенные интегралы отличаются лишь на константу:

(z) (z) = C exp dz.

ln || = dz + ln C, Если при z = z0 = 0, то z (u) (z) = 0 exp du.

z 1.1. Скалярная плоская задача Для решения неоднородного уравнения с граничными условиями восполь зуемся методом изменения произвольной постоянной Лагранжа ( фиксиро вано):

(z) (z) = C(z) exp dz.

Производную d dC (z) (z) exp dz = C(z) dz dz подставляем в уравнение (1.6) и получаем уравнение для определения константы (z) dC F (z, ) exp dz =.

dz Пусть = + (0, 1]. Тогда z t 1 C(z) = + F (t,+ )exp + (u) du dt + C1 ;

0 z z z 1 1 (z,+ ) = + F (t,+ )exp + (u) du dt + C1 exp + (u) du.

t 0 Значение константы C1 находится из граничного условия на 0 при z = 0:

C1 = f0 (+ ). Следовательно, z z 1 F (t, + ) exp (u) du dt + (z, + ) = + + t z + f0 (+ ) exp + (u) du. (1.7) Пусть = [1, 0). Тогда z H 1 F (t, ) exp (u) du dt + C2.

C(z) = | | | | t H 16 Глава 1. Одномерные плоские задачи Из граничного условия на z = H определяется C2 = fH ( ) и H t 1 F (t, ) exp (u) du dt + (z, ) = | | | | z z H + fH ( ) exp (u) du. (1.8) | | z Формулы (1.7) и (1.8) являются результатом интегрирования по харак теристике = const дифференциального оператора одномерного плоского кинетического уравнения Dz.

Для построения конечно-разностного аналога формул (1.7) и (1.8) вводим разностную сеть: по пространственной переменной zk, k = 1 K, z1 = 0, zK = H, с шагом zk = zk+1 zk, по угловой переменной + (0, 1], j j = 1 J +, и [1, 0), j = 1 J, и сеточную функцию j + (zk, + ), (zk, ).

j j kj kj Для = + 0 на границе 0 (k = 1): + = f0 (+ ) = f0j, а для k = 2K, j j 1j пользуясь свойством аддитивности экспонент, приходим к представлению zk+1 zk+ 1 F (t,+ )exp (u) du dt + + = + + j k+1,j j j t zk+1 zk+ zk 1 1 + f0j exp (u) du = F (t,+ )exp (u) du dt + + + + j j j j t 0 zk+1 zk= 1 F (t,+ )exp (u) du dt + + + + j j j z t k zk+ zk 1 + f0j exp (u) du exp (u) du = + + j j z k zk+1 zk zk 1 1 = exp (u) du F (t,+ )exp + (u) du dt + + + j j j j zk t 1.1. Скалярная плоская задача zk+1 zk+ zk 1 1 + f0j exp (u) du + + F (t,+ )exp + (u) du dt = + j j j j z t k zk+1 zk+1 zk+ 1 1 = + exp (u) du + F (t,+ )exp (u) du dt.

+ + + j kj j j j zk zk t (1.9) На границе H для = 0 имеем = fH ( ) = fHj и аналогичные j j Kj преобразования для k = (K 1) 1 приводят к выражению zk+ du + = k+1,j exp | | kj j zk zk+1 t 1 F (t, ) exp (u) du dt. (1.10) + j |j | |j | zk zk Представления (1.7) и (1.9) для 0 и (1.8) и (1.10) для взаимно-однозначны на сети zk. Но вторые представления в виде рекуррентных соотношений позволяют ввести существенную экономию объема вычислений и организовать последовательность расчета с переходом к следующему уровню по z от предыдущего уровня, а не от самой границы. Отметим, однако, что в задачах с учетом селективного поглощения излучения, описываемого не коэффициентами поглощения, а спектральной функцией пропускания, зависящей от эффективной оптической толщины вдоль всего луча от точки пересечения им границы до расчетной точки, используются алгоритмы, основанные на представлениях (1.7) и (1.8).

Для построения конечно-разностного аналога формул интегрирования по характеристикам необходимо ввести аппроксимацию коэффициента ослабле ния (z) и функции источника F (z, ) на отрезках [zk, zk+1 ]. Наиболее часто применяется кусочно-постоянная аппроксимация F (zk+1, ) + F (zk, ) (z) = k, F (z, ) =.

В этом приближении интегралы в формулах (1.9) и (1.10) вычисляются аналитически и конечно-разностный аналог приводится к явной двухточечной схеме:

+ ++ ++ + + k+1,j = kj Dkj + Fkj Akj + Fk+1,j Ckj, = Dkj + Fkj Ckj + Fk+1,j A kj k+1,j kj 18 Глава 1. Одномерные плоские задачи с положительными коэффициентами +() 1 Dkj k +() +() +() = exp Dkj, Ckj = Akj =, k = k zk.

|+() | 2k Эта схема является монотонной и устойчивой.

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

zk+1 zk+ 1 exp + (u) du dt, A+ = Ckj = + 2+ kj j j t zk zk+1 zk+ t 1 exp (u) du dt, A = Ckj = k = (u) du, 2| | kj |j | j zk zk zk вычисляются численно с учетом конкретного вида функции (z) на ос нове квадратурных формул. При кусочно-линейной аппроксимации функции источника коэффициенты схемы также рассчитываются путем численного интегрирования:

zk+1 zk+ zk+1 t 1 exp + (u) du dt, A+ = + kj zk j j zk t zk+1 zk+ t zk 1 exp + (u) du dt, + Ckj = + zk j j zk t zk+1 t t zk 1 exp (u) du dt, A = | | kj |j | zk j zk zk zk+1 t zk+1 t 1 exp (u) du dt.

Ckj = | | |j | zk j zk zk Если (z) = k для z [zk, zk+1 ], то +() 1 Dkj +() +() Akj Ckj =, k +() +() |j | |j | +() +() Dkj 1+ Akj =.

k k k 1.1. Скалярная плоская задача При симметричной разностной сети, когда + = | |, j j + A+ = A, + Dkj = Dkj, Ckj = Ckj.

kj kj Для = j=0 = 0 расчетные формулы получаются предельным переходом в исходном уравнении с учетом возможного разрыва решения на границе zk.

Поэтому для направлений = 0 рассчитываются два значения функции – – одно отвечает подходу к границе со стороны направлений s +, а другое –– s :

F+ F + = k,j=0, = k,j=0. (1.11) k,j=0 k,j= k1 k Особое внимание обращается на расчет однократного рассеяния. В задачах с малой оптической толщиной вклад однократного рассеяния практически всегда существен. Некоторые инженерные методы приближенного расчета радиационных полей и исходные уравнения в ряде обратных задач основы ваются на приближении однократного рассеяния. Для исследования угловых и пространственных структур полей излучения в нашей практике прочно закрепился общий подход к численному решению всех уравнений переноса с разделением вкладов нерассеянной, однократно и многократно (начиная со второй кратности) рассеянной компонент излучения. При этом обычно для расчета двух первых компонент привлекаются явные выражения. С помощью дискретного представления формул интегрирования по характеристикам в случае задачи с внешним мононаправленным (например, солнечным) потоком получаем следующие выражения для компоненты однократного рассеяния 1 :

если + = 0, то j + + S 0 s,k+1 tk (j, 0 ) k 1+ = 1+ exp + + + 0 + k+1,j kj t,k+ j j (zk ) k k exp exp exp, + 0 j sk tk (, 0 ) k S + 0 j 1 = 1 exp | | kj k+1,j |j | + 0 tk j 1 (zk ) exp 1 exp + k ;

|j | 0 если + = 0, то j + S s,k+1 tk (0, 0 ) k (zk+1 ) 1+ 0 = 1+0 exp exp +.

k+1,j kj 0 0 Величины tk описаны в следующем разделе.

Для однородного по высоте слоя, когда обычно делается переход к пространственной переменной – оптической толщине и вводится альбедо акта – 20 Глава 1. Одномерные плоские задачи рассеяния = s /t, получаем (при ламбертовой подложке с альбедо q) S 0 (z) (z) exp exp (, 0, 0 ), 1+ (z,, ) = 0 S 0 (z) exp 1 (z,, ) = || + 0 1 1 exp (H (z)) (, 0, 0 ) + + || 1 (z) + qS 0 exp (H (z)), + || 0 S (z) (0, 0, 0 ) exp 1+ (z, 0, ) =.

0 Отсюда нетрудно получить выражения для интенсивности пропущенного ( (z) = (H), s + ) и отраженного ( (z) = 0, s ) излучения в приближении однократного рассеяния.

Аналитические выражения используются для тестирования более сложных вычислительных алгоритмов. Отметим, что при 0 = 1, т. е. при нормальном падении внешнего потока, задача является азимутально-симметричной и достаточно рассчитывать только нулевую азимутальную гармонику.

1.1.2. Учет азимутальной зависимости. Азимутальную зависимость реше ния обычно учитывают двумя способами. В первом способе задача сводится к нахождению гармоник разложения решения по азимуту в ряд Фурье. Про иллюстрируем этот способ на решении задачи с внешним мононаправленным потоком:

Dz = S + F (z,, ), (1.12) = S ( 0 )( 0 ), = R + fH (, ).

0 H Для общности учтем три типа отражения на границе z = H:

+ + = 1.

R = Rq + Rf r + RH, При ламбертовом законе 2 q Rq = d (H,, ) d ;

0 при зеркальном, френелевском отражении Rrf = Z()(H,, ) ;

1.1. Скалярная плоская задача при неортотропном отражении 2 RH = d (H,, )H (cos H ) d.

0 Решение ищем в виде суперпозиции = 0 +. В аналитическом виде выделяем S ( )( ) exp (z), s + ;

0 0 (z,, ) = S Z()( + )( ) exp (z) 2 (H), s, 0 – нерассеянную и зеркально отраженную компоненту с -особенностью и – сглаженную через интеграл столкновений переводим ее во внутренний источ ник B0 S0. Многократно рассеянная компонента (z,, ) определяется численно из краевой задачи (B S):

= 0, Dz = B + B0 + F, = fH + R0 + R.

0 H Разложением в ряды Фурье по азимутальным гармоникам проводится разделение угловых переменных:

c (z, ) cos m + s (z, ) sin m. (1.13) (z,, ) = 0 (z, ) + m m m= Аналогично разлагаются источники и граничные условия на z = H:

1 c s F (z,, ) = F0 (z, ) + Fm (z, ) cos m + Fm (z, ) sin m, m= 1 c s fH (, ) = fH0 () + fHm () cos m + fHm () sin m, m= 1 c s R0 = r00 () + r0m () cos m + r0m () sin m, m= 1 c s R = r0 () + rm () cos m + rm () sin m.

m= Коэффициенты разложения граничных условий определяются с учетом пред ставления -функции:

( 0 ) = 0 + c s m cos m + m sin m.

m= 22 Глава 1. Одномерные плоские задачи Индикатрисы рассеяния разлагаются либо по полиномам Лежандра (z, cos ) = (z) P (cos ) = с использованием теоремы сложения, либо по полиномам Чебышева первого рода степени m:

m (z,, ) cos m( ) = (z, cos ) = m (z,, ) Tm (y), m=0 m= где y = cos( ), Tm (y) = cos m( ). В этом случае c s B(z,, ) = B0 (z, ) + Bm (z, ) cos m + Bm (z, ) sin m, m= B0 (z, ) 2s (z) 0 (z, ) 0 (z,, ) d ;

Bm (z, ) s (z) c c (z, ) m (z,, ) d ;

m Bm (z, ) s (z) s s (z, ) m (z,, ) d ;

m c s B0 (z,, ) = B00 (z, ) + B0m (z, ) cos m + B0m (z, ) sin m, m= B00 (z, ) a0 (z) [0 (z,, 0 ) + b0 (z) 0 (z,, 0 )], (z) a0 (z) S s (z) exp, (H) (z) b0 (z) Z(0 ) exp, B0m (z, ) a0 (z) [m (z,, 0 ) cos m0 + b0 (z) m (z,, 0 ) cos m0 ], c B0m (z, ) a0 (z) [m (z,, 0 ) sin m0 + b0 (z) m (z,, 0 ) sin m0 ].

s 1.1. Скалярная плоская задача Подставляя разложения в краевую задачу и приравнивая выражения при одинаковых гармониках, находим краевую задачу для 0-гармоники Dz 0 = 2s (z) 0 (z, ) 0 (z,, ) d + F0 (z, ) + + a0 (z) [0 (z,, 0 ) + b0 (z) 0 (z,, 0 )], = 0, 0 0 = fH0 () + r0 () ;

0 H краевые задачи для нахождения амплитуд m-гармоник по косинусам Dz c = s (z) c (z, ) m (z,, ) d + Fm (z) + c m m + a0 (z) cos 0 [m (z,, 0 ) + b0 (z) m (z,, 0 )], c m c c c c = 0, = fHm () + r0m () + rm ();

m 0 H краевые задачи для определения амплитуд m-гармоник по синусам Dz s = s (z) s (z, ) m (z,, ) d + Fm (z, ) + s m m + a0 (z) sin 0 [m (z,, 0 ) + b0 (z) m (z,, 0 )], s m s s s s = 0, = fHm () + r0m () + rm ().

m 0 H При ламбертовой подложке в граничных условиях остаются только два члена: (H) r00 = qS 0 exp r0 = 2q 0 (H, ) d.

, Введем tm (, ) 2m (, ) с правой частью + (cos ) Tm (y) dy m = 0, 1, 2,... (1.14) m (, ) =, 1 y Тогда t0 (, ) tm (, ) = m, = 0, m 0.

0 (, ) = m (, ) = 2 24 Глава 1. Одномерные плоские задачи В результате уравнения для амплитуд гармоник записываются единообразно:

Dz 0 = s 0 t0 (, ) d + F0 (z, ), Dz c = s c tm (, ) d + Fm (z, ), c m m Dz s = s s tm (, ) d + Fm (z, ).

s m m Если 0 = 0, то имеет место азимутальная симметрия () = () = (2 ) (1.15) и все амплитуды синусоидальных гармоник являются нулевыми. Такая форма записи позволяет учитывать разложение индикатрисы рассеяния по полиномам Лежандра: ( m)! m m tm (, ) = P () P ( ).

( + m)!

=m Коэффициенты разложения индикатрисы рассеяния вычисляются с помощью квадратуры Чебышева-Меллера 1 N 2k 1 v(x) dx v cos.

N N 1 x2 k= Приближенные значения N m (, ) (cos k ) Tm (yk ) ;

cos k = + sin sin yk, N k= (2k 1) m(2k 1) yk = cos, Tm (yk ) = cos.

2N 2N В частности, для 0-гармоники N 2 + 1 P () P ( ) t0 (, ) = (cos k ).

2 N k= = Гармоники функции источника – интегралы по на отрезке [1, 1] – вычис – – ляются по квадратурной формуле (гауссовой, Лобатто или какой-либо другой) с разделением четырех квадрантов в области = [1, 1] [1, 1]:

J+ J B(z, ± ) + a+ + tk (±, + ) a tk (±, ).

= sk + sk j j j j kj j j kj j j =1 j = 1.1. Скалярная плоская задача Для сокращения объема вычислений сеточные значения гармоник инди катрис рассеяния tk (j, j ) рассчитываются вне итераций и запоминаются в архиве моделей с последующим вызовом при необходимости в оперативную память. Естественно, что на отрезках [1, 0] и [0, 1] можно брать квадратуры разного порядка. Но, как правило, включается дополнительно еще один узел, совпадающий с направлением внешнего потока 0. При большой анизотропии рассеяния для повышения точности расчета интегралов столкновений берутся гауссовы квадратуры, пересчитанные на отрезки [0, 0 ] и [0, 1]. Такая процедура позволяет «сгустить» разностную сеть около 0 – направления – сильной анизотропии нисходящего излучения. Применение квадратур со сплайнами повышает точность и позволяет более детально воспроизводить решение по окончании счета путем пересчета решения на другую угловую сеть с помощью интеграла столкновений.

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

m(n) Sn (z,, ) (z,, ), m(n) m=M где частичные суммы n Sn = m (z, ) cos m, M = mmin, n mmax, m(n) mmax, m=M с критерием сходимости ряда Фурье:

Sm (z,, ) Sm1 (z,, ) m = max.

Sm (z,, ) {z,,} Скорость убывания гармоник m (z, ) с ростом оптической толщины (H) тем выше, чем больше ее номер m. На больших оптических расстояниях от источника решение (z,, ) с достаточной точностью определяется несколькими первыми (или даже нулевой) гармониками. Этот результат близок к асимптотическому. По мере приближения к источнику (к верхней границе) вклад старших гармоник растет. Вблизи источников или в задачах со слоями малой оптической толщины при сильно анизотропном рассеянии сумма (1.13) содержит почти столько же членов, сколько получается гармоник в разложении индикатрисы рассеяния по полиномам Лежандра (1.3).

Зафиксируем некоторые значения z0 [0, H] и 0 [0, ] и проследим за установлением углового профиля gM () M (z0,, 0 ) по мере увеличения числа гармоник M в представлении решения конечным рядом, содержащим M 26 Глава 1. Одномерные плоские задачи гармоник Фурье. Можно указать такие значения M1 M0 M2, что функции gM1 () и gM2 () незначительно отклоняются от gM0 () в противофазе. Даль нейшее увеличение M приводит ко все большим отклонениям gM () от gM0 ().

Такая ситуация изображена на рис. 1.1. Каждая гармоника рассчитывалась с точностью = 103. Причина такого явления заключена в некорректности операции суммирования ряда Фурье: ошибки расчета каждой из амплитуд гармоник m (z, ) при вычислении суммы (при больших M ) складываются, а это приводит к большому отклонению частичной суммы M (z,, ) от точного решения (z,, ). За приближенное решение задачи следует принять частичную сумму с M = M0, при этом погрешность оценивается из сравнения решений при M = {M1, M0, M2 }. На практике M0 N, где N – наименьшее – число коэффициентов разложения индикатрисы рассеяния в ряд по полиномам Лежандра, позволяющее восстановить индикатрисы рассеяния с погрешностью не хуже 1%.

Если ввести отклонение W (z,, ) = (z,, ) M (z,, ) и из уравнения для компоненты кратного рассеяния (z) Dz = s (z) ds + S s (z) (z, 0 ) exp вычесть уравнение для решения M в виде конечного ряда Фурье (z) Dz M = s (z) M M ds + S s (z) M (z, 0 ) exp, то приходим к уравнению для отклонения (1.16) Dz W = s (z) W (z,, ) (z, ) ds + F1 (z,, ), (z) F1 (z,, ) S s (z) (z, 0 ) M (z, 0 ) exp, так как из-за ортогональности функций cos m с разными значениями m M ds M M ds = 0.

Для коррекции приближенного решения M бывает достаточно вычислить отклонение W1 (z,, ) посредством интегрирования по характеристике уравне ния (1.16) с правой частью – источником F1 (z,, ), отвечающим погрешности – приближения однократного рассеяния из-за разложения индикатрисы в конечный ряд Фурье.

Во втором способе учета азимутальной зависимости решения по азимуту вводится разностная угловая сеть и интеграл столкновений вычисляется 1.1. Скалярная плоская задача · r I 5 5,47 · 108 l I 1,54 · 0,1 0,6 0, 0,2 0, Рис. 1.1. Угловые профили интенсивности () на высоте 30 км для 0 = 60, M = I – в приближении однократного рассеяния при M = 9 (кривая 1), 40 (2), 48 (3), 56 (4), 64 (5) – и 70 (6);

II – с учетом многократного рассеяния при M = 5 (7), 7 (8), 30 (9), 40 (10), 50 (11), 60 (12), – 70 (13) и 90 (14) 28 Глава 1. Одномерные плоские задачи с помощью квадратурных формул на поверхности единичной сферы.

Прямолинейный расчет по квадратурной формуле J + +J I (1.17) Bkji = sk ai bj kjij i kj i, i =1 j = где bj и j – веса и узлы квадратуры по на отрезке [1, 1], ai и i – – – веса и узлы квадратуры по азимуту на отрезке [0, 2], приводит к большим затратам машинного времени.

1.1.3. Алгоритм расчета интеграла столкновений. Решение задачи (1.6) при fH = 0 симметрично относительно азимутальной плоскости = 0.

Обычно берется 0 = 0, и расчетные формулы заметно упрощаются. Легко показать, что в этом случае имеет место азимутальная симметрия (1.15) для решения задачи (1.6) и интеграл столкновений сводится к отрезку [0, ]:

2 (z,, ) (z,,, ) d = B(z,, ) = s (z) d (z,, ) [ + + ] d, (1.18) = s (z) d где (z,,,, ) (z,,, ), + (z,,,, ) (z,,, + ).

+() Вычисление сеточных значений интегралов B +() (zk, j, i ) можно проводить с помощью квадратур на единичной сфере:

J+ J I I + bi kjij i + i a bi kjij i i, (1.19) Bkji = sk + sk kj j kj j =1 i =1 j =1 i = +() +() где kjij i (zk, j, j, i i );

и – веса и узлы квадратурной – aj j формулы для + [0, 1] и [1, 0] (не обязательно симметричной и одинаковой размерности на [0, 1] и [1, 0], т. е. допускается J + = J, + = | |);

bi и i – веса и узлы квадратуры по азимуту на отрезке – j j [0, 2].

Если zk не является границей слоев с разными характеристиками s (z), + то sk = sk = s (zk.) В противном случае при интегрировании по полусфере + берется значение + = (z 0), т. е. соответствующее характеристике sk sk – слоя, прилегающего сверху, а по полусфере – sk = s (zk ), отвечающее нижнему слою.

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

Для повышения экономичности алгоритма без существенных ограничений на общность разностная сеть по азимуту вводится специальным образом, а именно: вводим узлы по азимуту m, m = 1 M, так что 1 = 0, m m+1, M = /2, а остальные узлы определяются из следующих соотношений:

i= 1M;

для 0 /2 i = m, i = m, i = m, i = 2M + 1 m, i = (M + 1) 2M ;

для /2 i = (2M + 1) 3M ;

для 3/2 i = 2M + m, i = + m, 2 i = 2 m, i = 4M + 1 m, i = (3M + 1) 4M.

для 3/2 С учетом свойства азимутальной симметрии (1.15) возможны два варианта организации расчета, которым соответствует I = {2M, если [0, ];

если [0, 2]}.

4M, Принимая во внимание то обстоятельство, что индикатриса рассеяния зависит от азимута только через cos( ), получаем следующие соотношения симметрии:

( ) = ( );

1) ( + ) = (( + )) = (2 ( + )) = (( + ) 2);

2) ( ( + )) = (( + ) );

3) ( + ( )) = ( ( )) = (( ) ) = (( ) ).

4) С помощью формул приведения для косинуса устанавливаем, что этим соотношениям соответствуют всего четыре значения угла рассеяния, вычис ляемые по формулам 1 2 1 2 cos( ), cos 1 = + 1 2 1 2 cos( + ), cos 2 = + cos 3 = 1 2 1 2 cos( + ), cos 4 = 1 2 1 2 cos( ). (1.20) Обозначим (cos ), = 1 4.

В расчетах интегралов столкновений большая часть объема вычислений и времени счета приходится на определение значений индикатрисы рассе яния (zk, j, j, i, i ), которые находятся либо суммированием рядов с полиномами Лежандра, либо интерполяцией между заданными табличными значениями, либо по формулам. Установленные выше соотношения позво ляют сократить объем вычислений в 4 раза путем расчета (z, cos ) при 30 Глава 1. Одномерные плоские задачи фиксированных наборах аргументов {zk, j, j, i, i } только для четырех углов рассеяния (1.20) вместо 16 значений, входящих в формулу (1.18).

Учитывая, что разностная сеть по азимуту вводится через сеть m на [0, /2], экономия времени может быть получена, если в итерационном процессе предварительно рассчитать и использовать массивы значений cos+ cos(m + m ), cos cos(m m ), m, m = 1 M.

mm mm Для расчета индикатрис рассеяния, заданных с помощью таблиц значений (cos ) в отдельных узлах cos, в целях сокращения объема вычислений целесообразно воспользоваться монотонностью изменения косинуса угла рассеяния cos в зависимости от аргументов и, принимающих значения m и m. Циклы по дискретным наборам азимутов m и m организуются так, чтобы значения индексов m и m возрастали от 1 до M.

При фиксированном m = 1 M всегда 0, возрастает с ростом m, пока m m, cos(m m ) равен 1 при m = m, убывает с ростом m, когда m m ;

убывает монотонно с ростом m, cos(m + m ) принимает значения в пределах от +1 до 1.

С ростом m возрастает, пока m m, cos 1 имеет максимум при m = m, убывает, когда m m ;

cos 4 убывает монотонно с ростом m ;

cos 3 возрастает монотонно с ростом m ;

убывает, пока m m, cos 2 имеет минимум при m = m, возрастает, когда m m.

Эти выводы позволяют для вычисления (cos ), = 1 4, ускорить выбор интервалов интерполяции [cos 1, cos ], на которых cos cos cos, посредством запоминания индексов интервалов, в которые попадают cos при предыдущем значении индекса m, и просмотра таблиц cos для m + 1, начиная с соответствующих индексов. Кроме того, при m = 1 просмотр таблиц cos целесообразно начинать с = 1, если 0 cos 0, или с номера = N, если cos 0 (номер = N соответствует первому в таблице значению cos 0).

Поскольку обычно индикатрисы рассеяния (z, cos ) от z зависят по зонам, содержащим по несколько узлов сетки zk, рекомендуется самый внутренний 1.1. Скалярная плоская задача цикл при организации алгоритма расчета делать по zk и вычислять при фиксированных j, i, j, i значения (cos ) один раз для всех таких zk.


При фиксированных значениях = j и = j расчет азимутальных квадратур рекомендуется проводить одновременно в виде следующих сумм:

i = m, i = m, i = 1 M, для 0 /2 M ± ± (zk, ±, m )[1 (zk j, ±, cos(m m )) + (j ) bm j j 1k m = + 2 (zk, j, ±, cos(m + m ))] + ± (zk, ±, m ) (1.21) j j 3 (zk, j, ±, cos(m + m )) + 4 (zk, j, ±, cos(m m )) ;

j j i = m, i = (M + 1) 2M, для /2 M ± ± (zk, ±, m )[3 (zk, j, ±, cos(m + m )) + (j ) bm j j 2k m = + 4 (zk, j, ±, cos(m m ))] + ± (zk, ±, m ) (1.22) j j 1 (zk, j, ±, cos(m m )) + 2 (zk, j, ±, cos(m + m )).

j j Интегралы по зенитному углу :

J± J± ± ± ± ± a± a± (1.23) (j ), (j ).

j j 3k 1k 4k 2k j =1 j = Интегралы по полусферам + и :

± ± ± ± ± ± B 1/2 (zk, j, m ) = sk B 1/2 (zk, j, m ) = sk (1.24),.

3k 4k Оптимальная организация последовательности циклов: j ;

m ;

l2 –зоны для индикатрис рассеяния;

j ;

m ;

расчет, = 1 4;

zk в пределах l2 –зоны.

С учетом свойств симметрии и анизотропии индикатрисы рассеяния оказывается целесообразным использовать квадратурную формулу J + +J I +() +() Bkji = sk ai j kjij i (kj i + kj,i 1 + k,j 1,i + k,j 1,i 1 ), i =2 j = в соответствии с которой при расчетах интеграла на сфере функция +() +() ] (z,, ) аппроксимируется кусочно-постоянной в областях [j1, j [i1, i ], а интегралы j i (zk, j,, i ) d kjij i = d i j 1 32 Глава 1. Одномерные плоские задачи вычисляются более точно с привлечением дополнительной разностной сети или сплайн-интерполяции на отрезках [i1, i ] и [j1, j ].

Для повышения точности расчета значений индикатрисы рассеяния (1.3), заданной коэффициентами k, k = 1 N, разложения в ряд по полиномам Лежандра, применяется процедура суммирования, отвечающая схеме обратной прогонки:

(cos ) = b0 P0 + b1 (P1 B0 P0 ), P0 = 1, P1 = cos, 2n + 1 n Cn (cos ) = Bn (cos ) = cos,, n+1 n+ bN = bN +1 = 0, k = (N 1) 0.

bk = bk+1 Bk + bk+2 Ck+1 + k, При прямом суммировании ряда (1.3) положительные и отрицательные его члены суммируются отдельно, а полиномы Лежандра вычисляются по рекуррентной формуле (иногда с двойной точностью) 2 1 cos P1 (cos ) 2.

P (cos ) = P2 (cos ), Рекомендуется также вычисление индикатрис рассеяния (1.3) суммиро ванием частичных сумм по Фейеру. Необходимо отметить, что реальные индикатрисы рассеяния плохо разлагаются в ряды по полиномам Лежандра.

Для «хороших» в среднем индикатрис рассеяния при 50–60 членах разложения погрешность представления около 5%, а при 100 членах – около 2%. Отсюда – следует, что методы, использующие разложение индикатрисы рассеяния (1.3), заведомо вносят ошибки в решение, поскольку работают с приближенными ин дикатрисами, отличающимися от реальных. При этом чаще всего искажается ореольная часть индикатрисы, из-за которой требуются длинные ряды. На современных многопроцессорных компьютерах для моделирования переноса излучения в средах типа облаков и океана для описания сильной анизотропии рассеяния используют разложения решений и индикатрис рассеяния по полиномам Лежандра, содержащие по 1000–1500 членов. Необходимы меры для подавления таких искажений.

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

1.1. Скалярная плоская задача 1.1.4. Метод итераций. Сходимость простых итераций – естественных по – следовательных приближений по кратности рассеяния, т. е. по схеме N +1 = Dz S(N + 0 ), (1.25) где N – номер итерации, установлена В. С. Владимировым с привлечением – метода Келлога. Еще ранее она исследовалась Е. С. Кузнецовым и Е. Хопфом.

Приближенное решение определяется сходимостью последовательных прибли жений N +1 = (Dz ) [S (N + 0 )], (1.26) где (Dz ) – конечно-разностный оператор обращения дифференциальной 1 – части уравнения переноса, а S – разностный аналог интегрального оператора – S, зависящий от выбора схемы численного интегрирования на единичной сфере или отрезке [1, 1].

В некоторых задачах (например, с однородным слоем оптической тол щины ) можно оценить норму интеграла столкновений B = S:

s B [1 exp( )], t что позволяет оценить ошибку итерационного метода:

BN точн N.

1 B С развитием метода дискретных ординат неоднократно изучались вопросы сходимости итераций. В консервативных задачах сходимость итераций мед леннее, чем в задачах с поглощением. При большой оптической толщине и сильной анизотропии рассеяния простой итерационный процесс является медленно сходящимся. Было сосредоточено внимание на разработке методов ускорения сходимости, когда вместо (1.26) используется итерационный процесс N +1 = (Dz ) (S N +1/2 + F ), N +1/2 = U N, где оператор U определяется выбором метода ускорения. В качестве критерия окончания итерационного процесса проверяется условие на всей дискретной сети:

= {(zk, j, i ) : zk [0, H], ± [1, 1], i [0, 2]}, j |N +1 N | kji kji (1.27) max, N + kji где N – приближенное решение, полученное на итерации с номером N, – kji – – относительная точность итераций. Обычно критерий проверяется для со ставляющей многократно рассеянного излучения kji, начиная со второй кратности. В оптически нетолстых средах такое условие усиливает критерий сходимости итераций. При медленной сходимости итераций величина, стоящая в правой части критерия (1.27), вообще говоря, не дает представления об 34 Глава 1. Одномерные плоские задачи истинной близости N +1 к точному значению (z,, ) решения краевой kji задачи на : при малом отклонении двух последовательных приближений на большом числе итераций может накопиться достаточно заметное приращение.

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

Известно, что для рассеивающего и поглощающего слоя должно выпол няться соотношение баланса r + (1 q)t + a = S 0, (1.28) где поток излучения, отраженного слоем, r |J (z = 0)|, |J (z)| = (z,, ) d d, поток пропущенного излучения t Ft (z = H) = J + (z) + J0 (z), (z) J + (z) = J0 (z) = S 0 exp (z,, ) d d,, + чистый поток излучения F (z) = Ft (z) |J (z)|, поток поглощенного излучения a F (z = 0) F (z = H), F (z = 0) = S 0 |J (z = 0)|, F (z = H) = J + (z = H) + J0 (z = H) |J (z = H)|, |J (z = H)| = q[J + + J0 (z = H)].

Из (1.28) получаем соотношение для погрешностей r, t, a – откло – нений приближенного решения от точных значений:

r + (1 q)t + a = 0.

Отсюда следует, что ошибки r, t, a не должны быть одного знака.

Если альбедо q = 0, то J (z = H) = 0.

r + t + a = S 0 ;

Для консервативного случая имеют место два соотношения – законы – сохранения:

1) |J (z = 0)| + Ft (z = H) = S 0, т. е. сумма потоков излучения, пропущенного и отраженного слоем, равна внешнему падающему потоку;

1.1. Скалярная плоская задача 2) Ft (z) = Ft (z) + |J (z)| = const, т. е. полный поток излучения Ft (z) на всех уровнях постоянен.

Консервативные задачи допускают еще один интеграл:

2 (z,, ) 2 d = Ft (z) K(z) = d (z) + const, куда входит средний косинус индикатрисы рассеяния 2 1 (z) = d (z, cos ) d.

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

Решение уравнения переноса обладает локальными особенностями, кото рые достаточно подробно всегда можно учесть в алгоритмах с интегрированием по характеристикам. В этом состоит большое преимущество ИМХ при расчетах угловых распределений излучения. Но для задач с максимальным собственным числом, близким к единице, характерна медленная сходимость простых итераций. Ускорению сходимости итерационных методов решения линейных уравнений, в том числе уравнения переноса, посвящено множество работ. Для задач с анизотропным рассеянием построены несколько прие мов ускорения: линейные – с релаксацией, с вариационным функционалом, – по «выбранным направлениям»;

нелинейные – «квазидиффузионное», K–P-, – метод средних потоков (СП-метод) и др.

Ускорение по «выбранным направлениям», т. е. с переходом по итерациям на более редкие угловые сетки и соответствующей перенормировкой функции источника, по сравнению с простыми итерациями дает ускорение в среднем в 3–5 раз. Процесс замедляется к концу счета, когда ускорение уже не оказывает заметного влияния и устанавливается фактически процесс с простыми итерациями. Этот подход учитывает асимптотические особенности установления угловых профилей для оптически толстых и анизотропно рассеивающих сред.

Пусть 1 где A = Dz S, f = Dz F, y = Ay + f, – неоднородное уравнение с линейным оператором A и пусть y – решение – – этого уравнения:

y = A + f.

y Если y s – некоторое приближение к y, то положим – s+1/ = A ys + f y и очередное приближение определим по правилу y s+1 = y s + y s+1/2.

36 Глава 1. Одномерные плоские задачи Определение y s+1 по такому правилу есть линейное релаксационное ускорение.

Пусть i, i, i = 1, 2,..., – собственные функции и собственные значения – основной задачи. Наличие дискретного спектра установлено В. С. Владими ровым:


0 i 1, max i = max = 1, 0.

Представим y s = y + s, где отклонение s-приближения от точного решения s = cs i.

i i В простейшем случае = 1, когда y s+1 = y s + (y s+1/2 y s ), (1.29) находим, что s+1 = s + (A s s ), или cs+1 [1 (1 i )] cs = i cs, i i i где i 1 (1 i ).

Для сходимости итерационного процесса (1.29) необходимо и достаточно, чтобы |1 (1 i )| 1 для всех i = 1, 2,..., или 0, 1 min где min = min i. Оценку для max можно получить из анализа i. Легко видеть, что = min 0 переходит в = min 1 ;

(0, 1 ) и (1, 1 );

= max 1 и = max 1.

При 2 (при min 0) область изменения (0, 1 ) отражается в область изменения (1 + 1, 1 2) c 1 0, но 1 0.

Для предотвращения «раскачки» итераций при i 1 потребуем, чтобы min max, т. е. (1 ) (1 ), и отсюда находим оценку max.

1+ Такие значения max легко получить в процессе решения уравнения переноса простыми итерациями, поскольку ns max 1 max,, ns 1.1. Скалярная плоская задача где интегральная плотность ns = y s ds.

Релаксационные параметры (1, max ) обеспечивают устойчивую сходимость итераций с экстраполяцией («верхняя релаксация»). При этом | ln s | s (max )s или s (1 )s и s, т. е. число итераций, необходимых для достижения точности s, сокращается в раз, в лучшем случае почти в два раза, если max 2 (обычно max 2). Для итерационного метода с релаксацией и интегрированием по характеристикам типично монотонное установление гладких углового и пространственного профилей решения по итерациям.

В отличие от релаксационного итерационного метода Янга для задач с эллиптическими операторами для уравнения переноса не существует оптимального параметра, при котором число итераций резко сократилось бы по сравнению с простыми итерациями: = 1, 5 обеспечивает ускорение сходимости в полтора раза практически для всех физических расчетов, но при = max ускорение сходимости в max 2 раз.

Если известен почти весь спектр собственных чисел основной задачи i (0 i+1 i 1) и на s-итерации взять s = 1/(1 s ), как в методе Келлога, то на первой итерации c1 = 0, далее c2 = c2 = 0, на s-итерации для 1 1 всех k = 1 s cs = 0. Тем самым на первых s-итерациях постепенно гасятся k все первые s гармоник i с большими i, а затем простые итерации с = имеют хорошую сходимость в силу малости коэффициентов cs, поскольку k s + 1 имеет место cs+p = p cs, p для k 1, но k 1 1, и поэтому k k k cs+p cs. Однако задача нахождения спектра собственных значений для k k уравнения переноса труднее расчета решения самого уравнения.

Если на s-итерации берем большое, выходящее за пределы [1, max ], но такое, что 2 1 i+1 1 i (i неизвестно), то для коэффициентов при соответствующих i получаются следующие оценки:

при i = 1 i, cs+1 cs cs+1 cs при i i, i i i i т. е. коэффициенты cs сходятся для первых i собственных функций i и i расходятся для i + 1, i + 2,... Но при больших i сами i являются сильно колеблющимися функциями. В результате наступает «раскачка» решения и сходимость не имеет места. Если перейти к меньшим [1, max ], то «авост»

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

38 Глава 1. Одномерные плоские задачи Один из возможных приемов выбора в зависимости от итерации – – использование вариационного принципа, когда параметр определяется из условия минимума квадратичного функционала F[y] = (y, y) (2y, y) 2(f, y), где в качестве y выступает y s+1 (s ) = y s + s (y s+1/2 y s ), а выражения типа (u, v) – скалярные произведения. Из такого условия можно – выбрать и два параметра релаксации:

y s+1 (s, s ) = s y s + s y s+1/2.

В общем случае линейного оператора A оценку сходимости таких ите рационных процессов получить не удается. Но если A – самосопряженный – оператор, то устанавливается сходимость в смысле нормы | y s |2 = ( y s A [ y s ], y y s ) 0.

y y y В случае изотропного рассеяния уравнение переноса эквивалентно урав нению Пайерлса для интегральной плотности n(x) H E1 (|x x |) n(x ) dx + f (x), n(x) = где H S exp(t) f (x) = 0 E1 (|x x |) exp(x ) dx, E1 (x) = dt.

2 t x Обозначим H A y(x) E1 (|x x |) y(x ) dx, A – линейный самосопряженный оператор. Построим итерации так, что – y s+1/2 = n(s+1) A ns + f.

y s = ns, y s+1 = ns+1, Из условия min F [y s+1 ] получаем (n(s+2) n(s+1), n(s+1) ns ) s = 1, (n(s+1) ns ) 1.1. Скалярная плоская задача а для процесса с двумя параметрами (f, ns ) c (f, n(s+1) ) b a (ns A ns, ns ), s =, D (f, n(s+1) ) a (f, ns ) b b (ns A ns, ns+1 ), s =, D c (ns+1 A n(s+1), n(s+1) ), D a c b2.

При y s (x) 0 определитель D = 0, а потому при двух параметрах нельзя брать нулевое приближение y 0 (x) 0.

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

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

Для наглядности основные идеи СП-метода изложим на примере одно скоростной задачи d + (x, ) = (x, ) (x,, ) d + f (x, ), dx 1 (1.30) = 0, = 0.

0 H Решение основной задачи L = S + f ищем итерационным методом:

Ls+1 = S s + f, s+1 s+ = 0, = 0, 0 H где s выражается через s в зависимости от способа ускорения сходимости.

Для простых итераций s s. В основе СП-метода лежит перенормировка углового профиля решения на s-итерации по правилу ns+ (x) s+ s+ (x, ), 0, s (x, ) = n (x) (1.31) ns (x) s (x, ), 0, s n (x) где s± (x, ) – результат интегрирования уравнения (1.30) по характеристи – кам на s-итерации:

1 (x) (x) s+ s+ s s (x, ) d, n (x, ) d, n 40 Глава 1. Одномерные плоские задачи а ns+ (x) и ns (x) – аналоги средних потоков нисходящего и восходящего – излучения в положительном и отрицательном направлениях – являются реше – нием следующей вспомогательной задачи. Проинтегрируем уравнение (1.30) по отдельно от 1 до 0 и от 0 до 1 с некоторой весовой функцией V (x, ) и получим систему двух уравнений:

+ + (x) d n + a+ (x) n+ (x) + b+ (x) n (x) = f + (x), n+ = 0, dx (1.32) d n (x) + a (x) n+ (x) + b (x) n (x) = f (x), n = 0, dx H в которых коэффициентами являются «средние косинусы»:

1 + (x) V (x, ) (x) s+ V (x, ) s (x, ) d, (1.33) (x, ) d, а другие коэффициенты и правые части – интегральные величины по :

– 1 ds+ + a (x) d s+ V (x, ) (x, ) d + V (x, ) dx 0 1 V (x, ) d s+ (x,, ) (x,, ) d ;

0 1 b+ (x) s (x, ) (x,, ) d, V (x, ) d 0 a (x) V (x, ) d s+ (x, ) (x,, ) d ;

1 0 ds b (x) d s V (x, ) (x, ) d + V (x, ) dx 1 0 s (x, ) (x,, ) d ;

V (x, ) d 1 1 + f (x) f (x) (1.34) V (x, ) f (x, ) d, V (x, ) f (x, ) d.

n+ (x) n (x) Определение и как решения системы (1.32) с коэффици ентами (1.33), (1.34) и вычисление s (x, ) по формуле (1.31) называем вспомогательной («ускоряющей») задачей. Одна полная итерация СП-метода состоит из решения следующих задач:

1.1. Скалярная плоская задача 1) интегрирование уравнения (1.30), вычисление «средних косинусов»

(1.33), коэффициентов и источников (1.34);

2) решение вспомогательной задачи;

3) вычисление следующего приближения (1.31).

Таким образом, решение линейного уравнения переноса сводится к нелиней ному итерационному процессу.

При симметричных по индикатрисах рассеяния и свободных членах уравнение переноса (1.30) с помощью замены du 2 u(x, ) = (x, ) + (x, ), 2 = (x, ) (x, ) dx приводится к самосопряженной форме (1.35) L0 u = Su + f (x, ) с операторами d 2u L0 u Su 2 u(x, ) d + u(x, ), dx и однородными граничными условиями du du u = 0, = 0. (1.36) u+ dx dx H Построение итерационного процесса для решения задачи (1.35)–(1.36) СП-методом проводится аналогично. Приближение us+1 (x, ) определяется интегрированием уравнения L0 us+1 = S us+1 + f (x, ) при краевых условиях (1.36). Поскольку u(x, ) – симметричная функция, – достаточно положить us (x, ) = ns (x) us (x, ), где ns (x) – поток, усредненный по профилю us (x, ), – находится из обык – – новенного дифференциального уравнения 1 1 1 d 2 (ns us ) + ns us V d = 2 s s V d n u d + V f d dx 0 0 0 при граничных условиях dns dns = 0, = 0.

dx dx H Аналогично вспомогательная задача строится для одномерной сферической геометрии, естественно, с учетом локальных особенностей. Многочисленные методические исследования показали, что объем вычислений в СП-методе 42 Глава 1. Одномерные плоские задачи сокращается в десятки раз по сравнению с простыми итерациями. СП-метод, как и «квазидиффузия», основан на балансном принципе, который приводит к консервативной итерационной схеме.

Для задач, допускающих приведение к самосопряженной форме, решение задачи (1.35)–(1.36) сообщает минимум квадратичному функционалу G[u] = (L0 u, u) (Su, u) 2(u, f ) в пространстве N симметричных по функций u(x, ) со скалярным произведением 1 H H u(x, ) v(x, ) d = (u, v) = dx dx u v d.

0 0 Эти функции удовлетворяют условиям абсолютной непрерывности вдоль характеристик уравнения (1.35), условиям интегрируемости весьма общего характера и граничным условиям (1.36). Всякая последовательность, ми нимизирующая функционал G[u], сходится по среднеквадратичной норме с некоторым весом к решению уравнения (1.35). Выберем в качестве допустимых функций u(x, ) = n(x) (x, ), где (x, ) – некоторая фиксированная функция из пространства N. Чтобы – найти функцию n(x), сообщающую минимум функционалу G [u], построим уравнение для функционала G [n], считая заданной функцией. Положим dn G n, dx 1 d(n) (n) (n) (x,, ) d 2(n) f (n)2 + 2 d, dx 0 1 M [n(H)] (n) d N [n(0)] (n)2 d,.

x=H x= 0 Уравнения Эйлера для задачи с естественными краевыми условиями запишутся так:

d G dn Gn = 0, G dn + Mn |x=H = 0, G dn Nn |x=0 = 0. (1.37) dx dx dx dx Подставив в (1.37) соответствующие выражения, придем к уравнению, совпадающему с уравнениями СП-метода (1.32), если положить V =.

Краевые условия (1.37) также совпадают с условиями (1.32) СП-метода.

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

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

При решении «ускоряющей» системы уравнений методом прогонки пере ходим к разностной сети xk, проинтегрировав уравнения по x от xk до xk+ и заменяя n+ (x) и n (x) на [xk, xk+1 ] их средними значениями, т. е. для получения дискретного представления используем интегро-интерполяционный метод: + + + + + ck nk+1 + dk nk + bk nk+1 + bk n k = fk, a n+ + a n+ + hk n + tk n k = fk, k k+1 kk k+ + n k = 0, k = 0 (K 1), n0 0, где xk+ + + + + + + a+ a± dk = a+ a± (x) dx, + k+1 k k+1 k ck =, =, k k k k k xk xk+ + + b b± tk = b b± (x) dx, + k+1 k k+1 k hk =, =, k k k k k xk xk+ ± ± = ± (xk ), n± = n± (xk ).

f ± (x) dx, fk = k k k xk Прогонку вводим в виде n+ = 0, q0 = 0 ;

p0 = 0, n+ qk n, k =1K;

= pk + k k n = 0, K n = k + k n, k = (K 1) 0.

k k+ Коэффициенты схемы находятся по следующим формулам:

для k = 0 (K 1) + mk lk h k l k b + mk Vk = ck mk a lk, k k k pk+1 =, qk+1 =, k Vk Vk для k = (K 1) a pk+1 hk + a qk+ k k k k = k =,, mk mk где lk = b+ + qk d k, mk = tk + qk a, = fk a pk, + = fk pk d k.

+ k k k k k 44 Глава 1. Одномерные плоские задачи n Решение (0, ) (H, ) 2 Решение 2 1 0,75 1 0,5 1 2 0,25 1 2 Решение 0 8 H 1 0,5 2 6 0 0, Рис. 1.2. Изменение профилей плотности n по итерациям в СП-методе Сплошные линии n(H )) – интегрирование по характеристикам, штриховые – результат – – после «ускоряющей» задачи. Цифры у кривых – номер итерации – Из анализа численного материала установлено, что коэффициент + (x) имеет ярко выраженный минимум вблизи границы x = 0, а (x) – максимум – вблизи x = H. Такое поведение обусловлено локальными особенностями решения уравнения переноса вблизи границ слоя. Для повышения точности СП-метода разностная сеть xk для «ускоряющей» задачи вблизи границ выбирается более подробной. Алгоритм СП-метода обобщается на N СП-метод, когда интегрирование уравнения по проводится на N + отрезках [n, n+1 ] для 0 и N отрезках для 0. Обычно достаточно взять N, N + 3.

Типичный пример установления интенсивности и плотности излучения по итерациям в СП-методе приведен на рис. 1.2.

СП-метод по сравнению с другими методами (при одинаковой достигаемой точности решения = 103 106 ) имеет ряд преимуществ:

1) метод одинаково применим для любых индикатрис рассеяния, не требует их разложения в ряд по полиномам Лежандра и выделения какой-то части, как в «квазидиффузии»;

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

3) внешние краевые условия ускоряющей системы ставятся точно (не зависят от итераций) в отличие от «квазидиффузии»;

1.1. Скалярная плоская задача 4) СП-метод является балансным: уже после первых итераций выполня ется соотношение H H ns (x) dx n(x) dx ;

0 5) СП-метод слабо зависит от начального приближения. Обычно для установления профиля решения требуется 3–5 итераций и за 7–12 итераций получаются 3–4 установившиеся цифры.

6) Сходимость СП-метода в задачах с большой оптической толщиной ( 10) в десятки раз быстрее простых итераций (6 итераций с ускорением и примерно 400 итераций без ускорения при относительной точности расчета = 105 ).

1.1.5. Учет сильной анизотропии рассеяния. Аэрозольные атмосферы, облачные и морские среды характеризуются сильно вытянутыми индика трисами рассеяния, когда отношение ( = 0 )/( = 180 ) достигает нескольких порядков. В таких случаях решение уравнения переноса имеет две особенности: одна обусловлена -источником мононаправленного излучения, вторая – сильной вытянутостью индикатрисы в направлении малых углов – рассеяния.

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

(z, cos ) = c(z) (1 cos ) + (z) (z, cos ), (1.38) где (z, cos ) – сглаженная индикатриса, полученная из исходной графическим – или другим способом «обрезания» вытянутости в области малых углов.

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

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

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

Изложим алгоритмы учета сильной анизотропии рассеяния на примере решения задачи d = S (s s0 ), (1.39) + (z,, ) = B, = R;

dz 0 H q B (z) ds, R = ds.

+ Выделив нерассеянное излучение 0 в аналитическом виде, интенсивность многократно рассеянного излучения находим из задачи ( = 0 + ) d = 0, + = B + B 0, = R + R 0.

dz 0 H Однократное рассеяние 1 ( = 1 + N2 ) определяется аналитически интегрированием по характеристикам z zz S z + (z, + (z ) (z, cos 0 ) exp +, ) = + dz, 1 z S z + (z, 0, ) = exp (1.40) (z ) (z ) dz, 1 0 H z z S z (z,, ) = (z ) (z, cos 0 ) exp dz + | | 1 | | z H z H + q S 0 exp.

| | Функция 1 пропорциональна индикатрисе (z, cos 0 ), а потому анизо тропия рассеяния наиболее сильно выражена в нисходящем излучении;

+, особенно вблизи направления s0 = {0, 0 } с cos 0 = 1. Функция 1 не будет точно повторять угловой ход (z, cos ) из-за наличия двух множителей, зависящих от :

1.1. Скалярная плоская задача а) дробные множители 1/+ и 1/| | заметно увеличивают 1 при малых значениях + и ;

б) экспоненциальный множитель уменьшает 1 при малых + и.

При фиксированных значениях + азимутальный ход 1 совпадает с зависимостью (z, cos 0 ) от. Это обстоятельство часто используют для на хождения индикатрисы по однократному рассеянию пропущенного излучения в альмукантарате. Экстремум в модифицированном виде может сохраняться в функции многократного рассеяния N2, удовлетворяющей задаче dN = 0, + N2 = B N2 + B 1, N2 N2 = R N2 + R 1.

dz 0 H Проанализируем разные алгоритмы учета анизотропии.

АЛГОРИТМ 1. Запишем источник B 0 в виде суммы двух источников:

B 0 = B 0 + B0, z B 0 (z) (z) 0 ds = S exp (z) (z) (z, cos 0 ), z B0 (z) c(z) 0 (1 cos 0 ) ds = S exp (z)c(z) (1 cos 0 ).

Выделим в 1 = 1 + 1 гладкую часть 1 и вытянутую в направлении s сингулярную часть 1 как решения задач d 1 = 0, + 1 = B 0, 1 1 = R 0, dz 0 H d = 0, = + 1 = B0, 1 dz 0 H аналитически по формулам (1.40). При этом используем представление:

1 = 1 (z, 0, 0 ) (1 cos 0 ), z S z 1 (z, 0, 0 ) exp (t) c(t) dt.

0 Для однородного слоя, в котором c(z) = c = const, (z) = = const, z cz 1 = S exp.

0 По индукции можно показать следующее. Если для некоторого последо вательного приближения с номером Nn 2 имеем краевую задачу d Nn + Nn = B Nn + B n1, dz (1.41) N = 0, Nn = RNn + Rn1, n 0 H 48 Глава 1. Одномерные плоские задачи в которой функция рассеяния (n 1)-й кратности имеет представление n1 = n1 + n1, (z, 0, 0 ) (1 cos 0 ), B n1 = Bn1 + Bn1,, Bn1 = B n1 + (z) c(z) n1 + (z) (z) n1, (z, cos 0 ), q Bn1, = (z) c(z) n1, (1 cos 0 ), Rn1 = Rn1 + 0 n1, (H), n z cz n1, (z, 0, 0 ) = S exp (1.42), (n 1)!

0 то такие же представления справедливы и для приближения Nn+1. Таким образом, суммарная интенсивность в приближении n итераций n n (1 cos 0 ) + Nn+1.

(1.43) n + n n =1 n = При n получаем = (1 cos 0 ), = +, lim Nn+1 = 0, = n, = n, n n=1 n= – где – гладкая составляющая полной функции многократного рассеяния, – ее -анизотропная составляющая в направлении s0.



Pages:   || 2 | 3 | 4 | 5 |   ...   | 15 |
 





 
© 2013 www.libed.ru - «Бесплатная библиотека научно-практических конференций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.