مدل سازی فرایند های شیمیایی و حل عددی معادلات

[حــامد]

پس از مدل کردن یک فرایند یا یک تجهیز به معادلاتی میرسیم که معمولا حل آنها از طریق تحلیلی بسیار مشکل یا کلا مقدور نیست بنابراین برای حل اینگونه مسائل که معمولا معادلات مشتقات جزئی هستند بهتر است از روشهای عددی استفاده نمود من خودم به شخصه علاقه زیادی به نرم افزار matlab دارم چون زبان سطح بالاست کار کردن با آن ساده تر از دیگر زبانهای برنامه نویسی سطح پایین مانند Cو فرترن و... هست و همچین رابط گرافیکی GUI باعث میشه که نتایج خروجی در قالب نمودارهای گرافیکی به نمایش در بیاد

در اینجا با یک مثال ساده از انتقال حرارت شروع میکنم که قابل تعمیم یه دیگر معادلات مانند جرم و مومنتوم میباشد
فرض کنید پس از مدل کردن یک تجهیز مثلا یک راکتور به معادله زیر رسیده ایم:

​

که یک معادله انتقال حرارت یک بعدی است و k ضریب انتقال حرارت هدایتی است برای حل ما احتیاج به دو شرط مرزی و یک شرط اولیه داریم که شرایط مرزی را فرضا یه این شکل در مظر میگیریم :
ساده ترین شرایط مرزی شرایط مرزی دریکله است که مقدار عددی متغیر مثلا دما در مرزها تعیین میشه مثلا ما در اینجا دما را در x=0 برر 15 درجه سانتیگراد و در x=end برابر 25 درجه سانتیگراد میگیریم طول قطع را 10 سانتیمتر و تغییرات دما در قطعه را در زمان 0 تا 4 ثانیه بدست میاوریم و شرط اولیه را هم به شکل زیر در نظر میگیریم:​

[COLOR=#000000][COLOR=#008000][FONT=Courier New][B]init = 20 + 5*sin(x)[/B][/FONT][/COLOR][/COLOR]

​

برای حل از چند روش توسط متلب میشود استفاده نمود:​

1-با استفاده از روشهای عددی
2-pdetool
3-تابع pdepe
4-simulink​

 

[حــامد]

1- روش حل عددی بوسیله تفاضلات متناهی finite-differences

با استفاده از تعریف روش تفاضلات متناهی معادله بالا به معادله زیر تبدیل میشود :

​

که پس از مرتب کردن بر اساس

و s = kDt/(Dx)2 معادله بالا به معادله زیر تبدیل میشود:​

​

که بصورت زیر در متلب وارد میشود:​

[LEFT]function u = heat(k, x, t, init, bdry)
% solve the 1D heat equation on the rectangle described by
% vectors x and t with u(x, t(1)) = init and Dirichlet boundary
% conditions u(x(1), t) = bdry(1), u(x(end), t) = bdry(2).[/LEFT]
 
[LEFT]J = length(t);
N = length(x);
dx = mean(diff(x));
dt = mean(diff(t));
s = k*dt/dx^2;[/LEFT]
 
[LEFT]u = zeros(N,J);
u(1, :) = init;[/LEFT]
 
[LEFT]for n = 1:N-1
u(n+1, 2:J-1) = s*(u(n, 3:J) + u(n, 1:J-2)) + (1 - 2*s)*u(n, 2:J-1);
u(n+1, 1) = bdry(1);
u(n+1, J) = bdry(2);
end[/LEFT]

 

[حــامد]

یه توضیح بدم این کدی که در بالا مشاهده میکنید و باید در m-file کپی کنید و در دایرکتوری جاری save کنید و هر با که خواستید معادلاتی از این قبیل را حل کنید یا دستور زیر call کنید

u = heat(k, x, t, init, bdry)

 

[حــامد]

دستورات زیر را انجام بدهید تا نتیجه مشاهده بشه

[LEFT]tvals = linspace(0, 4, 41);
xvals = linspace(-5, 5, 21);
init = 20 + 5*sin(xvals);
uvals = heat(0.02, tvals, xvals, init, [15 25]);
surf(xvals, tvals, uvals)
xlabel x; ylabel t; zlabel u 
title('Iran-eng by hamed')[/LEFT]

 

[حــامد]

همین مساله رو میتوانیم با دستور Pdepe حلش کنیم در زیر روش حل رو توضیح میدم:
حل معادلات دیفرانسیل پاره ای وابسته به زمان در یک بعد
فرض کنید u تابعی باشد که در معادله زیرصدق کند.

و شرایط مرزی آن به صورت زیر باشد.

این معادله می تواند معادله حاکم بر انتقال حرارت در یک تیغه , استوانه و یا کره باشد؛ در این صورت f عبارت مربوط به شار انتقال حرارت و s مربوط به تولید یا مصرف انرژی می باشد.
برای خل این معادله از دستور pdepe استفاده می شود.

sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan)

m : مشخص کننده هندسه مساله است. 0 برای تیغه، 1 برای سیلندر و 2 برای کره
pdefun : تابعی که معادله را تعریف می کند

[c,f,s] = pdefun(x,t,u,dudx)

c, f, s همان پارامترهای معادله دیفرانسیل پاره ای هستند
icfun : تابعی که شرایط اولیه را تعریف می کند

u = icfun(x)

bcfun : تابعی که شرایط مرزی را بیان می کند

[pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t)

اندیس l مربوط x0 , و اندیس r مربوط به xn (نقاط ابتدایی و انتهایی بردار xmesh )
xmesh : برداری شامل نقاط x1 تا xn
tspan : بردار زمان متناظر با بردار xmesh
مثال
یک لوله استوانه ای را در نظر بگیرید که از وسط با یک غشا به دو نیم تقسیم شده است.در یک طرف این لوله گاز A با فشار 10 بار و در طرف دیگر گاز B وجود دارد.اگر t=0 غشا پاره شود گاز A در B نفوذ می کند تغییرات فشار جزیی گاز A را در طول لوله حساب کنید.
معادله حاکم بر این سیستم به این شکل است

D ضریب نفوذ گاز A در B است و معمولا از مرتبه 1e-5 است.
شرایط اولیه

و شرایط مرزی

با مقایسه معادله3 با معادله 1 می بینیم که
m = 0
c = 1
s = 0
f = 1e-5 * DuDx
پس تابع pdefun به این صورت تعریف می شود.

function [c,f,s] = pdefun0(x,t,u,DuDx)
c = 1;
f = 1e-5*DuDx;
s = 0;
تابع pdeic 
function u0 = pdeic0(x)
if ((x >= 0) & (x <= .5))
u0=10;
elseif ((x >= 0.5) & (x <= 1))
u0=0;
end
و تابع pdebc 
function [pl,ql,pr,qr] = pdebc0(xl,ul,xr,ur,t)
pl = 0;
ql = 100000;
pr = 0;
qr = 100000;

حالا با استفاده از pdepe (برنامه زیر) می توانیم معادله را حل کنیم

 

[حــامد]

ابتدا معادله رو تعریف میکنیم با توجه به شکل روتین معادلات pde :

و معادله اصلی

مشخص میشود :
c=1
m=0
f=k*DuDx
s=0​

پس:​

[/B]
[LEFT][FONT=Courier New][FONT=Courier New][B]function [c,f,s] = pdex(x,t,u,DuDx)[/B][/FONT][/FONT][FONT=Courier New]
[FONT=Courier New][B]c = 1;[/B][/FONT]
[FONT=Courier New][B]f = (0.02)*DuDx;% flux is variable conductivity times u_x[/B][/FONT]
[FONT=Courier New][B]s = 0;[/B][/FONT][/LEFT]
[/FONT][B]

 

[حــامد]

حال باید شرایط مرزی و شرط اولیه مشخص شود:
شرایط مرزی:

شکل کلی بصورت زیر است:

پس:​

[LEFT][FONT=Courier New][FONT=Courier New]function [pl,ql,pr,qr] = pdexbc(xl,ul,xr,ur,t)[/FONT][/FONT]
[FONT=Courier New][FONT=Courier New]% q's are zero since we have Dirichlet conditions[/FONT][/FONT]
[FONT=Courier New][FONT=Courier New]% pl = 0 at the left, pr = 0 at the right endpoint[/FONT][/FONT]
[FONT=Courier New][FONT=Courier New]pl = ul-15;[/FONT][/FONT]
[FONT=Courier New][FONT=Courier New]ql = 0;[/FONT][/FONT]
[FONT=Courier New][FONT=Courier New]pr = ur-25;[/FONT][/FONT][FONT=Courier New]
[FONT=Courier New]qr = 0;[/FONT][/LEFT]
[/FONT]

 

[حــامد]

حال باید شرایط اولیه وارد شود:

[LEFT][FONT=Courier New]function u0 = pdexic(x)[/FONT]
[FONT=Courier New]% initial condition at t = 0[/FONT]
[FONT=Courier New]u0 = 20+5*sin(x);[/FONT][/LEFT]

 

[حــامد]

حال دستورات زیر زا اجرا میکنیم:

[LEFT][FONT=Courier New]%Solves a sample Dirichlet problem for the heat equation in a rod,[/FONT]
[FONT=Courier New]%this time with variable conductivity, 21 mesh points[/FONT]
[FONT=Courier New]m = 0;  %This simply means geometry is linear.[/FONT]
[FONT=Courier New]x = linspace(-5,5,21);[/FONT]
[FONT=Courier New]t = linspace(0,4,81);[/FONT]
[FONT=Courier New]sol = pdepe(m,@pdex,@pdexic,@pdexbc,x,t);[/FONT]
[FONT=Courier New]% Extract the first solution component as u.[/FONT]
[FONT=Courier New]u = sol(:,:,1);[/FONT]
[FONT=Courier New]% A surface plot is often a good way to study a solution.[/FONT]
[FONT=Courier New]surf(x,t,u)    [/FONT]
[FONT=Courier New]title('Numerical solution computed with 21 mesh points in x.')[/FONT]
[FONT=Courier New]xlabel('x'), ylabel('t'), zlabel('u')[/FONT]
[FONT=Courier New]% A solution profile can also be illuminating.[/FONT]
[FONT=Courier New]figure[/FONT]
[FONT=Courier New]plot(x,u(end,:))[/FONT]
[FONT=Courier New]title('Solution at t = 4')[/FONT]
[FONT=Courier New]xlabel('x'), ylabel('u(x,4)')[/FONT][/LEFT]

 

[eghbali66]

سلام حامد آقا. مرسی بازم. واقعاً مثه همیشه بازم گل کاشتی...

 

[msn_syst3m]

آفرین
مهشره
من شخصا خیلی matlab کار می کنم
تا حالا ندیده بودم یه نفر همه چیزو راجه به pde یه جا جمع کنه
آفرین صد آفرین
حامد جان image processing کار کردی ؟
تکنیک Eigenface رو چیزی نوشتی؟

 

[harasani]

حل معادله ديفرانسيل جزئي بيضوي سه بعدي حرارت با شرايط مرزي نوع دوم

حل معادله ديفرانسيل جزئي بيضوي سه بعدي حرارت با شرايط مرزي نوع دوم

با سلام
ميخوام با مطلب برنامه اي براي توزيع دماي يك مكعب (سه بعدي كارتزين) بنويسم كه از دو وجه داراي convection و 4 وجه ديگر عايق هستند. كسي چنين برنامه اي نوشته؟
لطفا راهنمايي و كمك كنيد.

 

[Hermione Granger]

سلام دوستان
من میخوام یه برنامه با متلب بنویسم ولی مشکلم اینه که نمیدونم چطوری باید یه ماتریس تعریف کنم که مثلا" i*j*k متغیر داشته باشه
یعنی در حقیقت آرایه هاش 3 بعد داشته باشه
هدفم نوشتن برنامه حل معادلات دیفرانسیل PDE هست که سه متغیر x , y , t داشته باشند

 

[payam8453]

سلام
من نمیدونم x,t که باید وارد بشه چیه یعنی باید بخش بخش کنم و به صورت بردار وارد کنم یا زمان و مکان کل را به صورت یه عدد وارد کنم
شرط همگرایی روش finite differenceهم در این کد رعایت شده؟

 

[hex6789]

دوست عزیز ابتدا یک ماتریس صفر (مثلاً ماتریس A) با ابعاد 3در3 ایجاد کنید: (A=zeros(3,3,3
حالا به راحتی می‏تونید آرایه های مورد نظرتون رو به صورت A(i,j,k)=x تعریف کنید

 

[nedabadali]

واقعا خوب بود سپاس