Approach to numerical implementation of drift-diffusion model of field effects induced by moving source.pdf Введение В настоящее время в числе объектов математического моделирования важное место занимают системы, в которых протекают процессы диффузионной природы [1, 2]. Уравнения типа «конвекция - реакция - диффузия» (или «адвекция - реакция - диффузия») используют для описания широкого класса явлений в самых различных предметных областях: в задачах экологического прогнозирования - при моделировании распространения загрязняющих веществ [3]; в медицине - для описания гемодинамики в крупных кровеносных сосудах [4]; в экономике - для моделирования поведения цен в европейских и азиатских опционах (уравнение Блэка - Шоулза) [5]; в физике - для квантово-механических расчетов функций плотности вероятности частиц (в постановке уравнений Фоккера - Планка и Колмогорова) [6] и для моделирования тепломассопереноса и транспортных механизмов (диффузии и дрейфа) заряженных частиц, например электронов в диэлектрических материалах [7, 8]. Математическая формализация моделей конвективно-реакционно-диффузионных процессов приводит к начально-граничным задачам для уравнений с частными производными параболического типа при реализации нестационарных режимов в моделируемых системах или к граничным задачам для уравнений эллиптического типа в стационарных режимах [7, 8]. Сложность математических постановок (формализации дифференциальной задачи и задания граничных условий) обуславливает востребованность развития численных методов решения данного класса задач. Одними из широко распространенных методов решения многомерных эволюционных задач математической физики являются конечно-разностные методы расщепления [7]. Методы этой группы при решении задач с несложной геометрией расчетных областей не уступают по точности методу конечных элементов, при использовании в последнем случае линейных базисных функций. При этом конечно-разностные методы легче программируются и дают выигрыш в плане времени машинных вычислений при решении ресурсоемких задач, в частности при решении задач с запаздыванием [9]. Специфика конечно-разностного решения конвективно-диффузионных задач заключается в построении монотонных вычислительных схем с требуемым порядком аппроксимации дифференциального уравнения и граничных условий. Присутствие дрейфового (конвективного) слагаемого накладывает определенные ограничения на устойчивость и монотонность рассматриваемых задач. Для получения численного решения конвективно-реакционно-диффузионного уравнения с порядком точности выше первого необходимо либо использовать ограниченную по монотонности схему с аппроксимацией производной центральными разностями, либо применять специальные схемы повышенной точности: кососимметрические разностные схемы, схемы высокого порядка с аппроксимацией конвективного слагаемого по большему числу узлов, схемы с искусственными монотонизаторами или схемы с использованием адаптивной искусственной вязкости и др. [10-13]. К рассматриваемым классам прикладных задач можно отнести и диффузионно-дрейфовую модель процесса зарядки полярных диэлектриков в условиях электронного облучения [14]. Развитию диффузионно-дрейфового подхода и численному моделированию процесса зарядки диэлектриков посвящен ряд современных работ [15-17]. В серии авторских публикаций (например, [18, 19]) описаны результаты разработки средств и методов физико-математического моделирования полевых эффектов инжектированных зарядов в сегнетоэлектриках. Модель описывается системой соотношений, включающей уравнение «конвекция - реакция - диффузия» для определения пространственного распределения объемной плотности зарядов, локально-мгновенное уравнение Пуассона для расчета распределения потенциала и уравнение, выражающее связь между напряженностью поля инжектированных зарядов и потенциалом. При моделировании процесса зарядки диэлектрика пучком электронов в условиях, приближенных к экспериментальным режимам зондирования, часто требуется описать отклик системы при воздействии на объект движущимся источником. Особенность моделирования динамики таких процессов связана с аналитическим заданием функции источника и с вычислительной реализацией алгоритма, определяющего перемещение источника и суперпозицию вкладов отклика среды на его воздействие. В связи с этим настоящая работа направлена на развитие математического, алгоритмического и программного обеспечения для реализации диффузионно-дрейфовой модели полевых эффектов, наблюдаемых в объектах в режиме динамического воздействия сосредоточенного источника. 1. Математическая постановка задачи моделирования диффузионно-дрейфовой системы Сформулируем концептуальную постановку задачи. Будем считать, что в некотором объекте, заданном двумерной прямоугольной областью G1 с характерными линейными размерами X и Y, действует внутренний сосредоточенный источник, зависящий от пространственных координат и занимающий в объекте геометрическую область G2. Источник движется с постоянной скоростью v в одном из координатных направлений. Требуется провести математическое моделирование отклика диффузионно-реакционно-дрейфовой системы на воздействие движущегося источника. Математическую постановку задачи введем в рассмотрение в виде одной из частных постановок, которая в общем виде соответствует математической модели оценки полевых эффектов инжектированных электронов в диэлектриках. Сформулируем начально-граничную задачу для двумерного уравнения типа «конвекция - реакция - диффузия»: , (1) где u(x, y, t) - искомая функция; D - коэффициент диффузии; vc > 0 - скорость конвективного переноса; с - реакционный коэффициент; f(x, y) - функция внутреннего источника; - конечный момент времени наблюдения процесса. Будем полагать, что внутренний источник задается геометрически половиной окружности с радиусом R. Источник начинает действовать в позиции x = 0, y = S в момент времени t0 и перемещается вдоль оси ОY с постоянной скоростью v. Схематическое представление взаимного расположения объекта и действующего источника приведено на рис. 1. Рис. 1. Геометрия расчетной области Для замыкания математической постановки задачи (1) введем в рассмотрение начальное условие (2) и граничные условия (3) Граничные условия II рода, заданные на границе y = 0, обеспечивают непрерывность решения задачи (позволяют, в силу симметрии задачи, строить решение для половины расчетной области), а на границе x = 0 - соответствуют физическому смыслу задачи (отсутствие потока через верхнюю грань объекта). Граничные условия I рода отвечают условию затухания переходных процессов и стабилизации полевых характеристик на границах объекта, значительно удаленных от градиентной зоны. При вариации значений характерных параметров или условий реализации модели приоритетное влияние на поведение диффузионно-дрейфовой системы может оказывать либо диффузия, либо конвекция [8]. Для оценки соотношения вкладов диффузии и конвекции используют безразмерный критерий Пекле: , (4) где vc - абсолютная величина скорости переноса, м/с; xsc - масштаб расстояния (линейный размер градиентной зоны), м; D - коэффициент диффузии, м2/с. Анализ (4) позволяет заключить, что при Pe > 1 процесс характеризуется доминированием процесса переноса и слабой диффузией. В первом случае задачи относят к виду регулярно возмущенных задач, во втором - к виду сингулярно возмущенных. Сингулярно возмущенные задачи могут отличаться наличием областей сильного изменения решения. Для определения характерного времени наблюдения динамики процесса можно использовать критерий подобия Фурье: , (5) где tsc - характерное время наблюдения динамики процесса. Поэтому при временах будет наблюдаться изменение характеристик системы с течением времени, при система будет соответствовать режиму, близкому к стационарному состоянию. 2. Вычислительная схема для реализации модели диффузионно-реакционно-дрейфового процесса Для построения конечно-разностной схемы решения задачи (1) - (3) был использован гибридный подход, состоящий в модификации неявно-явной конечно-разностной схемы расщепления для решения многомерной эволюционной задачи математической физики. Одной из вычислительных проблем, связанных с конечно-разностным решением уравнений вида (1), является аппроксимация дрейфового слагаемого. Использование центрально-симметричной аппроксимации приводит к потере монотонности решения (особенно для задач с «большими» значениями числа Пекле), а применение противопотоковой («upsream») схемы хотя и дает монотонное решение, но, вместе с тем, понижает порядок аппроксимации по координатам всей схемы до первого [8]. С целью аппроксимации дрейфовой компоненты была применена вычислительная схема Робертса - Вейсса. Ранее, в работе [13], эта схема была использована в части аппроксимации конвективного слагаемого для решения одномерного уравнения «конвекция - реакция - диффузия» при реализации модели опционов Блэка - Шоулза. Введем конечно-разностную сетку с шагами h1 и h2 по координатам и  - по времени: . Схема является трехслойной. Начальное условие аппроксимируется в виде . Первые две подсхемы («предиктор») являются неявными и получаются расщеплением разностного оператора по координатам х и y на дробном шаге /2, который делится на полушаг /4, в совокупности со схемой Робертса - Вейсса для аппроксимации конвективного слагаемого: (6) (7) Конечно-разностные уравнения (6) и (7) определены для всех внутренних узлов сетки . Подсхемы (6) и (7) дополняются аппроксимацией граничных условий I рода: , и , соответственно. Краевые условия II рода учитываются в модифицированных разностных уравнениях (6), (7) с использованием аппроксимации второго порядка точности и дополнительных слоев по координатам с фиктивными узлами: , Итоговые системы линейных алгебраических уравнений решаются методом прогонки на каждом временном подслое. Третья подсхема является явной («корректор») и на полном дробном шаге  позволяет вычислить значение искомой функции в сочетании с аппроксимацией конвективного слагаемого по схеме центральной разности, с учетом функции источника и аппроксимации реакционного слагаемого (для всех ): (8) В подсхеме (8) также учитываются аппроксимации граничных условий I и II рода. Можно показать, что итоговая схема является абсолютно устойчивой, имеет второй порядок аппроксимации по координатам и времени , а также монотонна. В разностных уравнениях (6) - (8) роль предиктора состоит в обеспечении устойчивости всей схемы, а функция корректора заключается в повышении порядка аппроксимации по времени дифференциальной задачи. Применение схемы Робертса - Вейсса обеспечивает монотонность реализуемых неявных подсхем. 3. Компьютерная реализация диффузионно-дрейфовой модели системы в режиме динамического воздействия источника 3.1. Инициализация и контроль параметров вычислительного эксперимента Для решения рассмотренного класса задач вычислительный алгоритм, соответствующий схеме (6) - (8), был реализован в виде программного приложения в ППП Matlab. Программа позволяет наблюдать динамику изменения искомой функции u(x,y) в заданный период времени. Проведение вычислительного эксперимента предусматривает инициализацию параметров математической модели: скорости конвективного переноса vс, м/с; коэффициента диффузии D, м2/с; линейного размера объекта X = Y = L и источника R, м; периода времени наблюдения , с (в соответствии с критерием Фурье). Полагая, что искомая полевая характеристика процесса имеет некоторую физическую размерность [u] (например, объемная плотность зарядов - Кл/м3), необходимо задать константу реакционной компоненты с, 1/([u]с), и генерационное слагаемое - функцию источника f(x,y), [u]/с. Для управления вычислительным процессом требуется задать шаги по координатам h1, h2 и времени . Выбор значений шагов по пространственным координатам и времени должен обеспечивать выполнение естественных (физических) условий устойчивости, отнесенных к диффузионному, дрейфовому и реакционному членам соответственно: , , , . (9) В прикладной программе предусмотрено, что источник перемещается из позиции (0, S) (например, S = L/2) с постоянной скоростью, равной v, которая задается через расстояние, пройденное за один временной шаг . Функция внутреннего источника инициализируется с помощью распределения Гаусса с параметрами f0, a1 и a2 (последние имеют размерность расстояния): . (10) 3.2. Тестирование вычислительной схемы Для верификации работы модулей прикладной программы проведено тестирование практической сходимости схемы, реализующей нестационарный режим. С этой целью использована схема двойного пересчета в двух вариантах. Моделирование проводилось для сетки с шагами h p = h и h p+1 = h/2 (h = h1 = h2) для одного и того же временного слоя с одинаковым шагом по времени. В тесте принято значение шага по времени  = 0.001 отн. ед., период времени наблюдения - 0.1 отн. ед., параметры D и vc нормированы к единичным значениям, реакционное слагаемое отсутствует, источник неподвижен, L = 10 отн. ед. Результат оценки погрешности с использованием относительной невязки для расчета с шагом h p: приведен на рис. 2, а. Можно сделать вывод, что допустимым представляется выбор числа разбиений N = M = 200-400 (соответственно при шагах сетки h1 = h2 = 0.05-0.025 отн. ед.). Также проводились вычисления для сетки с фиксированным значением шагов по координате (h = 0.05 отн. ед.) при варьировании шага по времени, период времени наблюдения составил 0.1 отн. ед. Результат приведен на рис. 2, б. Рис. 2. Оценка относительной погрешности решения при варьировании: шагов по координатам h = h1 = h2 (а), шага по времени dt =  (б) Варьирование шага по времени dt =  от 0.001 до 0.005 отн. ед. не дает существенного увеличения точности вычислений по порядку измеряемых величин. Подобный тест для каждого набора модельных параметров позволяет контролировать точность при определенных значениях шагов по координатам и времени. Основополагающим аргументом при выборе параметров вычислительного процесса является удовлетворение условиям (9). 3.3. Результаты вычислительных экспериментов Результат работы программы продемонстрируем на примере моделирования диффузионно-дрейфовой системы при следующих значениях параметров: D = 1, vc = 1, L = 10, с = 1, а1 = 1, а2 = 1, f0 = 1, R = 0.5. Источник перемещается вдоль координатной оси oy с постоянной скоростью v = 0.01 отн. ед. Фрагменты представления анимированного решения показаны на рис. 3. Рис. 3. Фрагменты анимации - модельное представление искомой функции u(x, y) в фиксированные моменты времени: t1 = 0.1 отн. ед. (а), t2 = 0.2 отн. ед. (б) Основной особенностью визуализированного численного решения задачи (1) - (3) является асимметричное смещение изолиний графика многомерной функции и изменение амплитудного значения в точке воздействия источника в течение всего периода воздействия. Спустя время, большее отн. ед. (по критерию Фурье (5)), состояние системы будет близко к стационарному. В данном случае рассмотрен пример с отсутствием доминирования в соотношении процессов диффузии и дрейфа (число Пекле Pe  1). В условиях конкретной прикладной задачи этот аспект требует детального рассмотрения, поскольку доминирование дрейфа (конвекции) может привести к сингулярному возмущению самой задачи. На рис. 4, а представлен результат моделирования в тех же условиях в случае доминирования диффузии (D = 10 отн. ед.). Как видно, меняется геометрическая форма распределения полевой характеристики (диффузия приводит к «размытию» профиля функции u(x,y) и уменьшению максимального значения в точке воздействия источника). Кроме того, в приложениях дополнительного контроля требуется и анализ значения реакционного коэффициента с. На рис. 4, б показан результат эксперимента при усилении реакционной компоненты (с = 1000 отн. ед.). В последнем случае также наблюдается изменение значений и профиля искомой функции. Рис. 4. Решение u(x, y) в фиксированные момент времени t1 = 0.1 отн. ед. при варьировании параметров: коэффициент диффузии D = 10 (а), константа реакции с = 1000 (б) Заключение Таким образом, работа обобщает результат построения и программной реализации модели диффузионно-реакционно-дрейфовой системы в режиме динамического воздействия. Вычислительный алгоритм сконструирован на основе гибридного подхода, сочетающего модификацию схемы расщепления типа «предиктор - корректор» для аппроксимации диффузионной части, а также схему Робертса - Вейсса и центрально-симметричную формулу для аппроксимации дрейфовой компоненты. Вычислительная схема является абсолютно устойчивой, имеет второй порядок аппроксимации по координатам и времени, монотонна. Разработана прикладная программа в ППП Matlab, позволяющая визуализировать анимацию изменения искомой функции с течением времени. Представлены результаты вычислительных экспериментов для решения модельной задачи в нормированном виде при варьировании значений управляющих параметров модели в условиях доминирования диффузии и усиления реакционной компоненты.

Скачать электронную версию публикации