Математическое моделирование диффузии элементов в поверхностном слое материала

Размещено на http://www.allbest.ru/

Размещено на http://www.allbest.ru/

Министерство образования и науки РФ

Федеральное государственное бюджетное образовательное учреждение высшего профессионального образования

«Национальный исследовательский

Томский политехнический университет»

Институт кибернетики

Направление - Прикладная математика и информатика

Кафедра - прикладной математики

ПОЯСНИТЕЛЬНАЯ ЗАПИСКА

к курсовой работе на тему: «Математическое моделирование диффузии элементов в поверхностном слое материала»

Выполнил студент гр. 8Б90 Чан Ми Ким Ан

Томск - 2012

Реферат

Ключевые слова: Уравнение диффузии, коэффициент диффузии, эффективная глубина активации, степень активации, pdede (Matlab)

Объект исследования: Железная подложка с медным покрытием .

Предмет исследования: решение дифференциального уравнения диффузии аналитически и с помощью математического пакета Matlab.

Цель работы: Аналитическое и численное решение уравнения диффузии. Математическое моделирование диффузии Cu в поверхностном слое материала (Fe), представление графического материала и анализ результатов численного эксперимента.

Введение

Уравнение диффузии или уравнение теплопроводности представляет собой частный вид дифференциального уравнения в частных производных.

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

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

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

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

1. Нестационарное уравнение диффузии

1.1 Общее определение диффузии. Законы Фика

Диффузия - процесс выравнивания концентрации частиц в среде. При наличии градиента концентрации ЃЮС частиц в веществе возникает поток этих частиц j, выравнивающий их концентрации. Связь между величиной потока и градиентом концентрации выражается первым законом Фика (первое уравнение Фика):

j = -DЃЮC. (1.1.1)

Здесь D - коэффициент диффузии. Знак «минус» означает, что поток направлен из области с большей концентрацией в область с меньшей концентрацией.

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

Коэффициент диффузии D - это константа, зависящая от природы растворителя и растворенного вещества. Его размерность - м²/с. В справочниках часто приводится D в см²/с. Основным источником информации о коэффициентах диффузии в твердых телах, как и в жидкости, является эксперимент.

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

(1.1.2)

где: D₀ - фактор диффузии, называемый также предэкспоненциальным множителем;

E - энергия активации,

R - постоянная Больцмана,

T - температура.

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

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

Изменение концентрации диффундирующего вещества в пространстве и по времени описывает второе уравнение (второй закон) Фика. Оно формируется следующим образом.

По закону сохранения вещества, который свидетельствует о том, что изменение плотности (или концентрации) вещества в данном элементе объема происходит за счет втекания (или вытекания) вещества в этот элемент (уравнение непрерывности, входящее в систему уравнений сплошной среды) имеем:

(1.1.3)

Подставляя в это уравнение (1.1.3) выражение для плотности потока (1.1.1), получим второе уравнение Фика с коэффициентом диффузии, зависящим от координат:

 (1.1.4а)

Если D не зависит от координат, то имеем:

(1.1.4б)

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

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

где -- коэффициент взаимной диффузии.

Если размер кристалла и его плотность в процессе опыта не изменяются, то j'_A = - j'_B и dC_A/dx = - dC_B/dx, тогда коэффициент взаимной диффузии одинаков, входит ли он в j'_A или в j'_в. Очевидно, что величина зависит как от подвижности обоих компонентов, так и от взаимодействия между ними. Поэтому коэффициент взаимной диффузии -- наиболее сложная для интерпретации диффузионная характеристика.

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

Тогда можно записать связь между потоками в неподвижной (лабораторной, j'_i) и движущейся (j'_i) системах координат в виде

где v_k -- скорость «течения» решетки, и определить собственный коэффициент диффузии i-того компонента (D_i) через поток относительно движущегося наблюдателя:

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

В отсутствие внешних сил, если D_A = D_B, то v_k = 0. Тогда = D_A = D_B.

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

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

Помимо собственных коэффициентов, подвижность i-того компонента можно охарактеризовать также парциальными коэффициентами диффузии (D_ik) которые вводят следующим образом:

Парциальные коэффициенты можно определить как для собственной, так и для взаимной диффузии; в последнем случае их надо вводить через поток j'_i, найденный в неподвижной системе координат.

Как следует из определения , надо различать диагональные парциальные коэффициенты (D_ii), которые отражают влияние на поток i-того компонента его же градиента концентрации, от перекрестных парциальных коэффициентов (D_ik), отражающих влияние «чужих» градиентов (dc_k/dx).

Сравнение выражений и позволяет установить связь между собственными и парциальными коэффициентами:

В бинарной системе dC_A/dx = -dC_B/dx и D_A=D_AA-D_AB.

Таким образом, если все D_ik малы, то D_i = D_ii, поэтому часто собственные коэффициенты мало отличаются от парциальных.

Все введенные до сих пор коэффициенты были коэффициентами гетеродиффузии (или химическими коэффициентами диффузии).

Проще для интерпретации случай, когда градиент концентрации отсутствует (c_i = 0), а имеется только градиент изотопной концентрации, концентрации радиоактивных атомов (). В этом случае

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

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

Существуют два вида уравнения диффузии: нестационарное и стационарное.

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

.

1.2 Нестационарное уравнение диффузии

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

а) Одномерный случай

В случае одномерного диффузионного процесса с коэффициентом диффузии (теплопроводности) D уравнение имеет вид:

Где f(x,t) - заданная функция.

При постоянном D уравнение приобретает вид:

где С(x,t) - концентрация диффундирующего вещества, af(x,t)- функция, описывающая источники вещества (тепла).

б) Трёхмерный случай

В трёхмерном случае уравнение приобретает вид:

а при постоянном D приобретает вид:

где-- оператор Лапласа.

в) Решение нестационарного уравнения диффузии в одномерном случае

В одномерном случае фундаментальное решение однородного уравнения (при начальном условии, выражаемом дельта-функцией С(x,0)=д(x) и граничном условии С(?,t)=0) есть

В этом случае С(x,t) можно интерпретировать как плотность вероятности того, что одна частица, находившаяся в начальный момент времени в исходном пункте, через время t перейдёт в пункт с координатой x. Тогда средний квадрат пройденного пути есть

Такому параболическому закону подчиняется, например, рост окисных плёнок на поверхностях металлов при высоких температурах. Необходимым условием при этом является достаточное количество кислорода в окружающей среде.

В случае произвольного начального распределения С(x,0) общее решение уравнения диффузии представляется в интегральном виде как свёртка:



Частные замечания

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

2. Постановка краевых задач о зависимости концентрации вещества от пространственной координаты и времени. Одномерный случай

2.1 Краевая задача без начальных условий

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

Найти решение уравнения диффузии в области

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

где u(x,t) - концентрация частиц в среде.

м₁(t)им₂(t)-- заданные функции.

l - длина стержня.

2.2 Краевая задача с начальными условиями (задача Коши)

а) Краевая задача для бесконечного стержня

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

Найти решение уравнения диффузии в области

удовлетворяющее условию u(x,t₀) = ц(x) ()

где ц(x) -- заданная функция.

б) Краевая задача для полубесконечного стержня

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

Найти решение уравнения диффузии в области

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

где ц(x)и м(t) -- заданные функции.

в) Краевая задача для ограниченного стержня.

В области рассмотрим следующую краевую задачу

- уравнение диффузии.

Если f(x,t)=0, то такое уравнение называют однородным, в противном случае - неоднородным.

u(x,0)=ц(x) (0?x?l) - начальное условие в момент времени t=0, концентрация в точке x задается функцией ц(x)

, 0 ? t ? T - краевые условия.

Функции м₁(t) и м₂(t) задают значение концентрации в граничных точках 0 и l в любой момент времени t.

В зависимости от рода краевых условий, задачи для уравнения диффузии можно разбить на три типа. Рассмотрим общий случай при б_i²+в_i²?0, (i=1,2).

Если б_i=0, (i=1,2), то такое условие называют условием первого рода.

Eсли в_i=0, (i=1,2)- второго рода.

Если б_i и в_i отличны от нуля, то условием третьего рода.

Отсюда получаем задачи для уравнения диффузии - первую, вторую и третью краевую.

3. Решение уравнения диффузии

3.1 Решение уравнения диффузии с однородными граничными условиями методом разделения переменных

а) Однородное уравнение диффузии с однородными граничными условиями

Рассмотрим следующую задачу

, 0 < x < l, 0 < t ? T

, 0 ? x ? l

; 0 ? t ? T

Требуется найти функцию u(x,t)для (x,t): .

Представим искомую функцию в виде произведения u(x,t)=X(x)T(t). Затем предполагаемую форму решения подставим в исходное уравнение, получим

X(x)T?(t)=a²X??(x)T(t).

Разделим выражение на a²X(x)T(t):

, л=const.

Так в левой части уравнения у нас находится функция, зависящая только от t, а в правой - только от x, то фиксируя любое значение x в правой части получаем, что для любого t значение левой части уравнения постоянно. Таким же образом можно убедиться, что и правая часть постоянна, т.е. равна некой константе ? л(минус взят для удобства). Таким образом, мы получаем два обыкновенных линейных дифференциальных уравнения:

Обратим внимание на граничные условия исходной задачи и подставим в них предполагаемый вид уравнения, получим: , откуда X(0) = X(l) = 0(T(t) ? 0, т.к. в противном случае мы имели бы решение u(x,t) = 0, а мы ищем только нетривиальные решения).

С учетом полученных граничных условий мы получаем задачу Штурма-Лиувилля:

Ее решение сводится к решению линейного дифференциального уравнения и рассмотрению трех случаев:

а.1. л < 0

В этом случае общий вид решения будет следующим:

Подставив граничные условия, мы убедимся, что решение будетX(x)?0, а мы ищем только нетривиальные решения, следовательно, этот случай не подходит.

а.2. л = 0

Общий вид решения

X(x) = C₁x+C₂.

Несложно убедиться, что этот вариант нам также не подходит.

а.3. л > 0

Общий вид решения

Подставим граничные условия:

X(0) = C₁ = 0

X(l) = = 0

Т.к. мы ищем только нетривиальные решения, C₂ = 0нам не подходит, следовательно

Отсюда:

C учетом найденных л, выведем общее решение линейного дифференциального уравнения:

Должен получиться ответ: D_n = const.

Теперь все готово для того, чтобы записать решение исходной задачи:

, n = 1, 2…

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

Осталось определить значение константы C (зависящей от n) из начального условия

u(x,0) = ц(x)

Для того, чтобы определить значение C_n, необходимо разложить функцию ц(x) вряд Фурье:

Получаем:

Откуда общее решение:

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

б) Неоднородное уравнение диффузии с однородными граничными условиями

Рассмотрим способ решения неоднородного уравнения:

, 0 < x < l, 0 < t ? T

, 0 ? x ? l

; 0 ? t ? T

Пусть

Тогда, пользуясь очевидным соотношением , перепишем исходное уравнение:



Решим последнее линейное неоднородное уравнение методом вариации постоянной. Сначала найдем общее решение однородного линейного уравнения

В общем решении заменим постоянную D на переменную D(t) и подставим в исходное уравнение.

Из начального условия получаем:

u_n(x,0) = X_n(x)T_n(0) = 0

T_n(0) = 0

С учетом условия для T, получаем:

Так как

то F_n(t), очевидно, является коэффициентом ряда Фурье, и равен

В результате, общая формула такова:

3.2 Решение уравнения диффузии с неоднородными граничными условиями методом разделения переменных

Во многих случаях удается решить уравнение теплопроводности с неоднородными краевыми и начальными условиями

С помощью методов, описанных выше и следующего несложного приема. Представим искомую функцию в виде суммы:



Найдем функцию U(x,t):

Таким образом, исходная задача свелась к следующей:

После того, как мы найдем функцию , искомую функцию найдем по формуле:

Получено решение уравнения диффузии с неоднородными граничными условиями методом разделения переменных.

4. Решение уравнения диффузии в среде Matlab. Функция ядра Matlab для решения одномерных нестационарных уравнений (PDE)

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

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

Синтаксис:

sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan) sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options) sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options,p1,p2...)

Аргументы:

+ m - параметр, соответствующий симметрии проблемы:

m = 0 - “плоская” симметрия (m=0 соответствует плоскопараллельной краевой задаче),

m = 1 цилиндрическая симметрия,

m = 2 - сферическая симметрия.

+ pdefun - функция пользователя, которая определяет коэффициенты PDE.

+ icfun - функция пользователя, которая определяет начальные условия. + bcfun - функция пользователя, которая определяет граничные условия.

+ xmesh - вектор [x0, x1, ..., xn] , задающий координатыxточек, в которых численное решение требуется для каждого значения в tspan (xmesh - массив координат узлов одномерной сетки по пространственной координатеx, см. уравнение (4.1)). Элементы xmesh должны удовлетворить условию x₀ < x₁ < ... < x_n, т.е. должны располагаться в порядке возрастания. Длина xmesh должна быть >= 3.

+ tspan - вектор [t₀, t₁, ..., t_f] , задающий моменты времени, в которые решение PDE требуется для каждого значения в mesh (tspan - массив координат узлов одномерной сетки по временной координатеt, см. уравнение (4.1)). Элементы tspan должны удовлетворить t₀ < t₁ < ... < t_f, т.е. должны располагаться в порядке возрастания. Длина tspan должна быть >= 3.

+ options - некоторые параметры поддерживающего решателя ОДУ, вызываемого из pdepe: RelTol, AbsTol, NormControl, InitialStep, и MaxStep. В большинстве случаев, значения по умолчанию для этих параметров порождают удовлетворительные по точности решения. См. odeset для более детального ознакомления.

+ p1, p2, ... - необязательные параметры, которые нужно передать к pdefun, icfun, и bcfun.

Описание sol = pdepe(m, pdefun, icfun, bcfun, xmesh, tspan) решает краевые задачи для систем параболических и эллиптических PDE, содержащих операции дифференцирования первого порядка по одной пространственной переменной и времени. Обыкновенные дифференциальные уравнения (ОДУ), следующие из дискретизации по пространственной координате, интегрируются, чтобы получить значения приближенного решения в моменты времени, точно установленные в tspan. Функция pdepe возвращает значения решения на сетке, заданной в xmesh.

pdepe решает PDE следующего вида:

(4.1)

PDE или система PDE решается для пространственно-временного интервала t₀?t?t_f и a ? x ? b. Интервал [a,b] должен быть конечен m может быть 0, 1, или 2, что соответствует плоской, цилиндрической, или сферической симметрии. Если m > 0, то a должно быть больше, чем 0.

В уравнении 4.1 f(x,t,u,u/x) - функция “потока” и s(x,t,u,u/x) - функция источника. Частные производные по времени от компонентов искомой величины умножаются на диагональную матрицу c(x,t,u,u/x) (см. 4.1). Диагональные элементы этой матрицы являются или тождественно нулями или положительными значениями. Элемент, который является тождественно нулем, соответствует эллиптическому уравнению, иначе параболическому уравнению. Должно иметься, по крайней мере, одно параболическое уравнение во всей системе уравнений. Элемент, с который соответствует параболическому уравнению, может обращаться в нуль в изолированных значения хx, если эти значения x - узлы сетки. Разрывы распределения c и/или s в заданном интервале пространственных координат разрешаются при условии, что точка сетки представлена в каждом подинтервале.

Для t = t₀и всех значений x компоненты решения удовлетворяют начальным условиям вида

u(x,t₀) = u₀(x) (4.2)

Для всех t, а также x = a или x = b компоненты решения удовлетворяют граничному условию вида

p(x,t,u) + q(x,t)Чf(x,t,u,u/x) = 0 (4.3)

Элементы q должны быть тождественно равны нулю или ни при каких значениях (x, t) не быть равными нулю. Также, из этих двух коэффициентов, только p может зависеть от u.

При вызове функции sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan)

+ m соответствуетm;

+ xmesh(1) и xmesh(end) соответствуют aиb;

+ tspan(1) и tspan(end) соответствуют t₀иt_f;

+ pdefun вычисляет значения коэффициентов PDE c, f, и s (уравнение 4.1). Эта функция вызывается по форме [c,f,s] = pdefun(x,t,u,dudx)

+ Входные параметры - скаляры x и t, и векторы u и dudx, которые приблизительно соответствуют решениюuи его частной производной поx, соответственно.

+ c, f, и s - столбцовые матрицы. c хранит диагональные элементы матрицыc(уравнение 4.1).

+ icfun вычисляет начальные условия. Эта функция вызывается по форме u = icfun(x). Когда вызывается с параметром x, icfun вычисляет и возвращает начальные значения компонентов решения в столбцовые матрице u.

+ bcfun вычисляет значения параметров p и q граничных условий (уравнение 4.3). Эта функция вызывается по форме [pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t)

ul - приближенное решение в границе левой стороны xl = a, и ur - приближенное решение в правой границе xr =b.

pl и ql - столбцовые матрицы, соответствующие p и q, вычисленному в xl, аналогично pr и qr соответствуют xr.

Когда m > 0 и a = 0, ограниченность решения (т.е. конечное значение решения) около x = 0 требует, чтобы поток f обратился в нуль в x = 0. pdepe налагает это граничное условие автоматически и тем самым игнорирует значения, возвращенные в pl и ql.

pdepe возвращает решение в многомерный массив sol.u_i=ui=sol(:,:,i) является приближением к i-му компоненту вектора решения u. Элемент ui(j,k)=sol(j,k,i) является приближением u_i при (t,x) = (tspan(j),xmesh(k)). ui = sol(j,:,i) содержит i-й компонент решения в момент времени tspan(j) и узлах сетки xmesh(:).

sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options, p1, p2, …) передаёт дополнительно параметры p1, p2, ... функциям pdefun, icfun, bcfun. В этом случае функции pdefun, icfun, bcfun могут быть оформлены в виде [c,f,s] = pdefun(x,t,u,dudx,p1,p2,…), u = icfun(x,p1,p2,…), [pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t,p1,p2,…). Можно задать options=[], если параметры режима устанавливать не нужно (см. odeset для более детального ознакомления).

Замечания:

Массивы xmesh и tspan играют различные роли в pdepe.

tspan : функция pdepe выполняет интегрирование по времени с решателем ОДУ, которое выбирает и шаг по времени, и формулу динамически. Элементы tspan просто определяют, в какие моменты времени нужно получить решение, и погрешность слабо зависит от длины tspan.

xmesh : аппроксимации второго порядка к решению сделаны на сетке, точно установленной в mesh. Вообще, лучше использовать близко расположенные узлы (точки) сетки, где решение изменяется быстро. pdepe не выбирает сетку автоматически. Нужно задавать соответствующую фиксированную сетку в xmesh. Погрешность зависит строго от длины xmesh. Когда m > 0, нет необходимости использовать мелкую сетку около x = 0, чтобы правильно смоделировать особенность решения вблизи этой точки.

5. Численные эксперименты. Диффузия алюминия в железную подложку при воздействии импульсным электронным пучком

5.1 Подготовка задачи

Уравнение диффузии

(5.1)

(5.2)

Граничные и начальные условия:

С(0,t)=1; ;

C(x,0)=0.

Выберем примесь алюминия в железе. Нужные численные значения переменных:

T = Tпл = 933 K - температура плавления Al в Fe .

D0 = 3.10^-6 см²/с - предэкспоненциальный множитель.

T = 0,02 c = 2.104 мкс - время счета.

L = 1,5.10^-4 см = 1,5 мкм - расстояние исследования.

- заданное значение.

k = 426 c^-1 - величина характеризует скорость активации и определяется параметрами внешнего воздействия.

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

а)

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

б)

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

+ без учета активации поверхностного слоя

+ с учетом активации поверхностного слоя.

5.2 Решение и анализ результатов

Решение: Обозначаем D₀ = D (для случая а) и (для случая б). Получено новое уравнение диффузии с заданными граничными и начальными условиями

> (5.3)

Решить это уравнение функцией пакета Matlab.

Анализ результатов:

а)

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

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

Получены графики:

Рис.1. Зависимость концентрации Al в Fe от расстояния диффузии в моменте t=0.02c

Рис.2. Зависимость концентрации Al в Fe от времени при x=0.01cм

На рисунках рис.1 и рис.2 представлены аналитическое решение и численное решение, получено из программы Matlab.

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

(5.4)

Результат исследования зависимости концентрация Al в Fe при требовательных условиях.

Рис.3. Концентрация Al в Fe по расстоянию в разных моментах

Рис.4. Концентрация Al в Fe по времени в разных точках.

Из рис.3 и рис.4 видно, что концентрация Al в Fe прямо пропорциональна времени и обратно пропорциональна расстоянию диффузии.

Задача 2:

Исследуем влияние активации среды на концентрацию, и сравнить диффузию при условиях:

+ без учета активации поверхностного слоя

+ с учетом активации поверхностного слоя.

Получены графики:



Рис.5. Концентрация Al в Fe по расстоянию при разных эффективных глубинах активации (t=0.01c).

Рис.6. Концентрация Al в Fе по времени при разных эффективных глубинах активации (x=0.75см).

На рис.5 и рис.6 представлено решение уравнения (5.1) с учетом что . Из рис.5 видно, что при расстояниях меньше чем xw, концентрация Al быстрее уменьшается, чем при расстоянии больше чем xw. А из рис.6 можно сказать, что чем более эффективная глубина активации, тем больше скорость увеличения концентрации Al в Fe. Al на границе подложки (Fe) уменьшается. Хотя на глубине эта зависимость можно измениться. Это можно видеть из рис.5.

Исследование влияния параметра, характеризующего степень активации, на концентрацию Al в Fe.

Получены графики:

Рис.7. Концентрация Al в Fe по расстоянию при разных степенях активации (t=0.01c).

Рис. 8. Концентрация Al в Fе по времени при разных степенях активации (x=0.75cм).

Как хорошо видно на рис.7 и рис.8, при наличии активированной активации распределение концентрации принципиально отличается от случая неактивированного активации. При активированной активации Fe, распределение концентрации Al уменьшается медленнее чем при ненеактивированной активации.

Заключение

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

Приобретены привычки решить дифференциального уравнения и построить графики в Matlab.

диффузия медь концентрация математический

Литература

1. Г.А. Блейхер, В.П. Кривобоков. Теоретические основы обработки материалов импульсными электронными материалами и ионными пучками. Учеб. Пособие. Издательство ТПУ, Томск 2009. - 227 с.

2. Б.С. Бокштейн. Диффузия в металлах. Издательство Металлургия, 119034, Москва, 1978, - 245 с.

3. Муллакаев М.С., Габитов Э.В. Особенности диффузии некоторых переходных металлов в сплавах никеля.

Приложение

Содержание программ в Matlab:

1. Текст файла для

a. нахождения числительных решений,

b. сравнения числительных решений с аналитическими решениями

c. построения графиков решений при требованных условиях.

function exercise_1

D00=3*10^-6; % [см^2/мкс] предэкспоненциальный фактор

D0=D00*exp(-19); % D0=D00*exp(-E/RT), E/RT=19,коэффициент

% объемной диффузии в неактивированном материале;

m = 0; % “плоская” симметрия

n=120; % количество шагов по расстоянию

k=100; % количество шагов по времени

T=20000; % [мкс] длительность воздействия

L=1.5*10^-4; % [см] расстояние

x = linspace(0,L,n); % задающий координаты x точек

t = linspace(0,T,k); % задающий моменты времени

sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t); % решатель

u = sol(:,:,1); % численное решение

% построение массива аналитических решений

for i=1:k

for j=1:n

z(i,j)=exp(-x(j)^2/(4*D0*t(i))); % аналитическое решение

end

end

% построение графиков аналитических и численных решений от расстояния

figure

plot(x,u(end,:),x,z(end,:),'o');

title('Решение при t = 0.02 с');

legend('Численное решение','Аналитическое решение');

xlabel('Расстояние x (см)');

ylabel('C(x,t) (частиц/(cм^2.с))');

% построение графиков численных решений от расстояния в разных t

figure

plot(x,u(10,:),'r',x,u(40,:),'g',x,u(end,:),'b');

title('График зависимости решения от расстояния при разных моментах');

legend('(1) t=0.002 c','(2) t=0,008 c','(3) t=0.02 c');

xlabel('Расстояние x (см)');

ylabel('C(x,t) (частиц/(см^2.с)');

% построение графиков аналитического и численного решений от времени

figure

plot(t.*10^-6,u(:,10),t.*10^-6,z(:,10),'o');

title('Решение при x = 0.125 cм');

legend('Численное решение','Аналитическое решение');

xlabel('Время t (с)');

ylabel('C(x,t) (частиц/(см^2.с)');

% построение графиков численных решений от времени в разных точках

figure

T=t.*10^-6; % переписать t в секундах

plot(T,u(:,4),'g',T,u(:,8),'b',T,u(:,25),'k',T,u(:,40),'r');

title('График зависимости решения от времени при разных растояниях ');

legend('(1) x=0.05cм','(2) x=0.1cм','(3) x=0.3125cм','(4) x=0.5cм');

xlabel('Время t (с)');

ylabel('C(x,t) (частиц/(см^2.с)');

% --------------------------------------------------------------------

function [c,f,s] = pdex1pde(x,t,u,DuDx) % определяет коэффициенты PDE

D00=3*10^-6;

D0=D00*exp(-19);

c = 1/D0;

f = DuDx;

s = 0;

% --------------------------------------------------------------------

function u0 = pdex1ic(x) % определяет начальные условия

u0 = 0;

% --------------------------------------------------------------------

function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t) %граничные условия

L=1.5*10^-4; % расстояние исследования

% ul - приближенное решение в границе левой стороны xl = a,

% ur - приближенное решение в правой границе xr = b.

% pl и ql - столбцовые матрицы, соответствующие p и q,

% вычисленному в xl, аналогично pr и qr соответствуют xr.

pl = ul-1;

ql = 0;

pr = 0;

qr = L;

2. Текст файла для исследования влияния эффективной глубины активации xw на концентрацию:

function xw

m = 0; % “плоская” симметрия

n=120; % количество шагов по расстоянию

k=100; % количество шагов по времени

T=20000; % [мкс] длительность воздействия

L=1.5*10^-4; % [см] расстояние

x = linspace(0,L,n); % задающий координаты x точек

t = linspace(0,T,k); % задающий моменты времени

sol1 = pdepe(m,@pdex1pde1,@pdex1ic,@pdex1bc,x,t); % решатель 1

sol2 = pdepe(m,@pdex1pde2,@pdex1ic,@pdex1bc,x,t); % решатель 2

sol3 = pdepe(m,@pdex1pde3,@pdex1ic,@pdex1bc,x,t); % решатель 3

u1 = sol1(:,:,1); % 1-й массив решения

u2 = sol2(:,:,1); % 2-й массив решения

u3 = sol3(:,:,1); % 3-й массив решения

xlabel('Расстояние x, мкм');

ylabel('C(x,t) (частиц/(cм^2.с))');

hold on;

plot(x*10^4,u1(50,:),'red');

plot(x*10^4,u2(50,:),'blue');

plot(x*10^4,u3(50,:),'m');

legend('(1) xw=10e-5 см ','(2) xw=5*10e-5 см','(3) xw=10e-4 см');

title('Зависимость концентрации от расстояния при разных xw (t=0.01c) ')

figure;

xlabel('Время t, с');

ylabel('C(x,t) (частиц/(cм^2.с))');

hold on;

plot(t.*10^-6,u1(:,50),'red');

plot(t.*10^-6,u2(:,50),'blue');

plot(t.*10^-6,u3(:,50),'m');

legend('(1) xw=10e-5 см ','(2) xw=5*10e-5 см ','(3) xw=10e-4 см');

title('Зависимость концентрации от времени при разных xw (x=0.75см)');

function [c,f,s] = pdex1pde1(x,t,u,DuDx)% определяет коэффициенты PDE

D00=3*10^-6;

D0=D00*exp(-19);

c = 1/D0;

gamma=2.5; % коэффициент, характеризующий чувствительность

% коэффициента диффузии к активации

kw=1000*10^(-6); % tD - характерное время активации поверхностного слоя

xw=10^(-5); % cм - эффективная глубина активации

if x<xw;

fx=1-x./xw; % характеризует активацию поверхностного слоя по глубине

eta=1-exp(-kw*t.*fx); % характеризующий степень активации

% и изменяющийся в пределах

f = DuDx*exp(gamma*eta);

s = 0;

else

f = DuDx;

s = 0;

end;

function [c,f,s] = pdex1pde2(x,t,u,DuDx)% определяет коэффициенты PDE

D00=3*10^-6;

D0=D00*exp(-19);

c = 1/D0;

gamma=2.5; % коэффициент, характеризующий чувствительность

% коэффициента диффузии к активации

kw=1000*10^(-6); % tD - характерное время активации поверхностного слоя

xw=5*10^(-5); % cм - эффективная глубина активации

if x<xw;

fx=1-x./xw; % характеризует активацию поверхностного слоя по глубине

eta=1-exp(-kw*t.*fx); % характеризующий степень активации

% и изменяющийся в пределах

f = DuDx*exp(gamma*eta);

s = 0;

else

f = DuDx;

s = 0;

end;

function [c,f,s] = pdex1pde3(x,t,u,DuDx)% определяет коэффициенты PDE

D00=3*10^-6;

D0=D00*exp(-19);

c = 1/D0;

gamma=2.5; % коэффициент, характеризующий чувствительность

% коэффициента диффузии к активации

kw=1000*10^(-6); % tD - характерное время активации поверхностного слоя

xw=10^(-4); % cм - эффективная глубина активации

if x<xw;

fx=1-x./xw; % характеризует активацию поверхностного слоя по глубине

eta=1-exp(-kw*t.*fx); % характеризующий степень активации

% и изменяющийся в пределах

f = DuDx*exp(gamma*eta);

s = 0;

else

f = DuDx;

s = 0;

end;

% --------------------------------------------------------------------------

function u0 = pdex1ic(x) % определяет начальные условия

u0 = 0;

% --------------------------------------------------------------------------

function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t) % граничные условия

L=1.5*10^-4; % расстояние исследования

% ul - приближенное решение в границе левой стороны xl = a,

% ur - приближенное решение в правой границе xr = b.

% pl и ql - столбцовые матрицы, соответствующие p и q,

% вычисленному в xl, аналогично pr и qr соответствуют xr.

pl = ul-1;

ql = 0;

pr = 0;

qr = L;

3. Текст файла для исследования влияния степени активации на концентрацию.

function eta

m = 0; % “плоская” симметрия

n=120; % количество шагов по растоянию

k=100; % количество шагов по времени

T=20000; % [мкс] длительность воздействия

L=1.5*10^-4; % [см] расстояние

x = linspace(0,L,n); % задающий координаты x точек

t = linspace(0,T,k); % задающий моменты времени

sol1 = pdepe(m,@pdex1pde1,@pdex1ic,@pdex1bc,x,t); % решатель 1

sol2 = pdepe(m,@pdex1pde2,@pdex1ic,@pdex1bc,x,t); % решатель 2

sol3 = pdepe(m,@pdex1pde3,@pdex1ic,@pdex1bc,x,t); % решатель 3

u1 = sol1(:,:,1); % 1-й массив решения

u2 = sol2(:,:,1); % 2-й массив решения

u3 = sol3(:,:,1); % 3-й массив решения

xlabel('Расстояние x, мкм');

ylabel('C(x,t) (частиц/(cм^2.с))');

hold on;

plot(x*10^4,u1(60,:),'red');

plot(x*10^4,u2(60,:),'blue');

plot(x*10^4,u3(60,:),'green');

legend('(1) - eta=0','(2) - eta=0.5','(3) - eta=1');

title('Зависимость концентрации от расстояния при разных eta (t=0.01c)');

figure;

xlabel('Время t, с');

ylabel('C(x,t) (частиц/(cм^2.с))');

hold on;

plot(t.*10^-6,u1(:,50),'red');

plot(t.*10^-6,u2(:,50),'blue');

plot(t.*10^-6,u3(:,50),'g');

legend('(1) - eta=0','(2) - eta=0.5','(3) - eta=1');

title('Зависимость концентрации от времени при разных eta (x=0.75см)');

function [c,f,s] = pdex1pde1(x,t,u,DuDx)% определяет коэффициенты PDE

D00=3*10^-6; % [см^2/мкс] предэкспоненциальный фактор

D0=D00*exp(-19); % D0=D00*exp(-E/RT), E/RT=19,коэффициент

% объемной диффузии в неактивированном материале

gamma=2.5; % коэффициент, характеризующий чувствительность

% коэффициента диффузии к активации

eta1 = 0; % параметр, характеризующий степень активации (неактивированная)

D01=D0*exp(gamma*eta1);

c = 1/D01;

kw=1000*10^(-6); % tD - характерное время активации поверхностного слоя

xw=10^(-5); % cм - эффективная глубина активации

if x<xw;

fx=1-x./xw; % характеризует активацию поверхностного слоя по глубине

eta=1-exp(-kw*t.*fx); % характеризующий степень активации

% и изменяющийся в пределах

f = DuDx*exp(gamma*eta);

s = 0;

else

f = DuDx;

s = 0;

end;

function [c,f,s] = pdex1pde2(x,t,u,DuDx)% определяет коэффициенты PDE

D00=3*10^-6; % [см^2/мкс] предэкспоненциальный фактор

D0=D00*exp(-19); % D0=D00*exp(-E/RT), E/RT=19,коэффициент

% объемной диффузии в неактивированном материале

gamma=2.5; % коэффициент, характеризующий чувствительность

% коэффициента диффузии к активации

eta2 = 0.5; % параметр, характеризующий степень активации (средно-активированная)

D02=D0*exp(gamma*eta2);

c = 1/D02;

kw=1000*10^(-6); % tD - характерное время активации поверхностного слоя

xw=10^(-5); % cм - эффективная глубина активации

if x<xw;

fx=1-x./xw; % характеризует активацию поверхностного слоя по глубине

eta=1-exp(-kw*t.*fx); % характеризующий степень активации

% и изменяющийся в пределах

f = DuDx*exp(gamma*eta);

s = 0;

else

f = DuDx;

s = 0;

end;

function [c,f,s] = pdex1pde3(x,t,u,DuDx)% определяет коэффициенты PDE

D00=3*10^-6; % [см^2/мкс] предэкспоненциальный фактор

D0=D00*exp(-19); % D0=D00*exp(-E/RT), E/RT=19,коэффициент

% объемной диффузии в неактивированном материале

gamma=2.5; % коэффициент, характеризующий чувствительность

% коэффициента диффузии к активации

eta3 = 1; % параметр, характеризующий степень активации (активированная поверхность)

D03=D0*exp(gamma*eta3);

c = 1/D03;

kw=1000*10^(-6); % tD - характерное время активации поверхностного слоя

xw=10^(-5); % cм - эффективная глубина активации

if x<xw;

fx=1-x./xw; % характеризует активацию поверхностного слоя по глубине

eta=1-exp(-kw*t.*fx); % характеризующий степень активации

% и изменяющийся в пределах

f = DuDx*exp(gamma*eta);

s = 0;

else

f = DuDx;

s = 0;

end;

% --------------------------------------------------------------------------

function u0 = pdex1ic(x) % определяет начальные условия

u0 = 0;

% --------------------------------------------------------------------------

function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t) % определяет граничные условия

L=1.5*10^-4; % расстояние исследования

% ul - приближенное решение в границе левой стороны xl = a,

% ur - приближенное решение в правой границе xr = b.

% pl и ql - столбцовые матрицы, соответствующие p и q,

% вычисленному в xl, аналогично pr и qr соответствуют xr.

pl = ul-1;

ql = 0;

pr = 0;

qr = L;

Размещено на Allbest.ru