ثبت نام در مدرسه تابستانه
روز
ساعت
دقیقه
ثانیه
سبد خرید

حل عددی معادله ریچادز در محیط متلب

حل عددی معادله ریچادز در محیط متلب

معادله ریچاردز برای شبیه سازی حرکت آب در محیط متخلخل تحت شرایط غیراشباع استفاده می‌شود و رطوبت خاک را به صورت تابعی از مکش و هدایت هیدرولیکی غیراشباع در عمق و زمان‌های مختلف شبیه‌سازی می‌نماید. معادله مذکور به دو صورت مکش (h-base) و رطوبت (Q-base) مورد استفاده قرار می‌گیرد. این معادله دارای مشتقات جزئی زمانی و مکانی می‌باشد و از این رو حل تحلیلی آن دشوار است و برای حل آن باید از روش‌های عددی بهره گرفت. یکی از روش‌های عددی تفاضل محدود (finite difference) می‌باشد. با استفاده از این روش می‌توان جملات دیفرانسیلی موجود در معادله باز و حل نمود. برای حل سریع و دقیق معادله نیاز به برنامه نویسی می‌باشد. در فایل ضمیمه معادله ریچاردز یک بعدی با استفاده از الگوریتم توماس در محیط متلب برنامه نویسی شده و با تغییر ورودی‌های مسئله می‌توان به شرایط دیگر بسط داد.

close all clear all clc % UNIT : cm & hr & % dz=input('dz='); L=input('L='); dt=input('dt='); T=input('T='); ks=input('ks='); Teta_r=input('Teta_r='); Teta_s=input('Teta_s='); n=input('n='); Alpha=input('Alpha='); tic; nn=round(T/dt); t1=round(0.25*nn); t2=round(0.5*nn); t3=round(0.75*nn); t4=nn; r=dt/(dz^2); z=0:dz:L; % initial Condition InitialTeta=0.2; N=L/dz+1; m=1-1/n; for i=1:N teta(i)=InitialTeta; Se=(teta(i)-Teta_r)/(Teta_s-Teta_r); h(i)=-1/Alpha*((1/Se)^(1/m)-1)^(1/n); hn(i)=h(i); end % Boundary Condition Ho=10; hn(1)=Ho; hn(N)=h(N); % Start calculations for t=1:nn for kk=1:1000; hnp=hn; % hydroulic conductivity and soil water retention calculate for i=2:N-1; hm=(hn(i)+h(i))/2; cm1=Alpha^n*(Teta_s-Teta_r)*m*n*abs(-hm)^(n-1); cm2=abs(1+abs(-Alpha*hm)^n)^(m+1); cm(i)=cm1/cm2; hm=(hn(i)+hn(i+1))/2; if hm>=0 km(i)=ks; else Se=1/(1+(-Alpha*hm)^n)^m; km(i)=ks*Se^.5*(1-(1-Se^(1/m))^m)^2; end hm=(hn(i)+hn(i-1))/2; if hm>=0 kn(i)=ks; else Se=1/(1+(-Alpha*hm)^n)^m; kn(i)=ks*(Se)^.5*(1-(1-Se^(1/m))^m)^2; end end % Coefficients matrix for i=2:N-2 c(i)=-r*km(i)/cm(i); end c(N-1)=0; for i=2:N-1 a(i)=-r*kn(i)/cm(i); end for i=2:N-1 b(i)=1-a(i)-c(i); end f(2)=h(2)-dz*(a(2)-c(2))-a(2)*hn(1); f(N-1)=h(N-1)-dz*(a(N-1)-c(N-1))-c(N-1)*hn(N); for i=3:N-2 f(i)=h(i)-dz*(a(i)-c(i)); end % Tomas algorithm for i=2:N-1 if i==2 alfa(i)=b(i); else alfa(i)=b(i)-a(i)*beta(i-1); end beta(i)=c(i)/alfa(i); end for i=2:N-1 if i==2 y(i)=f(i)/alfa(i); else y(i)=(f(i)-a(i)*y(i-1))/alfa(i); end end for i=N-1:-1:2 if i==N-1 hn(i)=y(i); else hn(i)=y(i)-beta(i)*hn(i+1); end end er=abs(hn-hnp); if er<0.0000001 break end end % Specify graph data if t==t1 hn1=hn; for i=1:N if hn1(i)>=0 teta1(i)=Teta_s; else teta1(i)=(Teta_s-Teta_r)/(1+abs(-1*hn1(i))^n)^m+Teta_r; end end end if t==t2 hn2=hn; for i=1:N if hn2(i)>=0 teta2(i)=Teta_s; else teta2(i)=(Teta_s-Teta_r)/(1+abs(-1*hn2(i))^n)^m+Teta_r; end end end if t==t3 hn3=hn; for i=1:N if hn3(i)>=0 teta3(i)=Teta_s; else teta3(i)=(Teta_s-Teta_r)/(1+abs(-1*hn3(i))^n)^m+Teta_r; end end end if t==t4 hn4=hn; for i=1:N if hn4(i)>=0 teta4(i)=Teta_s; else teta4(i)=(Teta_s-Teta_r)/(1+abs(-1*hn4(i))^n)^m+Teta_r; end end end h=hn; end z=z*-1; plot(teta4,z,teta3,z,teta2,z,teta1,z),xlabel('water content(%)'),ylabel('depth(cm)') xlim([0.2 0.45]) ylim([-100 0]) toc; % subplot(3,1,2),plot(teta4,z),xlabel('water Content'),ylabel('depth(cm)') % subplot(3,1,1),plot(hn4,z),xlabel('Pressure'),ylabel('depth(cm)') % subplot(3,1,3),plot(teta4,hn4),xlabel('water content'),ylabel('Pressure')

 

 

امتیاز: post
دیدگاه‌ها ۱