Журнал вычислительной математики и математической физики, 2022, T. 62, № 10, стр. 1696-1706

1. ВВЕДЕНИЕ. ПОСТАНОВКА КРАЕВОЙ ЗАДАЧИ

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

Детерминированный подход приводит к фундаментальному описанию моделей конвективно-реакционно-диффузионных процессов в постановке начально-граничных задач для уравнений с частными производными параболического типа в неустановившихся режимах и к краевым задачам для уравнений эллиптического типа в случае описания стационарных состояний. К указанному классу математических моделей можно отнести диффузионно-дрейфовую модель процесса зарядки диэлектриков в неравновесных внешних условиях (см. [1]–[8]). В терминологии данной предметной области вводится в рассмотрение “дрейф” носителей заряда (электронов и дырок), представляющий некоторый аналог процессов конвекции или адвекции, рассматриваемых в других приложениях. Классическая диффузионно-дрейфовая модель зарядки диэлектриков представлена во многих фундаментальных обзорах по теоретической физике полупроводников и диэлектриков (см., например, [9]).

Одной из важнейших частных задач является развитие диффузионно-дрейфового подхода для моделирования процесса зарядки полярных диэлектриков, индуцированного электронным облучением. Практический интерес данная область вызывает в связи с необходимостью прогнозировать состояние функциональных диэлектрических материалов при диагностике и модификации их свойств методами растровой электронной микроскопии. Разработке фундаментальных основ, развитию математических моделей, созданию математического и программного обеспечения для исследования процессов электронно-стимулированной зарядки посвящен широкий ряд современных работ [1]–[7] и [10]–[16].

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

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

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

Математическая модель процесса зарядки полярных диэлектриков при достаточно длительном электронном облучении может быть представлена следующей краевой задачей, рассматриваемой в ограниченной области $\Omega \subset {{\mathbb{R}}^{3}}$ с границей $\Gamma $:

Здесь $\rho $ – объемная плотность заряда, ${\mathbf{E}}$ – вектор-функция напряженности электрического поля, $d$ – коэффициент диффузии электронов, ${{\mu }_{n}}$ – дрейфовая подвижность электронов, $\varepsilon $ – диэлектрическая проницаемость материала, ${{\varepsilon }_{0}}$ – электрическая постоянная, $f$ – генерационное слагаемое, отвечающее за действие объемного источника зарядов в объекте. Ниже на задачу (1.1)–(1.3) будем ссылаться как на задачу 1.

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

Отметим работы [17]–[23], посвященные исследованию конвективно-реакционно-диффузионных моделей, а также статьи [24]–[27] по исследованию близких моделей сложного теплообмена.

2. РАЗРЕШИМОСТЬ КРАЕВОЙ ЗАДАЧИ

При анализе краевой задачи будем использовать функциональные пространства Соболева ${{H}^{s}}(D)$, $s \in \mathbb{R}$. Здесь $D$ обозначает область $\Omega $, либо некоторую подобласть $Q \subset \Omega $, либо границу $\Gamma $. Через ${{\left\| {\, \cdot \,} \right\|}_{{s,Q}}},\;{{\left| {\, \cdot \,} \right|}_{{s,Q}}}$ и ${{( \cdot , \cdot )}_{{s,Q}}}$ будем обозначать норму, полунорму и скалярное произведение в ${{H}^{s}}(Q)$. Нормы и скалярные произведения в ${{L}^{2}}(Q)$ и ${{L}^{2}}(\Omega )$ будем обозначать соответственно через ${{\left\| {\, \cdot \,} \right\|}_{Q}}$ и ${{( \cdot , \cdot )}_{Q}}$, ${{\left\| {\, \cdot \,} \right\|}_{\Omega }}$ и ${{( \cdot , \cdot )}_{\Omega }}$.

Введем функциональные пространства

где ${\text{ker}}{\kern 1pt} ({\text{rot}})$ – ядро оператора ротор, и произведение пространств $X = H_{0}^{1}(\Omega ) \times \tilde {H}_{N}^{1}(\Omega )$.

Предположим, что выполняются следующие условия:

(i) $\Omega $ – ограниченная область в ${{\mathbb{R}}^{3}}$ с границей $\Gamma \in {{C}^{{0,1}}}$;

(ii) $f \in {{L}^{2}}(\Omega )$.

Напомним, что в силу теоремы вложения Соболева пространство ${{H}^{1}}(\Omega )$ вкладывается в пространство ${{L}^{s}}(\Omega )$ непрерывно при $s \leqslant 6$, компактно при $s < 6$ и с некоторой константой ${{C}_{s}}$, зависящей от $s$ и $\Omega $, справедлива оценка

Ниже будем использовать следующие формулы Грина (см. [28], [29]):

Справедлива

Лемма 2.1 (см. [28]). При выполнении условий (i), ${\mathbf{E}} \in {{H}^{1}}{{(\Omega )}^{3}}$ существуют положительные константы ${{C}_{0}},\;{{\delta }_{1}},\;\gamma _{1}^{'}$ и ${{\gamma }_{1}}$, зависящие, соответственно, от $\Omega $, такие, что

Если функции ${\mathbf{E}} \in {{H}^{1}}{{(\Omega )}^{3}}$ и $\rho \in H_{0}^{1}(\Omega )$ связаны вторым соотношением в (1.2), то справедливо равенствопринимающее при $h = \rho $ следующий вид:

Доказательство. Оценки (2.4) вытекают из неравенства Гёльдера и (2.1). Неравенство (2.5) представляет собой неравенство Фридрикса-Пуанкаре.

Для доказательства (2.6) сначала заметим, что $h{\mathbf{E}} \in {{L}^{3}}{{(\Omega )}^{3}}$ и справедливо равенство

Из (2.8) вытекает, что $\operatorname{div} (h{\mathbf{E}}) \in {{L}^{{3/2}}}(\Omega )$. Это позволяет применить к $({\mathbf{E}} \cdot \nabla \rho ,h)$ формулу Грина (2.3). Применив, приходим к равенству С учетом второго уравнения в (1.2) и (2.8) формула (2.9) принимает вид (2.6), частным случаем которой является формула (2.7).

Справедлива

Лемма 2.2 (см. [28]). При выполнении условия (i) для любой функции $\sigma \in {{L}^{2}}(\Omega )$ существует единственное решение ${\mathbf{E}} \in \tilde {H}_{N}^{1}(\Omega )$ задачи

для которого справедлива оценка где ${{C}_{N}}$ – положительная константа, зависящая от $\Omega $.

Пусть $(\rho ,{\mathbf{E}}) \in ({{C}^{2}}(\Omega ) \cap {{C}^{0}}(\bar {\Omega })) \times ({{C}^{1}}{{(\Omega )}^{3}} \cap \tilde {H}_{N}^{1}(\Omega ))$ – классическое решение задачи 1. Умножим уравнение в (1.1) на функцию $h \in H_{0}^{1}(\Omega )$ и проинтегрируем по $\Omega $, применяя формулу Грина (2.2). Приходим к слабой формулировке задачи 1

Разрешимость задачи (2.10), (2.11) докажем с помощью теоремы Шаудера. Для этого построим отображение $G:H_{0}^{1}(\Omega ) \to H_{0}^{1}(\Omega )$, действующее по формуле

Здесь функция $\rho \in H_{0}^{1}(\Omega )$ является решением линейной задачи

Из леммы 2.2 вытекает, что для любой функции $r \in H_{0}^{1}(\Omega )$ существует единственное решение ${\mathbf{E}} \in H_{N}^{1}(\Omega )$ задачи (2.13), для которого справедлива оценка

Из леммы 2.1 вытекает, что форма $a(\rho ,h)$ непрерывна. Так же из леммы 2.1 и равенства (2.13) следует, что

При $h = \rho $ равенство (2.15) принимает вид

Положим $h = \rho $ в (2.12). С учетом (2.16) приходим к равенству

Поскольку то из (2.17) и леммы 2.1 вытекает, что форма $a$ коэрцитивна на пространстве $H_{0}^{1}(\Omega )$ с константой $d{{\delta }_{1}}$. Тогда из теоремы Лакса-Мильграма следует, что при любом $r \in H_{0}^{1}(\Omega )$ существует единственное решение $\rho \in H_{0}^{1}(\Omega )$ задачи (2.12), для которого справедлива оценка

Из вышесказанного вытекает, что отображение $G:H_{0}^{1}(\Omega ) \to H_{0}^{1}(\Omega )$ определено корректно и переводит шар ${{B}_{R}}$ радиуса $R = {{C}_{*}}{{\left\| f \right\|}_{\Omega }}$ в себя.

Докажем, что $G$ непрерывен и компактен на ${{B}_{R}}$. Пусть $\{ {{r}_{k}}\} $, $k = 1,2...,$ – произвольная последовательность из ${{B}_{R}}$. В силу рефлексивности пространства ${{H}^{1}}(\Omega )$ и компактности вложения ${{H}^{1}}(\Omega ) \subset {{L}^{p}}(\Omega )$, $p < 6$, существует подпоследовательность последовательности $\{ {{r}_{k}}\} $, $k = 1,2,...$, которую мы опять обозначим через $\{ {{r}_{k}}\} $, и функция $r \in {{B}_{R}}$ такие, что

Пусть ${{\rho }_{k}} \in H_{0}^{1}(\Omega )$ – решение задачи

получающейся из (2.12), (2.13) заменой $r$ на ${{r}_{k}}$, $\rho $ на ${{\rho }_{k}}$ и ${\mathbf{E}}$ на ${{{\mathbf{E}}}_{k}}$.

Покажем, что

Поскольку последовательность $\{ {{r}_{k}}\} $ содержится в шаре ${{B}_{R}}$, то из (2.21) и леммы 2.2 вытекает, что

Вычитая (2.12) из (2.20), будем иметь

Полагая $h = {{\rho }_{k}} - \rho $ и учитывая (2.21), получаем, что С учетом (2.23) приходим к равенству В силу леммы 2.1 и неравенства Гёльдера из (2.24) выводим, что

Из (2.25) в силу (2.19) вытекает, что ${{\rho }_{k}} \to \rho $ сильно в ${{H}^{1}}(\Omega )$.

В таком случае оператор $G$ является компактным и непрерывным. Тогда из теоремы Шаудера вытекает, что оператор $G$ имеет неподвижную точку $\rho = G(\rho )$, которая является решением уравнения (2.10). В свою очередь, пара $(\rho ,{\mathbf{E}})$ является искомым решением задачи (2.10), (2.11). При этом для функции $\rho $ справедлива оценка (2.18). Тогда для электрического поля ${\mathbf{E}}$ из (2.13) в силу леммы 2.2 и с учетом (2.18) получаем оценку

Установим достаточные условия единственности решения задачи (2.10), (2.11). Обозначим через $({{\rho }_{1}},{{{\mathbf{E}}}_{1}}) \in X$ и $({{\rho }_{2}},{{{\mathbf{E}}}_{2}}) \in X$ любые два ее решения.

Несложно проверить, что разности

удовлетворяют соотношениям

Поскольку

то с учетом леммы 2.1 получаем, что

Равенство (2.27) при $h = \rho $ с учетом (2.29) принимает вид

Для левой части (2.30) в силу (2.5) справедлива оценка

Оценки для правой части (2.30) вытекают из неравенства Гёльдера, соотношения (2.4), и с учетом (2.1), (2.18) и (2.14) принимают следующий вид:

Окончательно приходим к неравенству

Пусть исходные данные задачи 1 таковы, что выполняется условие

Тогда из (2.31) вытекает, что ${{\left\| \rho \right\|}_{{1,\Omega }}} = 0$ или ${{\rho }_{1}} = {{\rho }_{2}}$ в $\Omega $. В таком случае из (2.28) получаем, что

Последнее в силу леммы 2.2 означает, что ${\mathbf{E}} = {\mathbf{0}}$ или ${{{\mathbf{E}}}_{1}} = {{{\mathbf{E}}}_{2}}$ в $\Omega $.

Сформулируем полученные результаты в виде следующей теоремы.

Теорема 2.1. При выполнении условий (i), (ii) существует слабое решение $(\rho ,{\mathbf{E}}) \in X$ задачи 1 и справедливы оценки (2.18) и (2.26). Если к тому же выполняется условие (2.32), то слабое решение задачи 1 единственно.

3. ПРИНЦИП МИНИМУМА И МАКСИМУМА

Положим ${{f}_{{{\text{max}}}}}$ – положительная константа.

Пусть в дополнение к (i)–(ii) выполнено условие

(iii) $0 \leqslant f \leqslant {{f}_{{{\text{max}}}}}$ п.в. в $\Omega $.

Рассуждая как в [30], установим принцип минимума и максимума для плотности заряда $\rho $.

Лемма 3.1. При выполнении условий (i)–(iii) для слабого решения $\rho \in H_{0}^{1}(\Omega )$ задачи 1 справедлив следующий принцип минимума и максимума:

Доказательство. Сначала докажем, что $\rho \geqslant 0$ п.в. в $\Omega $. С этой целью введем функцию $w = \min \{ \rho ,0\} $. Ясно, что оценка $\rho \geqslant 0$ выполняется тогда и только тогда, когда $w = 0$ п.в. в $\Omega $.

Через $Q \subset \Omega $ обозначим измеримое подмножество $\Omega $, в котором $w < 0$. Ясно, что $w \in {{H}^{1}}(\Omega )$, а из [31], [32] вытекает, что

Тогда $w \in H_{0}^{1}(\Omega )$.

Справедливо следующее равенство (см. [31]):

Из вышесказанного вытекает, что

С учетом (2.7) из второго равенства в (3.3) получаем, что

При этом

Полагая в (2.10) $h = w$, будем иметь

С учетом (3.2)–(3.5) соотношение (3.6) принимает вид Поскольку $w < 0$ в $Q$ и $f \geqslant 0$ п.в. в $\Omega $, то ${{(f,w)}_{Q}} \leqslant 0$. Тогда в силу неравенства Фридрихса–Пуанкаре (2.5) из (3.7) приходим к оценке из которой вытекает, что ${{\left\| w \right\|}_{{1,\Omega }}} = 0$, а следовательно, $w = 0$ п.в. в $\Omega $. Тогда $\rho \geqslant 0$ п.в. в $\Omega $.

Докажем принцип максимума. Для этого введем функцию $\tilde {\rho } = \max \{ \rho - M,0\} $, которая, как и вспомогательная функция $w$, принадлежит ${{H}^{1}}(\Omega )$. Ясно, оценка $\rho \leqslant M$ п.в. в $\Omega $ выполняется тогда и только тогда, когда $\tilde {\rho } = 0$ п.в. в $\Omega $.

Из [31], [32] вытекает, что

В таком случае, $\tilde {\rho } \in H_{0}^{1}(\Omega )$.

Через ${{Q}_{M}} \subset \Omega $ обозначим открытое измеримое подмножество $\Omega $, в котором $\rho > M$. Как и выше, из [31] вытекает, что

Из сказанного выше вытекает, что

С учетом леммы 1.2 приходим к равенству

Важную роль играет следующее равенство:

Подставляя $h = \tilde {\rho }$ в (2.10), с учетом (3.10), (3.11) получаем Вычитая $({{\mu }_{n}}{\text{/}}\varepsilon {{\varepsilon }_{0}})({{M}^{2}},\tilde {\rho }{{)}_{{{{Q}_{M}}}}}$ из обеих частей (3.12), будем иметь Применив к левой части (3.13) неравенство (2.5), приходим к оценке из которой выводим, что $\tilde {\rho } = 0$ п.в. в $\Omega $, если параметр $M$ определен в (3.1).

4. ВЫЧИСЛИТЕЛЬНЫЙ ЭКСПЕРИМЕНТ

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

В качестве примера представим результаты компьютерной реализации математической модели процесса стационарной зарядки сегнетоэлектрического кристалла ниобата лития ${\text{LiNb}}{{{\text{O}}}_{{\text{3}}}}$ в условиях облучения электронным пучком средних энергий (1–40 кэВ). Рассматривая электронный зонд как сфокусированный внутренний источник зарядов, можно провести расчет характеристик процесса электронно-индуцированной зарядки сегнетоэлектрика в рамках исследуемой математической модели [1], [4], [5], [8], [16].

Задачу 1 сформулируем в следующем виде:

где $0 < x < L,$ $0 < y < L,$ $0 < z < L$; $\rho (x,y,z)$ – объемная плотность зарядов, Кл/м³; параметр $L$ определяет геометрические размеры объекта, м; $\varepsilon $ – диэлектрическая проницаемость материала; ${{\varepsilon }_{0}} = 8.85 \times {{10}^{{ - 12}}}$ – электрическая постоянная, Ф/м; $f$ – генерационное слагаемое, Кл/(м³ с); $\varphi $ – распределение потенциала, В; ${\mathbf{E}}$ – напряженность поля, В/м; $D$ – коэффициент диффузии электронов, м²/с; ${{\mu }_{n}}$ – дрейфовая подвижность электронов, м²/(В ⋅ с).

В данном случае будем считать, что сфокусированный внутренний источник зарядов действует перпендикулярно плоскости поверхности верхней грани кристалла $z = 0$. Время экспозиции в экспериментах составляет от 10^–3 с до 10 с, что, с учетом оценки характерного времени переходного процесса (от 1 нс до 1 мкс), позволяет принять стационарный режим соответствующим физическому смыслу задачи.

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

Проведем инициализацию параметров вычислительного эксперимента. Модельный образец представляет собой куб с линейным размером $L = {{10}^{{ - 4}}}$ м. Позиция источника определяется координатами $(L{\text{/}}2,L{\text{/}}2,0)$. Установлены следующие параметры моделирования для кристалла ниобата лития: коэффициент диффузии электронов $D = 1.915 \times {{10}^{{ - 6}}}$ м²/с, подвижность электронов в материале ${{\mu }_{n}} = 74 \times {{10}^{{ - 6}}}$ м²/(В с), диэлектрическая проницаемость материала $\varepsilon = 30$. Генерационное слагаемое:

где ${{f}_{0}}$ – нормировочный коэффициент; $\eta $ – коэффициент вторичной электронной эмиссии; ${{d}_{1}}$, ${{d}_{2}}$ – параметры аппроксимации, устанавливаемые методом наименьших квадратов.

В данном вычислительном эксперименте: ${{d}_{1}} = 0.76$ мкм, ${{d}_{2}} = 0.75$ мкм, $\eta = 0.244$, ${{f}_{0}} = 1.045 \times {{10}^{{16}}}$ Кл/(м³ с).

На фиг. 1–3 представлены результаты моделирования основных характеристик процесса стационарной зарядки ниобата лития при заданных параметрах вычислительного эксперимента: распределение потенциала $\varphi (x,y,z)$ (фиг. 1), распределение плотности зарядов (фиг. 2) и распределение абсолютного значения напряженности электрического поля, индуцированного инжектированными зарядами (фиг. 3).

Максимальное значение функции $f$ достигается в точках сферы ${{(x - L{\text{/}}2)}^{2}} + {{(y - L{\text{/}}2)}^{2}} + $ $ + \;{{z}^{2}} = \delta _{1}^{2}$ и равно ${{f}_{{{\text{max}}}}} = {{f}_{0}}(1 - \eta ) = 7.9002 \times {{10}^{{16}}}$ Кл/(м³ с). Тогда согласно (3.1) получаем следующую теоретическую оценку максимального значения объемной плотности заряда

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

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

5. ЗАКЛЮЧЕНИЕ

В настоящей работе обоснована корректность реакционно-диффузионно-дрейфовой модели (1.1)–(1.3), позволяющей прогнозировать эффекты последействия электронного облучения на полярные диэлектрические материалы по данным численных оценок объемной плотности распределения зарядов и потенциала в зондируемом материале. Приведены результаты вычислительных экспериментов, подтверждающих теоретически установленный принцип максимума для плотности зарядов $\rho $. При этом мы предполагали, что вся поверхность сегнетоэлектрика покрыта тонким проводящим слоем. В следующих работах будут рассмотрены ситуации, когда только часть границы является проводником, а другая часть остается диэлектриком (см. [33]).

Отметим, что в отношении сегнетоэлектрических материалов самостоятельный научный интерес представляет обратная задача, состоящая в определении условий воздействия электронным зондом, при которых возможны инверсия поляризации и локальное переключение доменов. Указанная и другие обратные задачи, а также задачи управления для модели (1.1)–(1.3) будут исследованы с использованием методов [22], [34], [35] в последующих работах авторов.