CC BY
0
Вычислительные технологии, 2022, том 27, № 2, с. 4-18. © ФИЦ ИВТ, 2022 Computational Technologies, 2022, vol. 27, no. 2, pp. 4-18. © FRC ICT, 2022
ISSN 1560-7534 elSSN 2313-691X
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ
D01:10.25743/ICT.2022.27.2.002
Численное стохастическое моделирование рассеяния оптического излучения ледяными кристаллами нерегулярных случайных форм
Б. А. Каргин1'2'*, Е. Г. Каблукова2, Ц. Му1
1 Новосибирский государственный университет, 630090, Новосибирск, Россия 2Институт вычислительной математики и математической геофизики СО РАН, 630090, Новосибирск, Россия
*Контактный автор: Каргин Борис Александрович, e-mail: bkargin.osmfOmail.ru
Поступила 2 марта 2022 г., доработана принята в печать 10 марта 2022 г.
Перистая облачность является существенным источником неопределенности при построении численных моделей радиационного баланса Земли. В данной работе вычисляются индикатрисы рассеяния оптического излучения кристаллическими ледяными частицами нерегулярной формы на основе геометрической оптики с использованием метода трассировки лучей, что позволит дополнить оптические модели перистой облачности новыми данными. В качестве моделей ледяных частиц нерегулярной формы предлагается использовать многогранники, являющиеся выпуклыми оболочками случайного множества точек, распределенных по некоторому закону в ограниченном объеме. Для оценки оптических характеристик ледяных кристаллов использованы модели с гладкими и шероховатыми поверхностями. Показано, что свойства рассеяния излучения кристаллами в виде гексагональных призм, типичных для большинства теоретических и прикладных работ, и предложенными выпуклыми многогранниками имеют значительные различия. Предложенные модели дают новую интерпретацию индикатрис рассеяния ледяных частиц, наблюдаемых в атмосфере.
Ключевые слова: перистые облака, геометрическая оптика, трассировка лучей, индикатриса рассеяния, шероховатая поверхность, выпуклая оболочка.
Цитирование: Каргин Б.А., Каблукова Е.Г., My Ц. Численное стохастическое моделирование рассеяния оптического излучения ледяными кристаллами нерегулярных случайных форм. Вычислительные технологии. 2022; 27(2):4-18. D01:10.25743/ICT.2022.27.2.002.
Введение
Известно, что перистая облачность оказывает существенное влияние на радиационный баланс атмосферы и изменение альбедо климатической системы Земли (см. работы [1, 2] и библиографию в них). Разработка численной радиационной модели перистой облачности, которая позволит оценить потоки приходящего и уходящего видимого и ближнего инфракрасного излучения, требует прежде всего знания таких базовых оптических характеристик кристаллических облаков, как индикатрисы рассеяния и сечения ослабления излучения облачными кристаллическими ледяными частицами. Для расчета этих характеристик, в свою очередь, необходима детальная информация
о структуре микрофизических параметров облачности, таких как концентрация и распределение частиц по формам, размерам и пространственной ориентации (для частиц несферических форм).
Согласно экспериментальным данным большая часть ледяных облачных частиц имеет неправильную и сильно изменчивую форму [2]. В частности, в [3] показано, что только 3% ледяных кристаллов арктических облаков имеют идеально правильную форму (в виде гексагональных призм, пластинок, игл, дендритов) и до 90% частиц атмосферного льда в слоистых облаках могут иметь неправильную форму во всех температурных интервалах [4]. Эти исследования объясняют, почему при наблюдении перистых облаков гало и другие оптические проявления, характерные для правильных многогранных кристаллов льда, наблюдаются довольно редко. Большинство численных расчетов оптических характеристик ледяных кристаллов выполнены для частиц правильной формы [1, 2, 5], но в работе [6] приведены результаты вычисления индикатрис рассеяния для частиц с шероховатой поверхностью или частиц случайной формы, со случайно варьирующимся углом между гранями, В работе [7] представлены индикатрисы рассеяния для моделей кристаллов, полученных из частиц правильной формы путем отсекания от них частей случайно наклоненными плоскостями,
В настоящей работе предложена новая модель кристалла случайной формы в виде выпуклого многогранника с треугольными гранями, геометрические параметры которого подчиняются заданным распределениям вероятностей. Для таких моделей ледяных частиц, крупных в сравнении с длиной волны, вычисляются индикатрисы рассеяния и сечения ослабления излучения с применением геометрической оптики,
В первой части работы кратко описан алгоритм компьютерного построения случайных реализаций трехмерных выпуклых тел. Более подробно алгоритм изложен в работе [8]. Алгоритмы компьютерного моделирования случайных реализаций таких тел являются предметом исследований и приложений таких разделов вычислительной информатики и теории вероятностей, как вычислительная геометрия [9] и геометрические вероятности [10]. Кроме того, такие алгоритмы применяются во множестве задач физики, астрономии, биологии, кристаллографии, петрографии, компьютерной томографии, компьютерной анимации и т.д. Во второй части работы описаны алгоритмы метода Монте-Карло для расчета первичных оптических характеристик ледяных кристаллов с гладкой и шероховатой поверхностью: индикатрис рассеяния и сечений ослабления излучения, необходимых при моделировании переноса оптической радиации в кристаллической облачности. Приведены примеры вычисленных индикатрис рассеяния для отдельных реализаций выпуклых случайных тел с гладкими и шероховатыми гранями, а также представлен вид индикатрисы рассеяния, осредненной по множеству таких реализаций,
1. Модели трехмерных выпуклых частиц несферической формы
В настоящем исследовании моделями кристаллических частиц являются выпуклые многогранники, а именно выпуклые оболочки, содержащие набор заданных или случайно сгенерированных точек. Для их построения использовались инкрементальный алгоритм и алгоритм направленных ребер [11]. Начальный многогранник является тетраэдром, вершины которого — случайные точки, распределенные в заданном объеме по некоторому закону распределения. Дальнейшее построение многогранника происходит путем моделирования следующей случайной точки и определения, является она внутренней
Рис. 1. Примеры случайных выпуклых многогранников, построенных программой ConvexHull Fig. 1. Examples of the random convex polyhedra constructed by the program ConvexHull
или внешней точкой для построенного многогранника. Если точка лежит за пределами многогранника, она становится повой его вершиной и т. д. Детали построения трехмерных выпуклых тел, являющихся многогранниками произвольных неправильных форм, содержатся в [8, 11]. На основе описанного в [8, 11] алгоритма разработана программа ConvexHull, которая в результате выдаст заданное или случайное число граней конечной выпуклой оболочки, координаты вершин, "правильные" очереди вершин для каждой грани, которые удобно использовать при исследовании свойств рассеяния излучения песферическими частицами методом трассировки лучей. Правильная очередь вершин задаст направления ребер па грани таким образом, что их векторное произведение соотвстствовуст правилу правой руки и определяет направление внешней нормали. Программа написана па языке С--, в пей реализована функция визуализации многогранников с помощью библиотеки OpenGL.
На рис. 1 приведены примеры построенных случайных выпуклых многогранников, для которых случайные точки, используемые в построении многогранников, моделируются равномерно в кубе. Выбором подходящей плотности распределения координат случайных точек, а также путем растяжения или сжатия полученных тел вдоль выбранных направлений можно получать реализации многогранников с заданными свойствами. Соотношение числа вершин и граней моделируемых многогранников однозначно определяется формулой Эйлера
Nv — Ne + Nf = 2,
а для многогранника с треугольными гранями верна формула 2Ny — Np = 4. Здесь Nv — число вершин многогранника, Ne — число ребер и Np — число граней.
2. Трассировка лучей
Для вычисления индикатрис рассеяния и сечений ослабления излучения кристаллическими частицами, крупными относительно длины волны излучения, используется приближение геометрической оптики, в частности алгоритм трассировки лучей. В данном алгоритме траектория оптического излучения представляется в виде ломаной липни, вершины которой — это точки столкновения с гранями кристалла либо с псодпородпос-тями рассеивающей среды. Вначале введем несколько известных понятий и утверждений, важных для выполнения процедуры трассировки световых лучей, попадающих па поверхность кристалла, и дальнейших операций, предназначенных для расчета сечений и индикатрис рассеяния кристаллической среды. Положение кристалла в пространстве удобно задавать углами Эйлера (рис. 2) [12]. Эти углы, определяющие поворот абсолютно твердого тела в трехмерном евклидовом пространстве, позволяют описать любое угловое преобразование системы в глобальной (лабораторной) системе координат. Обозначим оси лабораторной системы координат через (x,y,z), оси преобразованной системы координат — как (X,Y,Z). Пересечение координатных плоскостей ху и XY
I У
Рис. 2. К определению углов Эйлера Fig. 2. On the definition of Euler angles
называется линией узлов Ж. Введем обозначения: a Е Да = [0,2п] — угол между осью ж и линией узлов, ß Е Дß = [0,^] — угол между осями z и Z, 7 Е Д7 = [0, 2п] — угол между осью X и линией узлов.
Все ориентации кристаллов относительно фиксированной системы координат или вращения падающих лучей относительно фиксированного положения кристалла можно однозначно задать тройкой углов [a, ß, 7], определяющих поворот оси кристалла или луча. Совокупность индикатрис рассеяния для кристаллической частицы несферической формы представляет собой трехмерную матрицу, в которой для каждого направления падающего излучения, заданного углами Эйлера, задаются объемные индикатрисы рассеяния (в декартовой системе координат — трехмерные, а в полярной — двумерные) — непрерывные или дискретные вектор-функции по угловым переменным, описывающие плотности распределения случайных направлений движения световых квантов после рассеяния. В ряде публикаций по фотометрии в качестве объемных индикатрис рассеяния используется термин "тело яркости". Для вычисления индикатрис рассеяния излучения случайно ориентированными частицами требуется разыгрывать или задавать различные ориентации частиц. Вычисление индикатрис рассеяния частицами со случайной ориентацией можно осуществлять путем изотропного или иного вращения кристаллических частиц внутри некоторой заданной опорной сферы или вращения падающих на опорную сферу световых лучей при фиксированном положении частицы. В данной работе используется второй вариант.
Через щ,а обозначим показатель преломления льда относительно воздуха. Обозначим через n вектор внешней нормали к поверхности грани кристалла в точке падения светового кванта на эту грань, выраженный в глобальной (лабораторной) системе координат. Предполагаем, что кристалл является полностью прозрачным для световой волны, т. е. свет внутри кристалла не поглощается. В случае необходимости с помощью простейшей техники весового моделирования учет поглощения не представляет сложности. Для описания воздействия плоской грани кристалла на излучение примем модель, в которой отражение и преломление света гранью происходят по законам лучевой оптики (законам Снеллиуса и Френеля).
В таком случае при попадании луча, имеющего направление на поверхность грани с вероятностью R(u, n) происходит зеркальное отражение в направлении шreß = ш — 2(ш, n)n, а с вероятноетью 1 — R(u, n) происходит преломление в направлении urefr = vu — Dn. Здесь R(u, n) коэффициент отражения Френеля, который для удобства вычислений запишем в следующей преобразованной форме:
ш (И-В )2(А2В2+С2)
П)= (И|+ В)2(\А\В+С)2 ' ^ =("' П)' В = <
, А < 0,
^К«-1+А2, А > - 1/п2^ к 0, 0 < А < ^1- 1/<а.
С =1 - А2, Я = А - Sgn(A)B, И =[1/Пг,а, А < 0,
6 1 ; 1 П^а, А > 0.
Далее полезно обсудить некоторые простые, но весьма важные детали, связанные с вычислением случайных оценок целевых функционалов (сечений и индикатрис рассеяния) на множестве моделируемых траекторий световых квантов, испытавших взаимодействие с кристаллическими частицами, а также удобную для дальнейшего использования систему хранения полученной численной информации. Выполним дискретизацию пространства углов Эйлера, т.е. разобьем интервалы Да, Д^ и Д7 изменения углов а, ¡3 и 7 сответственно на па, п^ и п1 равных частей размером 5а, и 61 и построим сетку узлов:
[аА = {г5а,г =1, 2,...,па], Щ } = ЦЬц,] = 1, 2,...,щ}, {7к} = [Ы1 ,к = 1, 2,...,п1}.
Результаты моделирования будем регистрировать и хранить в виде некоторой "трехмерной блочной прямоугольной матрицы" С размером па х п^ х п-у. Номера блоков этой матрицы для суммирования случайных оценок будем определять по значениям углов а., ¡3 и 7, а именно, табор индексов г, А к реализуется при условии аг е [г6а,г = 1,2,...,па}, & е А = 1,2,...,щ}, е [к81 ,к = 1,2,...,п1}. "Элементом" каждого блока матрицы С является своя двумерная прямоугольная матрица д размером х п^. Здесь через е [0,^] и р е [0, 2^] обозначены соответственно широтный и азимутальный углы единичного вектора р (направление движения светового кванта после рассеяния) в полярной системе координат, связанной с исходной лабораторной декартовой системой координат ж, у, г, а именно $ — угол между вектором р и осью Ог, <р — угол между проекцией вектора р на плоскость хОу и осью Ох, отсчитываемый по часовой стрелке. Через и п^ обозначено количество интервалов, на которые разбиты области изменения углов § и <р.
Алгоритм трассировки лучей. Процедура численного оценивания па х пр х щ матриц (по числу блоков матриц О) тип а д па ансамбль моделируемых траекторий рассеянных кристаллами квантов света состоит в следующем,
1, Программа СопуехНиП в зависимости от постановки задачи выдает координаты вершин заданного или случайного выпуклого многогранника. Размеры частиц контролируются вычислением так называемого форм-фактора /, равного отношению длины главной оси многогранника к среднему диаметру основания. Случайное значение параметра f выбирается из заданного распределения вероятностей, входящего в число исходных данных, характеризующих микрофизическую структуру кристаллического облака. Размеры полученного случайного многогранника путем сжатия или растяжения главной оси подгоняются под значение параметра f.
2, Построенный кристалл помещается в условную опорную сферу таким образом, чтобы главная ось кристалла лежала па оси Ог фиксированной лабораторной системы координат, а ее средняя точка совпадала с началом координат. По координатам вершин многогранника можно выписать уравнения плоскостей, в которых лежат грани, а также угловые координаты векторов нормалей к этим плоскостям в системе координат (х, у, г).
3. Последовательно, соответственно узлам сетки = {г5а,г = 1, 2,... ,па}, {^} = ЪЧз ,3 = 1,2,...,пр}, {7^} = {к5у ,к = 1, 2,... ,п1}, строятся векторы направления падающего излучения ш (обозначим его
4. В круге с центром в начале координат и вектором нормали равномерно выбираем Да:,^ точек, го которых оказались внутри кристалла. Моделирование случайных точек, равномерно распределенных в круге, удобно проводить методом исключения [13, 14, с. 157]. В случае если точка оказалась внутри кристалла, из данной точки строится прямая линия I с направляющим век тором (рис. 3).
5. Проверяется, пересекает ли прямая I грани кристалла. Для этого вычисляем координаты точки Рг пересечения прямой I и плоскости, содержащей г-ю грань кристалла. Обозначим направляющие векторы двух соседних ребер г-й грани кристалла, выходящие из одной вершины, через а вектор, направленный из этой вершины в точку Ри через к, тогда, если
(С,V) < (Ск) (С^) < (^к) 1СЫ 1С 1|к|, 1СМ|к|,
то точка Рг — внутренняя для рассматриваемого угла. Если аналогичные неравенства выполнены для всех вершин г-й грани, то точка Рг лежит на г-й грани кристалла, т. е. является точкой пересечения светового луча и кристалла. Поскольку многогранник выпуклый, такая точка, если она существует, будет единственной. В счетчик вы-
числяющий количество квантов, попавших па поверхность кристалла, добавляется 1. Далее полагаем, что световой луч с направлением — ш попадает на кристалл извне, и продолжаем моделирование траектории луча, учитывая законы отражения и преломления.
6. При аналоговом моделировании в точке Pi осуществляется случайный выбор события, которое в соответствии с выбранной моделью воздействия поверхности кристалла на световые лучи может произойти. Обозначим через п вектор внешней нормали к поверхности кристалла в точке Р^ в лабораторной системе координат. В случае зеркального отражения, реализующегося с вероятностью К(—ш, п), световой квант движется в направлении рге^ = 2(ш, п)п — ш. Для блока г, к, матрицы С по координатам вектора рге/ь выраженным в переменных и определяется ячейка матрицы с/, в которую попадает это направление. В соответствующий счетчик суммируется 1. Если ре-
Рис. .3. К алгоритму моделирования облучения кристалла световым потоком в заданном направлении (зеленые точки лучи, падающие на частицу)
Fig. .3. То the modeling algorithm of the crystal irradiation with a light flux in a given direction (green dots are rays incident on a particle)
ализуется преломление, то световой квант проникает внутрь кристалла в направлении Рге^г = Оп - уш и движется прямолинейно до пересечения с соответствующей гранью. Внутри кристалла, в соответствии с принятой моделью, квант может испытать несколько отражений (в том числе и полное внутреннее отражение) от граней и затем после преломления вылететь наружу. По координатам вектора рГерг определяются индексы матрицы д блок а г, А К матрицы С, куда и суммируется 1,
В результате моделирования Ма,рг/ случайных траекторий весь совокупный вклад в искомый функционал, равный ^ 7, будет распределен в соответствии со случайными реализациями вектора р от каждой случайной траектории по х п^ клеткам матрицы д, хранящейся в блоке С с индексами г, А соответствующими индексам узла аг,
3. Модель рассеяния световых волн кристаллами с шероховатыми поверхностями
До сих пор речь шла о рассеянии световых лучей идеальными плоскими границами раздела двух сред. Однако многочисленные лабораторные и авиационные наблюдения показывают, что дня реальных атмосферных частиц и кристаллов льда характерны не идеально плоские, а шероховатые поверхности (см, рис, 4 из работы |15|),
Построение адекватной шероховатой поверхности является нетривиальной задачей в первую очередь из-за многообразия форм шероховатостей, а также из-за отсутствия падежных экспериментальных статистических данных о таких поверхностях. Рассмотрим далее две модели шероховатых поверхностей, заимствованные из области фотометрии, акустики и статистической структуры ветрового волнения морской поверхности, укладывающиеся в парадигму прямого моделирования трассировки оптических лучей.
Первая модель основана па представлении поверхностей ледяных частиц как абсолютно матовых. Такие поверхности отражают свет равномерно во всех направлениях, а отражение света от таких поверхностей называют диффузным (идеально рассеивающим), Отражение или рассеяние света такими поверхностями подчиняется закону Ламберта, который в вероятностной модели описывается выражением Р(р) = 2р, где р — косинус угла между направлением отраженного луча и нормалью к поверхности, а Р(р) — плотность распределения случайной величины р. Характерной особенностью такой модели является отсутствие зависимости случайного направления отраженного луча от направления падения луча на поверхность. Случайное значение р вычисляется по формуле р = у/^, оде £ — случайное число, распределенное равномерно в интер-
Рие. 4. Микроскопические изображения поверхностей .ледяных кристаллов Fig. 4. Microscopic images of ice crystal surfaces
вале (0,1) Азимутальный угол вектора отражения <р разыгрывается из равномерного распределения в интервале [0, 2^], Особенностью этой модели является то, что в ней поверхность предполагается плоской в отличие от реальных поверхностей. Поэтому для приближенного учета таких явлений, как переотражение и затенение одних участков поверхности другими в случае, когда случайные значения ^ близки к нулю, можно использовать простую процедуру исключения. При ^ < £ это значение ^ отбрасывается и разыгрывается новая случайная величина ^.Параметр £ зависит от свойств случайно неоднородной поверхности и подлежит предварительной оценке.
Этот же недостаток свойствен и второй модели шероховатости, применяемой для описания взволнованной поверхности океана в так называемом фацетном приближении (см., например, [16]), В этой модели поверхность представляется как совокупность случайно ориентированных микроплощадок (фацетов) с центрами, лежащими на гладкой поверхности, а случайные нормали к площадкам (уклоны пх и пу в локальной системе координат с осью OZ.J направленной вдоль перпендикуляра к поверхности) распределены по нормальному закону с плотностью
Р (пх,пу) = 1 ехЛ — 1 [ 2 +( х' у 2жахау ) 2 \ах) \&у)
где ах и ау — дисперсии уклонов. Сложность применения этой модели заключается в отсутствии каких-либо обоснованных данных о дисперсиях ах и ау. В расчетах их приходится выбирать по косвенным данным. Во второй модели, в сравнении с первой, с помощью приема, изложенного в [17], можно приближенно учесть вышеназванные эффекты переотражения и затенения без оценки параметра е. В работе использовались обе схемы для оценки влияния случайной неоднородности поверхностей кристаллов на их рассеивающие характеристики.
2
4. Вычисление геометрического поперечного сечения рассеяния световых волн кристаллическими частицами методом Монте-Карло
Для моделирования переноса лучистой энергии Солнца в кристаллических облаках наряду с индикатрисами рассеяния необходимо знать сечения рассеяния , г). Зависимость сечения рассеяния от вектора направления движения светового кванта связана с тем очевидным фактом, что кристаллическая среда является анизотропной вследствие отсутствия круговой симметрии кристаллов. Зависимость сечения рассеяния от радиуса-вектора г в глобальной системе координат х, у, г означает, что среда может быть неоднородной по пространственным координатам,
В данной работе процесс переноса не рассматривается, поэтому зависимость сече-г
чей, представленная в разд. 2, позволяет оценивать ) одновременно с расчетом
индикатрис рассеяния, С этой целью для каждого построенного многогранника вычисляется величина N1^ которая затем осредпяется по распределению всех форм многогранников. Эта величина умножается на пг2. Полученное значение, которое обозначим через является средней площадью проекций частиц па плоскость, перпендикулярную направлению ПРИ а = @ = Р], 1 = кинетической теории пе-
реноса величину Si^^ называют микроскопическим сечением рассеяния. После выполнения этой процедуры должна быть учтена дифракционная составляющая рассеяния. По теореме оптической экетинкции [1] дифракционная составляющая равна лучевой Si,j,k- В итоге имеем ) = 2Si>j>k при a = a^ ft = 7 = 7fc, В результате получа-
ем таблицу размером [na], в каждой ячейке которой храпим значения Т,а(шафу1) для всего набора узлов {ai} = {iSa,i = 1, 2,...,na}, [ftj} = {j5p ,j = 1, 2,...,np}, ilk} = {k51, k =1, 2,... ,n1}.
Заметим, что макроскопическое сечение рассеяния ), являющееся наряду
с индикатрисами рассеяния основной характеристикой среды и входящее в кинетическое уравнение Больцмана переноса частиц, получается умножением величины ) на объемную концентрацию частиц Nc в рассматриваемой среде,
5. Результаты некоторых численных экспериментов
Вычисления индикатрис рассеяния проводились для частиц случайной формы с гладкими и шероховатыми поверхностями. Для всех частиц был задан показатель преломления льда относительно воздуха Пг,а = 1.3078, Количество ячеек для регистрации рассеянных квантов по каждому углу па = пр = 180, п7 = 90, Исследовались характеристики случайно ориентированных частиц, поэтому значения матриц д осреднялись по направлениям входящего излучения. Распределение азимутального угла рассеяния <р в данном случае равномерно. Для оценки индикатриссы рассеяни отдельной частицы использовалось 108 траекторий фотонов. Относительная погрешность вычисления индикатрис рассеяния в этом случае не превосходила 3 % и была наибольшей для углов рассеяния 178-180°, Координаты всех точек, по которым строились исходные выпуклые многогранники, генерировались равномерно в интервале [-1,1], Проводились также вычисления для деформированных частиц, которые получались на основе исходных выпуклых многогранников путем растягивания или сжатия по определенным направлениям, На рис, 5 приведены результаты вычислений индикатрис рассеяния выпуклыми многогранниками с гладкими и шероховатыми поверхностями.
Из рис, 5 видно, что в индикатрисах рассеяния кристаллическими частицами случайной формы с гладкими гранями наблюдается несколько небольших пиков в интервале [0, 40°], Полученные значения индикатрис рассеяния в направлениях вперед и назад значительно меньше значений индикатрисы рассеяния правильной гексагональной призмы (рис, 6, о). Это объясняется так называемым явлением дельта-функции пропускания [1, 2, 18], вызванным двукратным преломлением (строго вперед), которое наблюдается, когда в качестве моделей берутся частицы с параллельными гранями.
Сходство между индикатрисами многогранника и его растянутого варианта наблюдается редко. Из расчетов видно, что многогранники, вытянутые вдоль главной оси,
пропускают большее количество энергии в направлении вперед (для углов рассеяния °
связаны с нерегулярной формой кристалла. Этот вывод можно сделать из того, что при увеличении числа моделируемых квантов излучения и уменьшении погрешности расчета форма индикатрисы рассеяния не изменяется.
На рис, 6, б представлены осредпеппые индикатрисы рассеяния для 10 000 реализаций выпуклых многогранников с гладкими и шероховатыми поверхностями. При оценке индикатрис рассеяния для каждой реализации частицы строились 10 000 траекторий фотонов. Число случайных точек N, по которым строился выпуклый многогранник,
100
10
0.1
0.01
Roughened
Smooth
100 Roughened
100
0 20 40 60 80 100 120 140 160 180 Scattering angle (degrees)
в
0.01
0 20 40 60 80 100 120 140 160 180 Scattering angle (degrees)
0.001
О 20 40 60 80 100 120 140 160 180 Scattering angle (degrees)
1000
100
10
0.1
0.01
(bughened I
[A oth
0 20 40 60 80 100 120 140 160 180 Scattering angle (degrees)
a
г
Рис. 5. Индикатрисы рассеяния PhF для двух реализаций случайных выпуклых многогранников с гладкой (Smooth) и шероховатой (Roughened) поверхностью при Ny = 12, Np = 20 (а) и Ny = 25, Np = 46 (в) и для тех же многогранников, вытянутых вдоль оси Oz с коэффициентом растяжения f = 10 (5 и г соответственно)
Fig. 5. The scattering phase functions PhF for two realizations of a random convex polyhedron with smooth (Smooth) and rough (Roughened) surfaces for Ny = 12, Np = 20 (a) and Ny = 25, Np = 46 (e) and for the same polyhedron stretched along the Oz axis with a stretching factor f =10 (6 and г respectively)
моделировалось равномерно в интервале значений 5 < N < 30. Число вершин моделируемых многогранников находилось в интервале 4 < Ху < 24.
Аналогичные выводы о форме индикатрис рассеяния для частиц неправильной формы представлены в работе [7], в которой методом физической оптики численно исследуется сигнал обратного рассеяния при дистанционном зондировании. Но выводы о поведении индикатрис рассеяния в окрестности направления рассеяния назад в наших расчетах расширяются на весь диапазон углов рассеяния [0,180°]: индикатрисы рассеяния, вычисленные для отдельной крупной частицы неправильной формы, представляют собой сложные функции с множеством пиков и провалов, различных по положению и ширине. Однако в индикатрисах рассеяния, осредненных по форме частиц, эти пики и провалы сглаживаются. Согласно расчетам, представленным на рис. 5 и 6, 6. в слу-
10 000
20 40 60 80 100 120 140 160 180 Scattering angle (degrees)
20 40 60 80 100 120 140 160 180 Scattering angle (degrees)
Рис. 6. Индикатрисы рассеяния PhF для правильной гексагональной призмы с гладкой (Smooth) и шероховатой (Roughened) поверхностями (а), индикатрисы рассеяний, осреднен-ные по 10 ООО реализациям выпуклых многогранников с гладкими (Smooth) и с шероховатыми поверхностями (Roughened) (б)
Fig. 6. The scattering phase functions (PhF) for a regular hexagonal prism with smooth (Smooth) and rough (Roughened) surfaces (a), the scattering phase functions averaged over 10 000 realizations of convex polyhedra with smooth (Smooth) and rough surfaces (Roughened) (б)
чае хаотической ориентации частиц поведение индикатрис рассеяния для кристаллов с шероховатыми поверхностями и кристаллов нерегулярной формы с гладкими поверхностями аналогично. Этот вывод совпадает с результатами из статьи [6].
Заключение
На основе прямого моделирования методом Монте-Карло представлен алгоритм трассировки световых лучей для исследования характеристик рассеяния оптического излучения ледяными частицами крупного размера в перистых облаках. Модели таких частиц разработаны на основе алгоритма построения выпуклых оболочек. Кратко описаны методы построения моделей многогранников с шероховатыми поверхностями.
Предложен эффективный способ определения пересечения грани кристалла и прямой при трассировке лучей. Приведены сравнительные данные нескольких численных экспериментов по расчету индикатрис рассеяния выпуклыми многогранниками с гладкими и шероховатыми поверхностями. Численные исследования показывают, что шероховатость поверхности играет весьма существенную роль в формировании рассеивающих свойств кристаллических частиц.
Рассмотренный алгоритм трассировки лучей прямым методом Монте-Карло хорошо распараллеливается не только по отдельным траекториям, но и посредством разделения на независимые расчеты по узловым значениям сетки направлений вращения лучей а, Д 7. Это весьма полезное свойство позволяет решать задачу на мультипроцессорных вычислительных кластерах. Тривиальные оценки объема вычислений даже при довольно грубом разбиении углового пространства по 1° по каждому из углов а, Д 7 показывают, что при моделировании всего ста случайных реализаций выпуклых оболочек на один узел сетки углов а, Д 7 и ста траекторий квантов на одну такую реализацию
для расчетов матрицы G с вложенными матрицами g потребуется моделировать порядка 1014 траекторий световых квантов.
Нужно заметить, что для некоторых приложений такой объем информации может быть сильно избыточным. Например, для задач наземного и аэрокосмического лазерного зондирования перистой облачности [7, 19-21], в которых интерес представляют оценки отраженного в направлении назад лазерного сигнала, В этом случае можно ограничиться существенно меньшим объемом вычислений и информации о рассеивающих свойствах кристаллов. При решении задач радиационного переноса в кристаллической облачности [1, 2, 5, 22] такая информация нужна в полном объеме и с приемлемой точностью.
Благодарности. Работа выполнена при финансовой поддержке CSC (China Scholarship Council) и в рамках проекта госзадания ИВМиМГ СО РАН № 0251-2021-0002,
Список литературы
fl] Liou K.N., Yang P. Light scattering by ice crystals: Fundamentals and applications. Cambridge University Press; 2016: 461.
[2] Mishchenko M.I., Travis L.D., Lacis A.A. Scattering, absorption, and emission of light by small particles. Cambridge University Press; 2002: 560.
[3] Korolev A.V., Isaac G.A., Hallett J. Ice particle habits in Arctic clouds. Geophysical Research Letters. 1999; 26:1299-1302.
[4] Korolev A.V., Isaac G.A., Hallett J. Ice particle habits in stratiform clouds. Quarterly Journal of the Royal Meteorological Society. 2000; (126):2873-2902.
[5] Mishchenko M.I., Hovenier J.W., Travis L.D. Light scattering by nonspherical particles: Theory, measurements and geophysical applications. San Diego: Academic Press; 1999: 721.
[6] Liu C., Panetta R.L., Yang P. The effective equivalence of geometric irregularity and surface roughness in determining particle single-scattering properties. Optics Express. 2014; 22(19):23620-23627.
[7] Шишко В.А. Исследование оптических свойств атмосферных ледяных кристаллов неправильной формы. Канд. диссертация. Томск: ИОА СО РАН; 2020: 143.
[8] My Ц., Каргин Б.А., Каблукова Е.Г. Компьютерное построение трехмерных выпуклых тел произвольных форм. Вычислительные технологии. 2022; 27(2):") I (¡1. D01:10.25743/ICT.2022.27.2.005.
[9] Preparata F.P., Shamos M.I. Computational geometry: An introduction. Springer Science Business Media; 1985: 413.
[10] Сантало Л. Интегральная геометрия и геометрические вероятности. Наука; 1983: 360.
[11] De Berg М., van Krefeld М., Overmars М., Schwarzkopf О. Computational geometry algorithms and applications. 3rd rev. ed. Springer-Verlag; 2008: 386.
[12] Журавлев В.Ф. Основы теоретической механики. М.: Физматлит; 2008.
[13] Михайлов Г.А., Войтишек А.В. Численное статистическое моделирование. Методы Монте-Карло. Издательский центр "Академия"; 2006.
[14] Марчук Г.И., Михайлов Г.А., Назаралиев М.А., Дарбинян Р.А., Каргин Б.А., Елепов B.C. Метод Монте-Карло в атмосферной оптике. Новосибирск: Наука; 1976: 283.
[15] Liu С., Panetta R., Yang P. The effects of surface roughness on the scattering properties of hexagonal columns with sizes from the Ravleigh to the geometric optics regimes. Journal of Quantitative Spectroscopy and Radiative Transfer. 2013; (129) :169—185.
[16] Пирсон В., Ньюман Дж., Джеймс Р. Ветровые волны. М.-Л.: Гидрометеоиздат; 1962: 42-124.
[17] Rakimgulov К.В., Ukhinov S.A. Local estimates in Monte Carlo method for the ocean atmosphere system with a random interface. Russian Journal of Numerical Analysis and Mathematical Modelling. 1994; 9(6):547-564.
[18] Takano Y., Liou K.N. Solar radiative transfer in cirrus clouds. Part I: Single-scattering and optical properties of hexagonal ice crystals. Journal of Atmospheric Sciences. 1989; 46(1):3—19.
[19] Коношонкин А., Кустова H., Боровой А. Рассеяние света на гексагональных ледяных кристаллах перистых облаков. Саарбрюкен: LAP LAMBERT Academic Publishing; 2013: 156.
[20] Коношонкин А.В. Рассеяние света на атмосферных ледяных кристаллах при лазерном зондировании. Докт. диссертация. Томск: ИОА СО РАН; 2017: 283.
[21] Konoshonkin A.V., Kustova N.V., Shishko V.A., Borovoi A.G. The technique for solving the problem of light backscattering by ice crystals of cirrus clouds by the physical optics method for a lidar with zenith scanning. Atmospheric and Oceanic Optics. 2016; 29(03):252-262.
[22] Macke A., Johannes M., Ehrhard R. Single scattering properties of atmospheric ice crystals. Journal of the Atmospheric Sciences. 1996; 53(19):2813-2825.
Вычислительные технологии, 2022, том 27, № 2, с. 4-18. © ФИЦ ИВТ, 2022 ISSN 1560-7534
Computational Technologies, 2022, vol. 27, no. 2, pp. 4-18. © FRC ICT, 2022 elSSN 2313-691X
MATHEMATICAL MODELLING
DOI: 10.25743/ICT.2022.27.2.002 Numerical stochastic simulation of optical radiation scattering by ice crystals of irregular random shapes
Kargin, Boris A.1'2'*, Kablukova, Evgeniya G.2, Mu, Quan1
Novosibirsk State University, 630090, Novosibirsk, Russia
2
Russia
*
Received March 2, 2022, accepted March 10, 2022.
Abstract
Prom numerous publications it is currently well known that cirrus clouds have a significant impact on the radiation balance of the atmosphere and the albedo change of the Earth's climate system. The development of a numerical radiative model of cirrus cloudiness, which allows estimating incoming and outgoing visible and near-infrared radiation fluxes, requires, first of all, knowledge of such basic optical characteristics of cirrus as scattering phase functions and attenuation cross sections of radiation by ice particles. In this paper, we propose a new model of a randomly shaped crystal in the form of a convex polyhedron with triangular faces, the geometric parameters of which obey given probability distributions. For such models of ice particles, which are large in comparison with the wavelength, the scattering phase functions and the radiation attenuation cross sections
are calculated using geometric optics. An efficient method for determining the intersection of a crystal face and a straight line in the ray tracing method is proposed. Models with smooth and rough surfaces are used to evaluate the optical characteristics of the ice crystals. Comparative data are presented for several numerical experiments which calculate the scattering phase functions by convex polvhedra with smooth and rough surfaces. The calculation results show that several small peaks are observed in the scattering phase functions of crystals with irregular shape and smooth faces. However, in the scattering phase functions averaged over the particle shapes and orientations, these peaks and dips are smoothed out. A similar behavior is demonstrated by the scattering phase functions of particles with rough surfaces. It is shown that the properties of radiation scattering by crystals in the form of hexagonal prisms, which are typical for most theoretical and applied studies, and by the proposed convex polvhedra have significant differences. The proposed models provide a new interpretation of the observed scattering phase functions of ice particles in the atmosphere.
Keywords: cirrus clouds, geometric optics, ray tracing, scattering phase function, rough surface, convex hull.
Citation: Kargin B.A., Kablukova E.G., Mu Q. Numerical stochastic simulation of optical radiation scattering by ice crystals of irregular random shapes. Computational Technologies. 2022; 27(2):4—18. D01:10.25743/ICT.2022.27.2.002. (In Russ.)
Acknowledgements. The work is supported by China Scholarship Council and the State Contract with ICM&MG SB RAS No. 0251-2021-0002.
References
1. Liou K.N., Yang P. Light scattering by ice crystals: Fundamentals and applications. Cambridge University Press; 2016: 461.
2. Mishchenko M.I., Travis L.D., Lacis A.A. Scattering, absorption, and emission of light by small particles. Cambridge University Press; 2002: 560.
3. Korolev A.V., Isaac G.A., Hallett J. Ice particle habits in Arctic clouds. Geophysical Research Letters. 1999; 26:1299-1302.
4. Korolev A.V., Isaac G.A., Hallett J. Ice particle habits in stratiform clouds. Quarterly Journal of the Royal Meteorological Society. 2000; (126):2873-2902.
5. Mishchenko M.I., Hovenier J.W., Travis L.D. Light scattering by nonspherical particles: Theory, measurements and geophysical applications. San Diego: Academic Press; 1999: 721.
6. Liu C., Panetta R.L., Yang P. The effective equivalence of geometric irregularity and surface roughness in determining particle single-scattering properties. Optics Express. 2014; 22(19):23620-23627.
7. Shishko V.A. Issledovanie opticheskikh svoystv atmosfernykh ledyanykh kristallov nepravil'noy formy [Investigation of the optical properties of irregularly shaped atmospheric ice crystals. Cand. Thesis]. Tomsk: IOA SO RAN; 2020: 143. (In Russ.)
8. Mu Q., Kargin B.A., Kablukova E.G. Computer-aided construction of three-dimensional convex bodies of arbitrary shapes. Computational Technologies. 2022; 27(2):">! (¡1. D01:10.25743/ICT.2022.27.2.005.
9. Preparata F.P., Shamos M.I. Computational geometry: An introduction. Springer Science Business Media; 1985: 413.
10. Santalo L. Integral geometry and geometric probability. USA: Addison-Wesley Publishing Company; 1976: 413.
11. De Berg M., van Krefeld M., Overmars M., Schwarzkopf O. Computational geometry algorithms and applications. 3rd rev. ed. Springer-Verlag; 2008: 386.
12. Zhuravlev V.F. Osnovy teoreticheskoy mekhaniki [Fundamentals of theoretical mechanics]. Moscow: Fizmatlit; 2008: 304. (In Russ.)
13. Mikhailov G.A., Voytishek A.V. Chislennoe statisticheskoe modelirovanie. Metody Monte-Karlo [Numerical statistical modelling. Monte Carlo Methods]. Moscow: IzdatePskiy tsentr "Akademiya"; 2006: 368. (In Russ.)
14. Marchuk G.I., Mikhailov G.A., Nazaraliev M.A., Dacbinjan R.A., Kargin B.A., Elepov B.S. Monte Carlo methods in atmospheric optics. Berlin: Springer-Verlag; 1980: 218.
15. Liu C., Panetta R., Yang P. The effects of surface roughness on the scattering properties of hexagonal columns with sizes from the Rayleigh to the geometric optics regimes. Journal of Quantitative Spectroscopy and Radiative Transfer. 2013; (129):169-185.
16. Pearson V., Newman G., James R. Vetrovye volny [Wind waves]. Moscow-Leningrad: Gidrometeoizdat; 1962: 42-124. (In Russ.)
17. Rakimgulov K.B., Ukhinov S.A. Local estimates in Monte Carlo method for the ocean atmosphere system with a random interface. Russian Journal of Numerical Analysis and Mathematical Modelling. 1994; 9(6):547-564.
18. Takano Y., Liou K.N. Solar radiative transfer in cirrus clouds. Part I: Single-scattering and optical properties of hexagonal ice crystals. Journal of Atmospheric Sciences. 1989; 46(1):3-19.
19. Konoshonkin A., Kustova N., Borovoy A. Rasseyanie sveta na geksagonal'nykh ledyanykh kristallakh peristykh oblakov [Light scattering by hexagonal ice crystals of cirrus clouds]. Saabryuken: LAP LAMBERT Academic Publishing; 2013: 156. (In Russ.)
20. Konoshonkin A.V. Rasseyanie sveta na atmosfernykh ledyanykh kristallakh pri lazernom zondirovanii [Light scattering by atmospheric ice crystals during laser sensing. Dr. thesis]. Tomsk: IOA SO RAN; 2017: 283. (In Russ.)
21. Konoshonkin A.V., Kustova N.V., Shishko V.A. Borovoi A.G. The technique for solving the problem of light backscattering by ice crystals of cirrus clouds by the physical optics method for a lidar with zenith scanning. Atmospheric and Oceanic Optics. 2016; 29(03):252-262.
22. Macke A., Johannes M., Ehrhard R. Single scattering properties of atmospheric ice crystals. Journal of the Atmospheric Sciences. 1996; 53(19):2813-2825.
