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

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

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

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