کدمتلب CFD یا دینامیک سیالات محاسباتی

بسمه تعالی

تمارین درس CFD به صورت برنامه نویسی با MATLAB

فهرست مطالب:
نفوذ حالت دائم یک بعدی ………………………………………………………………………………………….3

تمرین1و2………………………………………………………………………………………………………….4
الگوریتم تمرین 1و2……………………………………………………………………………………………….4
متن برنامه تمرین 1و2………………………………………………………………………..5
نتایج و نمودارهای تمرین 1…………………………………………………………………………………………6
نتایج و نمودارهای تمرین 2…………………………………………………………………………………………12
تمرین 3……………………………………………………………………………………………………………..16
الگوریتم تمرین 3…………………………………………………………………………………………………….17
متن برنامه تمرین 3………………………………………………………………………………………………….17
خروجی برنامه تمرین 3……………………………………………………………………………………………..18

جابجایی و نفوذ یک بعدی دائم………………………………………………………………………………………..20

طرح اختلاف مرکزی…………………………………………………………………………………………………20
طرح اختلاف بالادست…………………………………………………………………………………………………21
طرح نمایی…………………………………………………………………………………………………………….21

جواب تمرین 4 به روش اختلاف مرکزی………………………………………………………………………………22
الگوریتم………………………………………………………………………………………………………………..22
متن برنامه………………………………………………………………………………………………………………22
خروجی برنامه………………………………………………………………………………………………………….24

جواب تمرین 4 به روش اختلاف بالادست………………………………………………………………………………26

الگوریتم…………………………………………………………………………………………………………………26
متن برنامه……………………………………………………………………………………………………………….26
خروجی برنامه…………………………………………………………………………………………………………..27

جواب تمرین 4 به روش طرح نمایی……………………………………………………………………………………..29
الگوریتم…………………………………………………………………………………………………………………..29
متن برنامه…………………………………………………………………………………………………………………29
خروجی برنامه…………………………………………………………………………………………………………….30

نفوذ حالت دائم یک بعدی
در حالت کلی معادله ی حرکت به صورت زیر است:
∂/∂t (ρ∅)+div(ρu∅)=div(Γgrad∅)+S (1)
که از چپ به ترتیب از چهار ترم:
وابسته به زمان، جابجایی، پخش و چشمه تشکیل شده است.
∅ می تواند هر خاصیتی همانند دما یا سرعت را شامل شود. که در این گزارش ما از خاصیت دما استفاده کرده ایم.
Γضریب پخش بوده که در صورت در مباحث مربوط به خاصیت دما تبدیل به k به عنوان ضریب هدایت حرارتی خواهد شد.
S جمله ی چشمه است.
در حالت نفوذ دائم یک بعدی به دلیل دائم بودن و نیز نبود ترم جابجایی طرف اول معادله حذف خواهد شد و با توجه به یک بعدی بودن، معادله به صورت ساده شده ی زیر تبدیل می گردد:
d/dx (Γ (d∅)/dx )+S=0 (2)
اما استفاده از معادله ی فوق در روش حجم محدود مشروط به تبدیل آن به حالت معادله ی انفصال می باشد که در زیر تشریح خواهد شد.
برای به دست آوردن معادله ی انفصال در روش حجم محدود ابتدا با به کار بردن دسته گره هایی اقدام به تولید شبکه می نماییم. به صورت شکل 1.
〖(δx)〗_e 〖(δx)〗_w

شکل (1)
خط کمرنگ نمایانگر وجوه حجم کنترل است که با e و w نمایش داده شده. به دلیل یک بعدی بودن فرض را بر ضخامت واحد در امتداد های y و z می گیریم.
در نهایت اگر از معادله ی (2) روی حجم کنترل انتگرال بگیریم خواهیم داشت:
(k dT/dx)_e-(k dT/dx)_w+∫_w^e▒〖S dx=0 (3)〗
در رابطه ی (3) به منظور محاسبه ی گرادیان dT/dx و ضریب نفوذ k در سطح مشترک به نظر می رسد که تقریب های خطی واضح بوده و ساده ترین راه برای محاسبه ی آنها می باشد. پس معادله ی زیر به دست خواهد آمد:
(k_e (T_E-T_P))/((δx_e))-(k_w (T_P-T_W ))/((δx_w ) )+S ̅∆x=0 (4)
S ̅∆x=S_u+S_P T_P (5)
با جایگذاری رابطه ی (5) در رابطه ی (4) و با توجه به اینکه ازین پس ضرایب T_W, T_E, T_P را به ترتیب با a_W,a_E,a_P نشان خواهیم داد بنابراین معادله ی (4) به صورت زیر نوشته می شود:
a_P T_P=a_W T_W+a_E T_E+S_u (6)
a_P=a_W+a_E-S_P (7)
a_W=k_w/(δx_w ) , a_E=k_e/(δx_e ) (8)
معادله ی انفصال (6) باید در هر یک از نقاط گرهی به منظور حل مساله به کار رود. برای حجم های کنترلی که در همسایگی مرزهای ناحیه قرار دارند معادله ی گسسته کلی (6) جهت وارد کردن شرایط مرزی اصلاح می شود. در نهایت دستگاه معادلات جبری برای تعیین توزیع T حل می شود.
تمرین1-
مساله ی هدایت حرارتی چشمه – آزاد در یک میله عایق با دماهای ثابت ابتدایی و انتهایی به ترتیب 〖100〗^℃ و 〖500〗^℃ را در نظر بگیرید. توزیع درجه حرارت حالت دائم میله را محاسبه کنید.
ضریب هدایت حرارتی k برابر 1000 w/(m.k) و سطح مقطع A برابر 〖10〗^(-2) m^2 می باشد.

تمرین 2-
یک صفحه ی بزرگ به ضخامت L=2cm با هدایت حرارتی ثابتw/(m.k) k=0.5 و تولید گرمای یکنواخت q=1000 kw/m^3 را نشان می دهد. سطوح A و B به ترتیب در درجه حرارت 〖100〗^℃ و 〖200〗^℃ قرار دارند. فرض می کنیم که ابعاد آن در جهت های y و z بزرگ باشد، بنابراین گرادیان دما فقط در راستای x مورد توجه می باشد. مطلوبست توزیع درجه حرارت دائم و مقایسه ی نتایج عددی با حل تحلیلی.

جواب: الگوریتم حل هر دو مساله یکسان بوده و تنها تفاوت آنها در وجود چشمه ی یکنواخت در طول میله می باشد.

الگوریتم حل:
1-وارد کردن اطلاعات اولیه شامل :
ضریب هدایت حرارتی (k) – دما در ابتدا و انتهای میله (T_A) و(T_B) – طول میله (L)- مساحت سطح مقطع (A)- چشمه ی یکنواخت (q)
2- شبکه بندی: وارد نمودن تعداد گره ها و موقعیت هر گره نسبت به مبدا.
3-محاسبه ی ضرایب a_W,a_E,a_P و S_u , S_P در هرگره با استفاده از روابط (6) تا (8). لازم به ذکر است که گره های مجاور مرز همواره شامل عبارت های S_u , S_P بوده و گره های داخلی فاقد این عبارت ها هستند مگر در صورت وجود چشمه در این نقاط که در اینصورت شامل عبارت S_u خواهند بود (همانند تمرین 2)
4- جایگذاری ضرایب فوق در رابطه ی (6) و تشکیل دستگاه معادلات جبری و حل آن به منظور یافتن دما در تمامی نقاط گرهی.
5- وارد کردن شرایط مرزی (وارد کردن دما در مرزها)
6-وارد کردن معادله ی حل تحلیلی و محاسبه ی مقادیر T در هرگره.
7-رسم نمودارهای حل عددی و تحلیلی و مقایسه ی آنها با یکدیگر.

متن برنامه:

% % importing problem information
k=input('k:');
TA=input('TA:');TB=input('TB:');
L=input('length:')
A=input('area:');
q=input('q:');

% % reticulation
n=input('nodes:');
nn=n+2;
dx=L/n;
x(1)=0; x(nn)=L;
for j=2:nn-1
x(j)=x(j-1)+dx;
if j==2
x(j)=dx/2;
end
end

% % calculation of Temprature in each node by analayctical equation
ANAeq=(((TB-TA)/L)+(q/(2*k)).*(L-x)).*x+TA;

% % calculation of cofficients in each node:a_W,a_E,a_P,S_u , S_P
Z=k*A/dx;
U=q*A*dx;
for i=1:n
if i==1
SP(i)=-2*Z; Su(i)=2*Z*TA+U;
aW(i)=0; aE(i)=Z;
elseif i==n
aW(i)=Z; aE(i)=0;
SP(i)=-2*Z; Su(i)=2*Z*TB+U;
else
SP(i)=0; Su(i)=U;
aW(i)=Z; aE(i)=Z;
end
aP(i)=aW(i)+aE(i)-SP(i);
end

% % Formation of Algebraic equtions and solving it
T=sym(zeros(1,n));
for k=1:n
T(k)=sym(sprintf('T%d',k));
end
tt=T;
for i=1:n
if i==1
T(i)=(aE(i)*T(i+1)+Su(i))/aP(i);
elseif i==n
T(i)=(aW(i)*T(i-1)+Su(i))/aP(i);
else
T(i)=(aW(i)*T(i-1)+aE(i)*T(i+1)+Su(i))/aP(i);
end
end
TT=T-tt;
ans=solve(TT);
ANStemp=structfun(@subs,ans);

% % Importing boundary condition
therm=ones(1,nn); therm(1)=TA; therm(nn)=TB;
therm(2:1:nn-1)=ANStemp;

% % plotting results
plot(x,therm,'k','LineWidth',2)
hold on
plot(x,ANAeq,'rs','LineWidth',2)
legend('Numerical solution','Exact solution')
xlabel('distance')
ylabel('temperature')
title('A comparison of numerical and analaytical solution')
hold off;

نتایج و نمودارهای مربوط به تمرین 1:
ابتدا با وارد نمودن اطلاعاتی که در کتاب امده نتایج بدست آمده را با نتایج کتاب مقایسه می کنیم تا از صحت برنامه اطمینان حاصل شود. سپس داده های k, A, Δx را به تنهایی و با ثابت نگه داشتن سایر پارامترها مطابق صورت مساله، تغییر داده و بر روی نتایج بحث خواهیم کرد.
نتایج بدست آمده به شرح زیر می باشد.
[█(T_1@T_2@T_3@T_4@T_5@T_6@T_7 )]=[█(100@140@220@300@380@460@500)]

شکل(2)
نتایج بدست آمده کاملا با نتایج موجود در کتاب منطبق است فلذا صحت برنامه تایید می شود.

تغییر شبکه بندی: با n=3 وn=5 و n=25 ب.به بررسی نتایج می پردازیم.
شبکه بندی درشت:
شبکه بندی n=5 در صفحه قبل بررسی شد.
n=3

شکل(3)

شبکه بندی ریز:
n=25

شکل(4)

n=100

شکل(4)

تغییر در مقدار k :
در صورت تمرین مقدار داده شده برای k برابر 1000 w/(m.k) می باشد ما نتایج مقادیر 100,500,2000 را مورد بررسی قرار خواهیم داد.

k=100

شکل (5)

k=500

شکل (6)

k=2000

شکل (7)

تغییر در A:
مقدار A در مساله 〖10〗^(-2) m^2 داده شده است. ما با داده های A=1, 0.1,0.005 به بررسی نتایج می پردازیم.
A=5

شکل (8)

A=1

شکل (9)

A=0.005

شکل (10)

نتیجه گیری:
تغییر در شبکه بندی و پارامترهای k و A تاثیر قابل ملاحظه ای روی نتایج نداشت.

نتایج و نمودارهای مربوط به تمرین 2:
مراحل مربوط به تمرین 1 را عینا در اینجا تکرار می کنیم و به بررسی نتایج می پردازیم.
نتایج حاصل از داده های صورت مساله:
[█(T_1@T_2@T_3@T_4@T_5@T_6@T_7 )]=[█(100@150@218@254@258@230@200)]

شکل (11)
نتایج و نمودار کاملا با کتاب منطبق است. بنابراین اطمینان کامل از صحت برنامه داریم. حال به تغییر تک تک پارامترها و بررسی نتایج می پردازیم.
تغییر شبکه بندی: با n=3 وn=5 و n=25 ب.به بررسی نتایج می پردازیم.
n=5 قبلا بررسی شد.
n=3

شکل (12)
n=9

شکل (13)
n=10

شکل (14)
n=25

شکل (15)
تغییر در k:
مقدار k در این مساله 0.5 لحاظ شده که ما بررسی ها را برای مقادیر 0.1,5 انجام می دهیم.
k=0.01

شکل (16)

k=0.1

شکل (17)

k=5

شکل (18)

k=15

شکل (19)
نتیجه گیری:
برخلاف مثال قبلی تغییر در پارامترهای k و شبکه بندی تاثیر قابل ملاحظه ای بر روی نتایج داشت.
هرچقدر میزان k افزایش می یابد توزیع دما به حالت خطی نزدیک تر می شود.
درست است که شبکه بندی ریز نتایج دقیق تری داد اما از یک میزان به بعد نتایج اختلاف شدیدی با حل دقیق پیدا کرد.

تمرین 3-
در این تمرین، خنک کاری یه پره دایره ای بصورت انتقال حرارت جابجایی در امتداد طولی بررسی می شود. جابجایی بصورت افت گرما یا بخش چشمه ی وابسته به درجه حرارت در معادله ظاهر می شود.
شکل(20) یک پره ی استوانه ای با سطح مقطع یکنواخت A را نشان می دهد. پایه ی پره در درجه حرارت T_B=〖100〗^℃ قرار دارد و انتهایش عایق می باشد. درجه حرارت محیط اطراف 〖20〗^℃ می باشد. انتقال حرارت یک بعدی در این حالت از معادله ی زیر پیروی می کند.
d/dx (kA dT/dx )-hp(T-T_∞ )=0 (9)
در اینجا h ضریب انتقال حرارت جابجایی،p محیط سطح مقطع پره، k قابلیت هدایت حرارتی ماده و T_∞ درجه حرارت محیط اطراف است. توزیع درجه حرارت را در طول پره محاسبه کنید و نتایج را با حل تحلیلی که توسط رابطه ی زیر داده می شود مقایسه کنید.
(T-T_∞)/(T_B-T_∞ )=cosh⁡[n(L-x)]/cosh⁡(nL) (10)

شکل(20)
n^2=hp/kA=25 m^(-2) , L=1m)

حل: معادله ی حاکم در این تمرین شامل یک بخش چاه حرارتی(اتلاف حرارت) بصورت جابجایی بوده و تابعی از درجه حرارت می باشد.
با در نظر گرفتن n^2=hp⁄kA و انتگرال گیری از رابطه ی (9) خواهیم داشت:
[(A dT/dx)_e-(A dT/dx)_w ]-[n^2 (T_P-T_∞ )Aδx]=0 (11)
با اعمال یک تقریب خطی معمولی برای گرادیان درجه حرارت خواهیم داشت:
[((T_E-T_P)/δx)_ -((T_P-T_w)/δx)_ ]-[n^2 (T_P-T_∞ )δx]=0 (12)
این رابطه را بصورت زیر باز نویسی می کنیم:
(1/δx+1/δx) T_P=(1/δx) T_W+(1/δx) T_E+n^2 δxT_∞-n^2 δxT_(P ) (13)
در نتیجه با توجه به رابطه کلی (7) و معادله ی (13) ضرایب a_W,a_E,a_P وS_u , S_P برای گره های داخلی به صورت زیر استخراج می شود:
a_W=1/δx ,a_E=1/δx ,a_P=a_W+a_E-S_P (14)
S_P=-n^2 δx , S_u=n^2 δx T_∞ (15)
با اعمال شرایط مرزی در رابطه ی(13) برای گره های مرزی نیز خواهیم داشت:
گره مرزی غربی:
a_W=0 ,a_E=1/δx ,a_P=a_W+a_E-S_P (16)
S_P=-n^2 δx- 2/δx , S_u=n^2 δx T_∞+2/δx T_B (17)
گره مرزی شرقی: با توجه به عایق بودن این گره ترم چشمه فقط ناشی از چشمه ی ناشی از افت گرما خواهد بود.
a_W=1/δx ,a_E=0 ,a_P=a_W+a_E-S_P (18)
S_P=-n^2 δx , S_u=n^2 δx T_∞ (19)
الگوریتم حل:
1-وارد کردن اطلاعات اولیه شامل :
دما در ابتدا و انتهای میله (T_A) و(T_B) و دمای پیرامون(T_∞)- طول میله (L)- پارامتر مربوط به چاه (n^2=hp/kA)
2- شبکه بندی: وارد نمودن تعداد گره ها و موقعیت هر گره نسبت به مبدا.
3- وارد کردن معادله ی حل تحلیلی و محاسبه ی مقادیر T در هرگره.
4-محاسبه ی ضرایب a_W,a_E,a_P و S_u , S_P در هرگره با استفاده از روابط (14) تا (19).
5- جایگذاری ضرایب فوق در رابطه ی (6) و تشکیل دستگاه معادلات جبری و حل آن به منظور یافتن دما در تمامی نقاط گرهی.
6- وارد کردن شرایط مرزی (وارد کردن دما در مرزها)
7-رسم نمودارهای حل عددی و تحلیلی و مقایسه ی آنها با یکدیگر.

متن برنامه:
% % importing problem information
TA=input('TA:');TB=input('TB:');TO=input('TO:');
L=1;
n2=25; % % n2=(h*p/k*A)

% % reticulation
n=input('nodes:');
nn=n+2;
dx=L/n;
x(1)=0; x(2)=dx/2; x(nn)=L;
Z=1/dx; U=n2*dx;
for j=3:nn-1
x(j)=x(j-1)+dx
end

% % calculation of Temprature in each node by analayctical equation
EXsolution=((TB-TO).*(cosh(sqrt(n2).*(L-x))/(cosh(sqrt(n2).*L))))+TO;

% % calculation of cofficients in each node:a_W,a_E,a_P,S_u , S_P
for i=1:n
if i==1
SP(i)=-2*Z-U; Su(i)=2*Z*TB+U*TO;
aW(i)=0; aE(i)=Z;
elseif i==n
aW(i)=Z; aE(i)=0;
SP(i)=-U; Su(i)=U*TO;
else
SP(i)=-U; Su(i)=U*TO;
aW(i)=Z; aE(i)=Z;
end
aP(i)=aW(i)+aE(i)-SP(i);
end

% % Formation of Algebraic equtions and solving it
aW(1)=[]; aW(n)=0;
aE=flipud(aE')
a=zeros(n,3);
a(:,1)=-aW;
a(:,2)=aP;
a(:,3)=-aE;
K = full(spdiags(a, [-1 0 1], n, n));
K(n,n)=aP(n);
T=Su/K';

% % Importing boundary condition
therm=ones(1,nn); therm(1)=TB; therm(nn)=TA;
therm(2:1:nn-1)=T;

% % plotting results
plot(x,therm,'–rs','LineWidth',2)
hold on
plot(x,EXsolution,'k','LineWidth',2)
legend('Numerical solution','Exact solution')
xlabel('distance')
ylabel('temperature')
title('A comparison of numerical and analaytical solution')
axis([0 L 20 100])
hold off;
خروجی برنامه:
خروجی را برای سه حالت n=5, 10, 40 بررسی می کنیم.
n=5حالت
[█(T_1@T_2@T_3@T_4@T_5@T_6@T_7 )]=[█(100@64.22@36.91@26.50@22.60@21.30@0)]

شکل(21)

:n=10حالت

شکل(22)

:n=40حالت

شکل(23)
نتیجه گیری:
ملاحظه می شود با ریز شدن شبکه حل عددی انطباق بهتری با حل دقیق پیدا می کند.

جابجایی و نفوذ یک بعدی دائم:
در غیاب چشمه ها، نفوذ و جابجایی طبیعی ∅ در یک میدان جریان یک بعدی u از رابطه زیر پیروی می کند:
d/dx (ρu∅)=d/dx (Γ (d∅)/dx) (20)
و همچنین جریان بایستی که پیوستگی را نیز ارضا کند.
d/dx (ρu)=0 (21)
〖(δx)〗_e 〖(δx)〗_w

∆x
شکل(24)

انتگرال معادلات (20) و(21) روی حجم کنترل به روابط زیر منجر می شود:
(ρuA∅)_e-(ρuA∅)_w=(ΓA (∂∅)/∂x)_e-(ΓA (∂∅)/∂x)_w (22)
(ρu∅)_e-(ρu∅)_w=0 (23)
برای نشان دادن شار جرم جابجا شده در واحد سطح و قابلیت نفوذ در سطوح سلول استفاده می شود: F, Dاز دو متغیر
F= ρu , D=Γ/dx (24)
در نهایت با فرض A_w=A_e=A از روش اختلاف مرکزی برای نشان دادن سهم عبارتهای نفوذ در سمت راست استفاده می کنیم:
F_e ∅_e-F_w ∅_w=D_e (∅_E-∅_P )-D_w (∅_P-∅_W ) (25)
F_e-F_w=0 (26)
داریم. برای این منظور از سه طرح زیر در این گزارشe,w منتقل شده از سطوح ∅برای حل معادله ی(25) نیاز به محاسبه ی خاصیت
استفاده شده است:
(Exponential)طرح نمایی – (upwind)- طرح الگوی بالادست (CD) طرح اختلاف مرکزی

(CD) طرح اختلاف مرکزی
تقریب اختلاف مرکزی برای بیان عبارت های نفوذ که در سمت راست معادله ی(25) ظاهر می شوند استفاده می شود و برای محاسبه ی مقادیر جابجایی سطح سلول در سمت چپ این معادله استفاده از یک درونیابی خطی منطقی به نظر می رسد. برای یک شبکه ی را بصورت زیر می نویسیم. ∅یکنواخت مقادیر مربوط به خاصیت
∅_e=(∅_P+∅_E)/2 , ∅_w=(∅_W+∅_P)/2 (27)
جایگذاری نتایج بالا در رابطه ی (25) نتیجه میدهد:
F_e/2 (∅_P+∅_E )-F_w/2 (∅_W+∅_P )=D_e (∅_E-∅_P )-D_w (∅_P-∅_W ) (28)

نتایج اختلاف مرکزی برای معادله ی انفصال جابجایی بصورت زیر بدست می آید:
a_P ϕ_P=a_W ϕ_W+a_E ϕ_E (29)
a_W=D_w+F_w/2 , a_E= D_e- F_e/2 , a_P= a_W+a_E+( F_e-F_w ) (30)

طرح اختلاف بالادست (Upwind)
یکی از بزرگترین نقایص طرح اختلاف مرکزی، ناتوانی آن در تعیین اثر جهت جریان می باشد. در روش اختلاف بالادست این نقص برطرف شده بطوریکه مقدار جابجا شده ی ∅ در سطح سلول برابر مقدار آن در گره بالا دست در نظر گرفته می شود.
با فرض اینکه جهت جریان از غرب به شرق می باشد خواهیم داشت:
∅_w=∅_(W ) , ∅_e=∅_P (31)
با جایگذاری رابطه ی(31) در (25) خواهیم داشت:
[(D_w+F_w )+D_e+(F_e-F_w )] ∅_P=(D_w+F_w ) ∅_w+D_e ∅_E (32)
a_P ϕ_P=a_W ϕ_W+a_E ϕ_E
a_W=D_w+F_w , a_E= D_e , a_P= a_W+a_E+( F_e-F_w ) (33)

طرح نمایی (Exponential)
شار کلی J در این روش به صورت زیر است:
J=ρu∅-Γ dϕ/dx (33)
از رابطه ی(20) و(33) نتیجه می شود:
dJ/dx=0 (34)
که وقتی روی حجم کنترل نشان داده شده در شکل (1) انتگرال گرفته شود می دهد:
J_e-J_w=0 (35)
از طرفی با در نظر گرفتن شرایط مرزی معادله ی (20) می تواند بصورت دقیق حل شود که خواهیم داشت:
(ϕ-ϕ_0)/(∅_L-∅_0 )=(exp(Px⁄L)-1)/(exp(P)-1) , P=ρuL/Γ (36)
حال می توان از جواب دقیق (36) به عنوان یک پروفیل بین نقاط P و E استفاده کرد. وارد کردن این پروفیل در رابطه ی (33) رابطه ی مربوط به J_e را می دهد.
J_e=F_e (ϕ_P+(∅_P-ϕ_E)/(exp(P_e )-1)) (37)
که در آن:
P_e=((ρu)_e 〖 (δx)〗_e)/Γ_e =F_e/D_e (38)
F_e و D_e توسط روابط (24) تعریف می شوند.
سرانجام، با جایگزینی از معادله ی( 37) و به کار بردن عبارت مشابهی برای J_w در معادله ی (35) خواهیم داشت:
F_e (ϕ_P+(ϕ_P-ϕ_E)/(exp(P_e )-1))-F_w (ϕ_W+(ϕ_W-ϕ_P)/(exp(P_w )-1))=0 (39)
در نهایت روابط زیر را برای حالت نمایی خواهیم داشت:
a_P ϕ_P=a_W ϕ_W+a_E ϕ_E
a_E=F_e/(exp(P_e )-1) (40-1)
a_W=(F_w exp(P_w ))/(exp(P_w )-1) (40-2)
a_P= a_W+a_E+( F_e-F_w ) (40-3)
تمرین4-
مساله ی زیر را برای سه حالت: C-D, Upwind, Exponential بررسی کنید.
یک خاصیت ϕ بصورت نفوذ و جابجایی از یک محدوده یک بعدی که در شکل(25) ترسیم شده است منتقل ی شود. معادله ی حاکم همان معادله (20) می باشد. شرایط مرزی عبارتند از: در x=0 داریم ϕ_0=1 و در x=L داریم ϕ_L=0. از پنج سلول با اندازه مساوی و نیز از روش اختلاف مرکزی، برای نفوذ و جابجایی استفاده می شود. ϕ را برای حالتهای زیر بصورت تابعی از x محاسبه کنید.
الف- u=0.1m⁄s
ب-u=2.5 m⁄s و نتایج را با حل تحلیلی مقایسه کنید.
(ϕ-ϕ_0)/(∅_L-∅_0 )=(exp(ρux⁄Γ)-1)/(exp(ρuL/Γ)-1)
ج-محاسبه را برای u=2.5 m⁄s با 20 گره شبکه مجددا تکرار کنید و نتایج را با حل تحلیلی مقایسه کنید.
داده های زیر را اعمال کنید:
L=1 m, ρ=1 kg⁄m^(3,) , Γ=0.1 kg⁄ms

شکل(25)

اختلاف مرکزی(C-D):
الگوریتم حل:
1-وارد کردن اطلاعات اولیه شامل :
Γ, ρ,u,L,ϕ_A,ϕ_B,D_A,D_B
2- شبکه بندی: وارد نمودن تعداد گره ها و موقعیت هر گره نسبت به مبدا.
3- وارد کردن معادله ی حل دقیق و محاسبه ی آن در تمام نقاط.
4-محاسبه ی ضرایب a_W,a_E,a_P و S_u , S_P در هرگره با استفاده از روابط (30).
=لازم به توضیح است، پارامتر های فوق در گره های مرزی با وارد کردن شاخصه های مرزی در معادله ی (28) بصورت زیر می باشد:
a_P ϕ_P=a_W ϕ_W+a_E ϕ_E+S_u
a_P= a_W+a_E+( F_e-F_w )-S_P
a_W=0 , a_E= D_e- F_e/2 , S_u=(D_A+F_A ) ϕ_A , S_P=-(D_A+F_e ) : غربی مرز گره
a_W=D_w+ F_w/2 , 〖 a〗_E= 0 , S_u=(D_B+F_B ) ϕ_B , S_P=-(D_B-F_w ) : شرقی مرز گره

5- جایگذاری ضرایب فوق در رابطه ی (6) و تشکیل دستگاه معادلات جبری و حل آن به منظور یافتن دما در تمامی نقاط گرهی.
6- وارد کردن شرایط مرزی (وارد کردن دما در مرزها)
7-رسم نمودارهای حل عددی و تحلیلی و مقایسه ی آنها با یکدیگر.

متن برنامه
% % importing problem information
gamma=0.1;
ro=1;
L=1;
u=input('velocity:');

% % reticulation
n=input('nodes:');
nn=n+2;
dx=L/n;
F=ro*u; D=gamma/dx;
DA=2*gamma/dx; DB=DA;
fiA=1; fiB=0;
Fw=F; Fe=F;
x(1)=0; x(2)=dx/2; x(nn)=L;
for j=3:nn-1
x(j)=x(j-1)+dx;
end

% % exact solution
EXsolution=(fiB-fiA).*(exp(((ro*u).*x)./gamma)-1)./(exp((ro*u*L)/gamma)-1)+fiA;

% % calculation of cofficients in each node:a_W,a_E,a_P,S_u , S_P
for i=1:n
if i==1
SP(i)=-(DA+F); Su(i)=(DA+F)*fiA;
aW(i)=0; aE(i)=D-F/2;
elseif i==n
aW(i)=D+F/2; aE(i)=0;
SP(i)=-(DB-F); Su(i)=(DB-F)*fiB;
else
SP(i)=0; Su(i)=0;
aW(i)=D+F/2; aE(i)=D-F/2;
end
aP(i)=aW(i)+aE(i)+(Fe-Fw)-SP(i);
end

% % Formation of Algebraic equtions and solving it
aW(1)=[]; aW(n)=0;
aE=flipud(aE');
a=zeros(n,3);
a(:,1)=-aW;
a(:,2)=aP;
a(:,3)=-aE;
K = full(spdiags(a, [-1 0 1], n, n));
fi=Su/K';

% % Importing boundary condition
FI=ones(1,nn); FI(1)=fiA; FI(nn)=fiB;
FI(2:1:nn-1)=fi;

% % plotting results
plot(x,FI,'–rs','LineWidth',2)
hold on
plot(x,EXsolution,'k','LineWidth',2)
legend('Numerical solution','Exact solution')
xlabel('distance')
ylabel('phi')
title('A comparison of numerical and analaytical solution')
hold off;
خروجی
الف- u=0.1m⁄s

شکل(26)
ب-u=2.5 m⁄s

شکل(27)
ج-u=2.5 m⁄s و n=20

شکل(28)
نتیجه:
در حالت اول سرعت کم بود و شبکه درشت سازگاری خوبی داشت. حال آنکه در حالت دوم و سوم ملاحظه شد که با افزایش سرعت بایستی شبکه بندی نیز ریز تر شود. و برای شبکه بندی درشت حل عددی دارای نوسان بود. بنابراین می توان گفت در نظر گرفتن نسبت F/D در تحلیل نتایج و ارایه شبکه بندی بهتر ضروری است. هرچه این نسبت پایین باشد نتایج دقیق تر خواهد بود.

طرح اختلاف بالادست (Upwind)
الگوریتم حل:
1-وارد کردن اطلاعات اولیه شامل :
Γ, ρ,u,L,ϕ_A,ϕ_B,D_A,D_B
2- شبکه بندی: وارد نمودن تعداد گره ها و موقعیت هر گره نسبت به مبدا.
3- وارد کردن معادله ی حل دقیق و محاسبه ی آن در تمام نقاط.
4-محاسبه ی ضرایب a_W,a_E,a_P و S_u , S_P در هرگره با استفاده از روابط (30).
=لازم به توضیح است، پارامتر های فوق در گره های مرزی با وارد کردن شاخصه های مرزی در معادله ی (28) بصورت زیر می باشد:
a_P ϕ_P=a_W ϕ_W+a_E ϕ_E+S_u
a_P= a_W+a_E+( F_e-F_w )-S_P
a_W=0 , a_E= D_e , S_u=(D_A+F_A ) ϕ_A , S_P=-(D_A+F_e ) : غربی مرز گره
a_W=D_w+F_w , 〖 a〗_E= 0 , S_u=D_B ϕ_B , S_P=-D_B : شرقی مرز گره

متن برنامه
% % importing problem information
gamma=0.1;
ro=1;
L=1;
u=input('velocity:');

% % reticulation
n=input('nodes:');
nn=n+2;
dx=L/n;
F=ro*u; D=gamma/dx;
fiA=1; fiB=0;
DA=2*gamma/dx; DB=DA;
Fw=F; Fe=F;
x(1)=0; x(2)=dx/2; x(nn)=L;
for j=3:nn-1
x(j)=x(j-1)+dx;
end

% % exact solution
EXsolution=(fiB-fiA).*(exp(((ro*u).*x)./gamma)-1)./(exp((ro*u*L)/gamma)-1)+fiA;

% % calculation of cofficients in each node:a_W,a_E,a_P,S_u , S_P
for i=1:n
if F>0
aW(i)=D+F;
aE(i)=D;
else
aW(i)=D;
aE(i) =D-F;
end
if i==1
SP(i)=-(DA+F); Su(i)=(DA+F)*fiA;
aW(i)=0;
elseif i==n
aE(i)=0;
SP(i)=-(DB); Su(i)=(DB)*fiB;
else
SP(i)=0; Su(i)=0;
end
aP(i)=aW(i)+aE(i)+(Fe-Fw)-SP(i);
end

% % Formation of Algebraic equtions and solving it
aW(1)=[]; aW(n)=0;
aE=flipud(aE');
a=zeros(n,3);
a(:,1)=-aW;
a(:,2)=aP;
a(:,3)=-aE;
K = full(spdiags(a, [-1 0 1], n, n));
fi=Su/K';

% % Importing boundary condition
FI=ones(1,nn); FI(1)=fiA; FI(nn)=fiB;
FI(2:1:nn-1)=fi;

% % plotting results
plot(x,FI,'rs','LineWidth',2)
hold on
plot(x,EXsolution,'k','LineWidth',2)
legend('Numerical solution','Exact solution')
xlabel('distance')
ylabel('phi')
title('A comparison of numerical and analaytical solution')
axis([0 L 0 1.03])
hold off;
خروجی
الف- u=0.1m⁄s

شکل(29)
ب-u=2.5 m⁄s

شکل(30)
ج-u=2.5 m⁄s و n=20

شکل(30)
نتیجه:
طرح اختلاف مرکزی در تولید یک نتیجه مناسب و معقول، با همان اندازه شبکه موفق نبوده است. در طرح بالادست نتیجه واقع بینانه تری نسبت به اختلاف مرکزی ارایه داده شده است. به هر حال در نزدیکی مرز B خیلی به حل دقیق نزدیک نمی باشد.

طرح نمایی (Exponential)
الگوریتم حل:
-وارد کردن اطلاعات اولیه شامل :
Γ, ρ,u,L,ϕ_A,ϕ_B,D_A,D_B
2- شبکه بندی: وارد نمودن تعداد گره ها و موقعیت هر گره نسبت به مبدا.
3- وارد کردن معادله ی حل دقیق و محاسبه ی آن در تمام نقاط.
4-محاسبه ی ضرایب a_W,a_E,a_P و S_u , S_P در هرگره با استفاده از روابط (30).
=لازم به توضیح است، پارامتر های فوق در گره های مرزی با وارد کردن شاخصه های مرزی در معادله ی (28) بصورت زیر می باشد:
a_P ϕ_P=a_W ϕ_W+a_E ϕ_E+S_u
a_P= a_W+a_E+( F_e-F_w )-S_P
a_W=0 , a_E= a_E=F_e/(exp(F_e/D_e )-1) , S_u=F_A (1+1/(exp(F_w/(2D_w ))-1)) ϕ_A , S_P=〖-F〗_A (1+1/(exp(F_w/(2D_w ))-1)) : غربی مرز گره
a_W=(F_w exp(P_w ))/(exp(F_w/D_w )-1) , 〖 a〗_E= 0 , S_u=(F_B/(exp(F_e/D_e )-1)) ϕ_B , S_P=-(F_B/(exp(F_e/D_e )-1)) : شرقی مرز گره

متن برنامه
% % importing problem information
gamma=0.1;
ro=1;
L=1;
u=input('velocity:');

% % reticulation
n=input('nodes:');
nn=n+2;
dx=L/n;
F=ro*u; D=gamma/dx;
fiA=1; fiB=0;
Fw=F; Fe=F;
DA=2*gamma/dx; DB=DA;
P=F/D; PP=F/DA;
x(1)=0; x(2)=dx/2; x(nn)=L;
for j=3:nn-1
x(j)=x(j-1)+dx;
end

% % exact solution
EXsolution=(fiB-fiA).*(exp(((ro*u).*x)./gamma)-1)./(exp((ro*u*L)/gamma)-1)+fiA;

% % calculation of cofficients in each node:a_W,a_E,a_P,S_u , S_P
for i=1:n
if i==1
SP(i)=(-F*exp(PP))/(exp(PP)-1); Su(i)=((F*exp(PP))/(exp(PP)-1))*fiA;
aW(i)=0; aE(i)=F/(exp(P)-1);
elseif i==n
aW(i)=(F*exp(P))/(exp(P)-1); aE(i)=0;
SP(i)=(-F/(exp(PP)-1)); Su(i)=(F/(exp(PP)-1))*fiB;
else
SP(i)=0; Su(i)=0;
aW(i)=(F*exp(P))/(exp(P)-1); aE(i)=F/(exp(P)-1);
end
aP(i)=aW(i)+aE(i)-SP(i)+(Fe-Fw);

end

% % Formation of Algebraic equtions and solving it
aW(1)=[]; aW(n)=0;
aE=flipud(aE');
a=zeros(n,3);
a(:,1)=-aW;
a(:,2)=aP;
a(:,3)=-aE;
K = full(spdiags(a, [-1 0 1], n, n));
fi=Su/K';

% % Importing boundary condition
FI=ones(1,nn); FI(1)=fiA; FI(nn)=fiB;
FI(2:1:nn-1)=fi;

شکل(31)
ب-u=2.5 m⁄s

شکل(31)
ج-u=2.5 m⁄s و n=20

شکل(32)
نتیجه:
طرح نمایی نسبت به دو طرح قبلی نتیجه بسیار دقیق تری حتی در نزدیکی مرز B ارایه داد. به هر حال به نظر می رسد به دلیل محاسبه ی دشوار نماها این طرح کاربرد کمتری داشته باشد.
2

تعداد صفحات : 30 | فرمت فایل : WORD

بلافاصله بعد از پرداخت لینک دانلود فعال می شود