Автор: Пользователь скрыл имя, 10 Марта 2013 в 21:45, дипломная работа
В данной работе рассматривается задача с подвижной границей, описывающая поведение испаряющегося слоя двухкомпонентной жидкости. При решении задачи необходимо находить функции концентрации и температуры, связанные нелинейной зависимостью на границе. При этом следует определить положение границы слоя, изменение которого происходит с течением времени в результате испарения растворителя из раствора под действием постоянного теплового потока.
ВВЕДЕНИЕ………………………………………………………………………3
1. ПОСТАНОВКА ЗАДАЧИ…………………………………………………5
1.1. Описание физического процесса ……...……………………………….5
1.2. Математическая постановка……………………………………………6
1.3. Приведение задачи к безразмерному виду……………………………..9
2. РАЗРАБОТКА АЛГОРИТМА ЧИСЛЕННОГО РЕШЕНИЯ…………..12
2.1. Построение системы решения уравнений. Метод сеток……………12
2.1.1. Разностные соотношения для уравнения теплопроводности...13
2.1.2. Разностные соотношения для уравнения диффузии…………..15
2.1.3. Метод прогонки………………………...………………………..16
2.1.4. Аппроксимация уравнения движения границы………………..18
2.1.5. Метод Эйлера…………………………………………………….18
2.2. Алгоритм решения задачи…………………………………………….21
3. Программная реализация предложенного решения задачи……….23
3.1. Разработка программы………………………………………………23
3.2. Тестирование программы…………………………………………...25
3.3. Численный анализ задачи……………………………………………27
Заключение……………………………………………………………………..34
Список используемых источников…..………………………………………35
Приложение……………………………………………………………………..36
Исследуем зависимость полей концентрации и температуры, а также положения границы слоя от параметров αD и Q. Количество шагов по времени возьмем равным 1000, т.е. τ=0,0001, число шагов по пространственной переменной – 500.
Возьмем αD=0.95 и посмотрим график численного решения задачи при значениях теплового потока Q = – 1, Q = –0.5, Q = –0.1
Рис. 3: αD=0.95, Q = – 1
Рис. 4: положение границы во времени при αD=0.95, Q = – 1
Рис. 5: αD=0.95, Q = – 0.5
Рис. 6: положение границы во времени при αD=0.95, Q = – 0.5
Рис. 7: αD=0.95, Q = – 0.1
Рис. 8: положение границы во времени при αD=0.95, Q = – 0.1
Посмотрим на распределение поля концентрации и температуры при αD=0.5 и αD=0.2 при тех же значениях Q.
Рис. 9: αD=0.5, Q = – 1
Рис. 10: αD=0.5, Q = – 0.5
Рис. 11: αD=0.5, Q = – 0.1
Рис. 12: αD=0.2, Q = – 1
Рис. 13: αD=0.2, Q = – 0.5
Рис. 14: αD=0.2, Q = – 0.1
Заметна зависимость положения границы и концентрации растворителя от теплового потока. При более маленьких значениях αD концентрация растворителя сильнее снижается у нагреваемой границы, чем при тех же значениях Q для αD=0.95. При этом снижение концентрации продолжается лишь на небольшую глубину, а оставшийся отрезок до нижней границы она сохраняет начальное значение.
ЗАКЛЮЧЕНИЕ
Проведено аналитическое решение исходной задачи, построены разностные схемы для уравнений теплопроводности и диффузии. Разработан алгоритм решения задачи, который реализован программно. Разработанная программа позволяет изучать влияние различных теплофизических параметров на характер исследуемого процесса диффузии при испарении растворителя из двухкомпонентной смеси жидкости.
СПИСОК ИСПОЛЬЗУЕМЫХ ИСТОЧНИКОВ
6. Самарский А.А., Гулин А.В. Численные методы: Учеб. пособие для
вузов – М.:Наука. Гл. ред. физ. – мат. лит., 1989. – 432с.
ПРИЛОЖЕНИЕ
Unit UnitMain;
interface
uses
Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
Dialogs, Menus, StdCtrls, ExtCtrls, TeeProcs, TeEngine, Chart,
ComCtrls, Grids, Series, Math;
const
xk_min=0.001;
dt_min=0.001;
type
mas=array of Extended; //будет двухразмерный, размерность от пользователя
mas2=array of array of Extended; //будет двухразмерный, размерность от пользователя
TFormMain = class(TForm)
Chart1: TChart;
Series1: TFastLineSeries;
Series2: TFastLineSeries;
GroupBox1: TGroupBox;
Button5: TButton;
Button6: TButton;
Button1: TButton;
GroupBox2: TGroupBox;
Label3: TLabel;
Label4: TLabel;
Label5: TLabel;
Edit1: TEdit;
Edit2: TEdit;
Label6: TLabel;
Label7: TLabel;
Edit3: TEdit;
Label8: TLabel;
Label9: TLabel;
Edit10: TEdit;
Label11: TLabel;
Edit11: TEdit;
Label12: TLabel;
Label13: TLabel;
Edit6: TEdit;
Label14: TLabel;
Edit7: TEdit;
Button7: TButton;
Label15: TLabel;
CheckBox1: TCheckBox;
Label16: TLabel;
Edit5: TEdit;
Label17: TLabel;
Edit12: TEdit;
Edit8: TEdit;
Label18: TLabel;
Edit4: TEdit;
Label19: TLabel;
Label20: TLabel;
Edit13: TEdit;
Series3: TFastLineSeries;
procedure FormCreate(Sender: TObject);
procedure Button1Click(Sender: TObject);
procedure Button5Click(Sender: TObject);
procedure Button6Click(Sender: TObject);
procedure Button7Click(Sender: TObject);
procedure CheckBox1Click(Sender: TObject);
procedure EdsChange(Sender: TObject);
private
{ Private declarations }
procedure OutChart(j:integer); //отображение на графике температуры в j- момент
procedure OutChartG; //отображение на графике положения границы во времени
procedure OutChartS(j:integer);
public
{ Public declarations }
end;
//////////// открытые Переменные модуля
var
FormMain: TFormMain;
x,t:mas;
Tm,Kn:mas2;
S:mas;
dx,dt: Extended; //шаг сетки по координате и по времени
nx,nt,old_nx:integer;
A,aD,Q,R :Extended; //параметры эксперимента
Tm0,Kn0,t0,x0 :Extended; //начальные условия (температура,концентрация,
Tmk,Knk,tk,xk :Extended; //краевые условия (температура,концентрация,
Ypmax, Ypmin: Extended;
//максимальные значения на
boolAll:boolean;
implementation
uses ShellAPI;
{$R *.dfm}
//////////////////////////////
function PointOrComma(s:string):string; //первая запятая заменяется на точку
var p:integer;
ch_from, ch_to:char; //преобразование от ... к ....
begin
ch_from:='.'; ch_to:=','; //можно поменять в зависимости от настройки системы
p:=pos(ch_from,s);
if p>0 then s[p]:=ch_to;
PointOrComma:=s;
end;
procedure SetChartRange(Xmx,Xmn,Ymx,Ymn:
begin //устанавливает диапазон, без ошибок типа Chart1.LeftAxis.Maximum < Chart1.LeftAxis.Minimum
if (Xmx + abs(Xmx-Xmn)*0.05) > formMain.Chart1.BottomAxis.
begin
formMain.Chart1.BottomAxis.
formMain.Chart1.BottomAxis.
end
else
begin
formMain.Chart1.BottomAxis.
formMain.Chart1.BottomAxis.
end;
if (Ymx + abs(Ymx-Ymn)*0.05)> formMain.Chart1.LeftAxis.
begin
formMain.Chart1.LeftAxis.
formMain.Chart1.LeftAxis.
end
else
begin
formMain.Chart1.LeftAxis.
formMain.Chart1.LeftAxis.
end;
end;
/////////////////////// Описание двух
внутренних классов для
type
TprogonK = class // для трехточечной прогонки концентрации
va,vb,vc,F:mas;
//vb-вектор главной диагонали
ka,kb,kc:Extended;
//коэф.при неизвестных ka*Tm[
public
constructor Create; //конструктор для квадратной матрицы
procedure Progon1step(j:integer); //для прогонки по узлам слоя в j-момент времени
end;
//--------------------------
Tprogon = class // для трехточечной прогонки температуры
progonK: TprogonK; //
va,vb,vc,F:mas;
//vb-вектор главной диагонали
ka,kb,kc:Extended; //коэф.при неизвестных ka*Tm[i-1,j] +kb*Tm[i,j] +kc*Tm[i+1,j] = f(x,t) в уравнениях аппроксимации
public
constructor Create; //конструктор для квадратной матрицы
procedure Progon1step(j:integer); //для прогонки по узлам слоя в j-момент времени
procedure Execute;
end;
///////////////////// реализация Tprogon , для температуры
constructor Tprogon.Create; //конструктор потока
var i,j:integer;
begin
inherited Create; //конструктор родителя
old_nx:=0; //слой пока полный, нет исключенных узлов
//шаг сетки по координате и по времени
if xk<=0 then xk:=xk_min; //слой не может быть <=0 ; xk_min заданная константа
dx:=(xk-x0)/nx;
if tk<=0 then dt:=dt_min else dt:=(tk-t0)/nt; //dt_min заданная константа
if dt<=dt_min then begin dt:=dt_min; nt:=round((tk-t0)/dt)+1; end;
setlength(Tm,nx+1,nt+1); //установка размерности массивов
for j:=0 to nt do //j-счетчик по времени
begin
Tm[0,j]:=1;
//изначально температура у
Tm[nx,j]:=Tmk/Tm0; //граничные условия температур
end;
for i:=0 to nx do //i-счетчик по координате
begin
Tm[i,0]:=1; //изначально вся среда (слой) заполняется 1, это равносильно
начальной температуре
end;
setlength(va,nx+1); setlength(vb,nx+1); setlength(vc,nx+1);
setlength(F,nx+1);
setlength(S,nt+1); //толщина слоя - функция от времени
Tm[0,1]:=Tm0/Tm0; Tm[nx,1]:=Tmk/Tm0; //граничные условия температур
F[0]:= Tm[0,1]; F[nx]:= Tm[nx,1]; //свободные члены уравнений матрицы
S[0]:=xk/xk;
//расчет коэффициентов
уравнений на основе
ka:= -(dt*A/dx/dx); //при неявной схеме
kb:= 1-2*ka;
kc:=ka;
progonK:=TprogonK.Create; //создает
объект вспомогательного
end;
procedure Tprogon.Progon1step(j:integer)
var i :integer;
begin
//для j-го момента времени
F[0]:= Tm[0,j]; F[nx]:= Tm[nx,j]; //свободные члены уравнений матрицы
for i:=1 to nx-1 do //i-счетчик по координате
begin
va[i]:=ka; vb[i]:=kb; vc[i]:=kc;
//заполнение разряженной
F[i]:=Tm[i,j-1]; //правые части уравнения помещаю в массив у[i], чтобы после преобразований ничего не перемещать
if i=1 then
begin
F[1]:= Tm[1,j-1]-ka*Tm[0,j]; //это свободный член F для уравнения верхней строки матрицы (левая область или дно)
va[1]:=0;
vb[1]:=1;
vc[1]:=kc/kb;
F[1]:=F[1]/kb;
end;
if i=nx-1 then // обязательно до ni-1, а не до ni
begin