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

معادله ریچاردز برای شبیه سازی حرکت آب در محیط متخلخل تحت شرایط غیراشباع استفاده می‌شود و رطوبت خاک را به صورت تابعی از مکش و هدایت هیدرولیکی غیراشباع در عمق و زمان‌های مختلف شبیه‌سازی می‌نماید. معادله مذکور به دو صورت مکش (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')