Имитационное моделирование формирования отклика системы "реакция - диффузия" на воздействие движущегося источника Текст научной статьи по специальности «Математика»

Математика.Прикладная математика.

Механика

УДК 519.63:519.688:004.94

М.В. Большаков, А.Г. Масловская

ИМИТАЦИОННОЕ МОДЕЛИРОВАНИЕ ФОРМИРОВАНИЯ ОТКЛИКА СИСТЕМЫ «РЕАКЦИЯ - ДИФФУЗИЯ» НА ВОЗДЕЙСТВИЕ ДВИЖУЩЕГОСЯ ИСТОЧНИКА

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

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

RESPONSE SIMULATION OF «REACTION - DIFFUSION» SYSTEM EXPOSED TO MOVING SOURCE

The paper presents the multidimensional mathematical model of «reaction - diffusion» process with perturbation under exposure of moving source. The computational scheme based on implicit splitting method was constructed to realize the model. The program implementation was also proposed to solve the problem. The computing experiment with dynamic response visualization of simulated system was performed.

Key words: «reaction - diffusion» process, simulation, parabolic type equation with perturbation, finite-difference scheme, stability, application program, computing experiment.

Введение

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

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

Можно отметить, что различные типы диффузионных систем описываются подобными законами. Математическое моделирование диффузионных процессов во многих случаях приводит к нелинейным и квазилинейным уравнениям с частными производными. Примерами приложений, использующих нелинейные эволюционные уравнения математической физики, являются: уравнение Шре-дингера, описывающее вероятность состояния и динамику частиц в квантовой механике; уравнение Нагумо, моделирующее распространение нервных импульсов в биологических системах; модели динамики цен в опционах; уравнение Гинзбурга - Ландау для формализации фазовых переходов, сверхпроводимости и процессов волоконной оптики и др. [2-3]. Аналитические решения нелинейных или квазилинейных уравнений с частными производными во многих случаях вызывают существенные затруднения, в связи с чем широкое распространение на практике получили приближенные мето-

Моделирование, приближенное к реальным условиям эксперимента, часто требует описания отклика системы на воздействие движущегося источника (тепла, зарядов, источника загрязнения среды и др.) [1; 3; 5]. В частности, в ряде авторских работ [6-8] таким источником является электронный зонд растрового электронного микроскопа, использующийся для анализа и модификации свойств полярных диэлектриков. При этом особенность моделирующих систем связана с аналитическим заданием функции источника и с реализацией алгоритма, определяющего перемещение источника и суперпозицию вкладов отклика среды на его воздействие. В связи с этим актуальным направлением является разработка специального математического и программного обеспечения для моделирования реакционно-диффузионных систем в условиях воздействия на объект движущегося источника.

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

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

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

ды [4].

Математическая постановка задачи моделирования

(1)

где и (х, у, ?) - определяемая функция; В - коэффициент диффузии; С > 0 - параметр при реакционном слагаемом; / ( х, у) - функция источника.

Рис. 1. Модельное представление геометрии взаимного положения объекта и источника, движущегося с постоянной скоростью.

Пусть источник начинает движение из позиции с координатами (0, L/2) и движется с постоянной скоростью V в положительном направлении оси OY. Для выполнения требования математической замкнутости модели дополним уравнение (1) начальным условием:

и ( х, у, tй ) = ы0. (2)

Определим границы Гь Г2, Г3, Г4 объекта G\ таким образом, что на границах Г3 и Г4 задаются условия I рода (устанавливается значение функции, соответствующее начальному и0):

' (3)

иГ =и0,

3

и |Г = и0,

4

а на границах Г и Г2 - условия II рода (нормальная компонента градиента функции и равна нулю и нулевым является поток через указанные границы): ди дх

= 0, ди

ду

= 0.

(4)

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

Построение конечно-разностной схемы решения задачи

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

®Ц = {х, = СТ = 1,N;у} = Л,] = 1,М;tk = кх, к = 1,К},

где Н\, С2 - шаги по пространственным координатам; т - шаг по времени.

Конечно-разностная схема конструируется в два этапа с введением аппроксимации для уравнения и граничных условий на двух временных слоях, первый из которых соответствует полушагу по времени т/2, второй - шагу т.

Первая подсхема получается расщеплением по координате х на дробном шаге т/2:

ик+Ш _ ик п

иц иц В I к+1/2 т к+1/2 . к+1/2 \ (г° к\ к+1/2 . г

—т—=С2( ' _ + '_( 'ир +

V

(5)

для всех 1 = 2,Ы _ 1, V = 1,М _ 1.

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

С

ик+у2 +

2 2, V

^ 2Вт к ^ 1++Сикт

иГГ = иI, + Ат, ищ = 0, V = .

(6)

Вторая подсхема получается аналогично расщеплением по координате у на полном дробном

шаге т:

ик+1 ик+1/2 и _ и _ В

(ик;} _ 2ик+1 + ик+ ), т = 1,М _ 1, V = 2,М _ 1,

1 +

2Вт

"АГ у

к+1 1л

2Вт

С

к+1 = . к+1/2 2 и 2 — и,д ,

и'1'--— и"2' = и~г ~, и,М = 0, Т = 1, N.

Схема имеет второй порядок точности по координате и первый - по времени. Итоговые системы линейных алгебраических уравнений имеют трехдиагональный вид и эффективно решаются методом прогонки на каждом временном слое. В уравнениях (5) и (6) для аппроксимации в правой части нелинейного слагаемого Си2 используется метод замороженных коэффициентов [12]. Расчет реакционного слагаемого Си2 происходит заменой квадрата функции на произведение значений дискретной функции на слоях ^ и ¿+т/2: ик • ик+1/2 . Такой прием позволяет избежать решения системы

и и

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

Устойчивость вычислительной схемы

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

и (X ^ 1) = ЕЕипш ехР(1 (V + кУ)).

п т

Необходимое условие устойчивости для конечно-разностного аналога исходной дифференци-

альной задачи имеет вид:

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

рую подсхемы, в итоге получим:

ик А Вт 1 „ . 2 кпС

—— = 1 + о • ¿1 + Ст, о =—, ¿1 = 4sln2^-}-ик+1/2 11 1 С2 2

к+1/2

и

—т—= 1 + сгД, = —-, Ь2 = 4sln

ик+1 2 2 2 С 2 2

Вт , . . 2 ктС2

Далее с помощью (11) - (12) составим требуемое соотношение

1 1

, к+1 , к+1 к+1/2

и и и

к к+1/2 к

и и и

1 + о2 • ¿2

1 + о1 • Ь1 + Ст

(11) (12)

(13)

поскольку всегда о2Ь2 > 0, о1Ь1 > 0, Ст > 0, то условие дробных шагов абсолютно устойчива.

< 1 выполняется, т.е. данная схема метода

2

С

т

2

и

к

и

и

к

и

Алгоритм решения задачи и режимы работы прикладной программы

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

Рис. 2. Блок-схема алгоритма реализации математической модели.

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

Рис. 3. Главная оконная форма интерфейса программы.

Программа требует инициализации следующих параметров: D - коэффициент диффузии; С -параметр при реакционном слагаемом; L - линейный размер объекта; Т - период времени наблюдения; а1, а2 - параметры функции источника. Параметры X и Y - геометрические параметры расчетной области (соответствуют L и определяются в основном программном модуле). Параметры С2, т -шаги по пространственным координатам и времени инициализируются в основном программном модуле и связаны с числом разбиений N расчетной области. Выходные данные программы: и(х,у,0 -графический вид и матричное представление искомой функции в динамике в соответствии с заданным временем наблюдения Т.

Вычислительный эксперимент и анализ результатов

Вычислительный эксперимент № 1. Цель его - проверка практической сходимости вычислительной схемы. Предположим, что по истечении длительного времени в системе установился стационарный режим. Аналитически оценить время, необходимое для установления стационарного

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

a1t

безразмерный критерий подобия Фурье [1] в теории теплопроводности Fo = ——, где а - коэффици-

х

ент тепловой диффузии материала, позволяет ввести масштаб времени по отношению к масштабу по

х2

координате. В нашем случае для характерного масштаба времени будем иметь: tsca¡e = хса1е , где -

характерный размер градиентной зоны.

Используя последнее, можно заключить, что по истечении достаточного большого промежут-

х2

ка времени t >> *са1е задача перейдет в стационарный режим. Тогда уравнение параболического типа

будет трансформировано в уравнение Пуассона, представимое в операторном виде следующим образом:

Ь[и ] + П = D (ди + V Си2 + П = 0.

^дх ду )

Разностный оператор, соответствующий дифференциальному L[u], будет представлен в фор-

ме:

Ли* = Ц-(ик . - 2ик + и*,.) + Ц-(и* - 2и к. + ик ,)- Сик ,

У •+У 1-11 ) С2 \ У+1 V V-1) У '

к-1

где С = и у - значение искомой функции, вычисляемое на предыдущем временном шаге.

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

II л к , х11 II к+1 к ||

щим нормам: Аи. + / - норма погрешности численного решения, Ш. -иЛ - норма отклонения

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

Оценка сходимости при равных шагах по пространству к\=к2 и разных по времени т показана в таблице.

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

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

и1(х,у) при С 1=^2=0.1, а затем для этого же момента времени рассчитано и2(х,у) на более мелкой сетке: С1=С2=0.05. После чего для крупной сетки установлена норма разности двух решений для указанного

момента времени: Цм1. — и2|| « 0.012 . Можно отметить, достаточный уровень точности, который достигается при реализации схемы.

Результат тестирования вычислительной схемы (стационарный аналог): оценка практической сходимости (¿*=250 усл. ед.)

С1=С2 X лМ + / ч у ик+1 — и* Т Т

0.1 0.1 0.01 3.72-10—4

0.1 0.05 0.01 3.68-10—4

0.05 0.1 0.0025 0.001

0.05 0.05 0.0025 8.97 -10—4

Вычислительный эксперимент № 2. Целью данного эксперимента является моделирование динамического отклика системы «реакция - диффузия» на воздействие статичного источника. Моделирование проведем в нормированных единицах. Инициализируем следующие значения параметров моделирования:

1) шаги по пространству - 0.1 усл. ед. и времени - 0.1 усл. ед.;

2) геометрические размеры расчетной области G\ соответствуют L=10 усл. ед.;

3) геометрические размеры области действия источника G2: X = 1 усл. ед., У = 1 усл. ед.;

4) коэффициент диффузии D=2 усл. ед.;

5) параметр при реакционном слагаемом С=1 усл. ед.;

6) функция источника задается в виде функции Гаусса:

/ ( х У ) = /0 ехР

Vх 2 +( у—5 )2 — а

с параметрами: а!=1 усл. ед., а2=2 усл. ед.,/0=1 усл. ед., 8=П2 усл. ед.; 7) период времени наблюдения Т=2 усл. ед.

Источник непрерывно действует в точке с координатами (0,5) в течение всего периода наблюдения Т. Результаты (фрагменты анимации) показаны на рис. 4 (а, б).

Рис. 4. Фрагменты анимации отклика системы на воздействие неподвижного источника

в фиксированные моменты времени.

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

а

2

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

Вычислительный эксперимент № 3. Целью данного эксперимента является моделирование динамического отклика системы «реакция-диффузия» на воздействие движущегося источника. Моделирование проведем в нормированных единицах. Инициализируем аналогичные предыдущему случаю значения параметров моделирования. Кроме этого допустим, что источник перемещается из позиции (0,5) с постоянной скоростью 1 усл. ед.

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

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

Заключение

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

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

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

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

1. Мартинсон, Л.К., Малов, Ю.И. Дифференциальные уравнения математической физики - М.: МГТУ им. Н.Э. Баумана, 2011. - 368 с.

2. Otten, D. Mathematical Models of Reaction Diffusion Systems, their Numerical Solutions and the Freezing Method with Comsol Multiphysics. - Bielefeld: Bielefeld University, 2010. - 77 p.

3. Kuttler, C. Reaction-Diffusion equations with applications. - Munich: Technical University Press, 2011. - 101 p.

4. Harwood, R.C. Operator splitting method and applications for semilinear parabolic partial differential equations.

- Washington: Washington State University, 2011. - 112 p.

5. Drake, J.B. Climate modeling for scientists and engineers. - Knoxville, Tennessee: University of Tennessee, 2014. - 46 p.

6. Масловская, А.Г. Анализ тепловых эффектов, возникающих при взаимодействии электронных пучков с сегнетоэлектрическими кристаллами // Известия вузов. Физика. - 2010. - № 1. - С. 34-40.

7. Maslovskaya, A., Pavelchuk, A. Simulation of heat conductivity and charging processes in polar dielectrics induced by electron beam exposure // IOP Conf. Series: Materials Science and Engineering. - 2015. - V. 81. - P. 012119 (6).

8. Pavelchuk, A.V., Maslovskaya, A.G. Simulation of internal charge distribution and spatial charge characteristics of ferroelectrics irradiated by focused electron beam // Proc. SPIE 10176, Asia-Pacific Conference on Fundamental Problems of Opto- and Microelectronics (December 14, 2016). - P. 101760P (12).

9. Яненко, H.H. Метод дробных шагов решения многомерных задач математической физики - Новосибирск: Наука, 1967. - 197 с.

10. Tan, X. A splitting method for fully nonlinear degenerate parabolic PDEs // Electron J. Prorab, 2013. - No 145.

- P. 1-24.

11. Формалев, В.Ф., Ревизников, Д.Л. Численные методы - М.: ФИЗМАТЛИТ, 2004. - 400 с.

12. Годунов, С.К., Рябенький, B.C. Разностные схемы (введение в теорию). Учебное пособие - М.: Наука, 1977. - 440 с.