حل معادله جریان غیرماندگار در اطراف لوله زهکش به صورت عددی با استفاده از روش ضمنی و صریح

حل معادله جریان غیرماندگار در اطراف لوله زهکش به صورت عددی با استفاده از روش ضمنی و صریح در محیط متلب

یکی از بحث‌های مهم در آب‌های زیرزمینی و مدل‌سازی زهکشی تعیین جهت جریان و بررسی نوسانات سطح ایستابی در اثر تخلیه و تغذیه می‌باشد. برای این منظور بایستی معادلات دیفرانسیلی حاکم بر جریان حل شود.

معادله یک بعدی غیرماندگار به صورت زیر تعریف می‌شود :

در این رابطه R مقدار تغذیه، S ضریب ذخیره و T ضریب انتقال (KD) می‌باشد.

معادله فوق با الگوی صریح و ضمنی قابل حل می‌باشد، که در فایل پیوست ارائه شده است. در مثال ارائه شده جریان غیر ماندگار در اطراف لوله زهکش به فاصله 30 متر بررسی و حل شده است. سطح ایستابی اولیه 3/3 متر و در 10 روز  به 2/8 متر می‌رسد. در صورت لزوم می‌توان داده‌های ورودی را تغییر داد و مسئله برای شرایط دیگر حل گردد. در صورت تغییر داده‌های مسئله رعایت طول گام‌های زمانی و مکانی مناسب و تعداد حلقه‌های به کار برده شده برای گرفتن جواب مناسب ضروری است.

در ادامه میتوانید کدهای مربوط به این بحث را دریافت نمایید:

clc
clear
dt=5;
dx=1.5;
T=0.0013;
s=0.07;
x=0:dx:15;
for i=1:31
    ho(i)=3.3;
end
a=0.29;
hn(1)=3.3;
hn(11)=2.8;
jt=0;
for t=1:14400
    jt=jt+dt;
    for i=2:10
        hn(i)=ho(i)+T*dt/(s*dx^2)*(ho(i+1)-2*ho(i)+ho(i-1));
    end
    ho=hn;
    hn(1)=hn(2);
    if jt==1440
        h1440=hn;
        q=500*s*(1-(0.8*exp(-a*1)));
    elseif jt==2880
        h2880=hn;
        q1=500*s*(1-(0.8*exp(-a*2)));
    elseif jt==4320
        h4320=hn;
        q2=500*s*(1-(0.8*exp(-a*3)));
    elseif jt==5760
        h5760=hn;
         q3=500*s*(1-(0.8*exp(-a*4)));
    elseif jt==7200
        h7200=hn;
         q4=500*s*(1-(0.8*exp(-a*5)));
     elseif jt==8640
        h8640=hn;
        q5=500*s*(1-(0.8*exp(-a*6)));
    elseif jt==10080
        h10080=hn;
        q6=500*s*(1-(0.8*exp(-a*7)));
    elseif jt==11520
        h11520=hn;
        q7=500*s*(1-(0.8*exp(-a*8)));
    elseif jt==12960
        h12960=hn;
        q8=500*s*(1-(0.8*exp(-a*9)));
        elseif jt==14400
        h14400=hn;
        q9=500*s*(1-(0.8*exp(-a*10)));
    end
end
hold on
plot(x,h1440)
plot(x,h2880)
plot(x,h4320)
plot(x,h5760)
plot(x,h7200)
plot(x,h8640)
plot(x,h10080)
plot(x,h11520)
plot(x,h12960)
plot(x,h14400)
hold off
Qsarih=[q;q1;q2;q3;q4;q5;q6;q7;q8;q9];
filename = 'discharge.xlsx';
sheet = 1;
xlRange = 'B1';
xlswrite(filename,Qsarih,sheet,xlRange)



clc
clear
dt=1;
dx=0.25;
T=0.0013;
s=0.07;
x=0:dx:12.5;
A=2+(s*(dx^2))/(T*dt)
for i=1:51
    ho(i)=3.3;
end
a=0.4176
hn(1)=3.3;
hn(51)=2.8;
jt=0;
for t=1:dt:14400
    jt=jt+dt;
    for i=2:50
       hn(i)=ho(i);
    end
    for k=1:100
    for i=2:50
        old=hn(i);
        hn(i)=(1/A)*((hn(i-1))+hn(i+1)+(A-2)*ho(i));
        er(i)=abs(hn(i)-old);
    end
    if er<0.00000000000001
        break
    end
    end
    ho=hn;
    hn(1)=hn(2);
       if jt==1440
        h1440=hn;
        q=500*s*(1-(0.8*exp(-a*1)))
    elseif jt==2880
        h2880=hn;
        q1=500*s*(1-(0.8*exp(-a*2)));
    elseif jt==4320
        h4320=hn;
        q2=500*s*(1-(0.8*exp(-a*3)))
    elseif jt==5760
        h5760=hn;
         q3=500*s*(1-(0.8*exp(-a*4)))
    elseif jt==7200
        h7200=hn;
         q4=500*s*(1-(0.8*exp(-a*5)))
     elseif jt==8640
        h8640=hn;
        q5=500*s*(1-(0.8*exp(-a*6)))
    elseif jt==10080
        h10080=hn;
        q6=500*s*(1-(0.8*exp(-a*7)))
    elseif jt==11520
        h11520=hn;
        q7=500*s*(1-(0.8*exp(-a*8)))
    elseif jt==12960
        h12960=hn;
        q8=500*s*(1-(0.8*exp(-a*9)))
        elseif jt==14400
        h14400=hn;
        q9=500*s*(1-(0.8*exp(-a*10)))
    end
end
hold on
plot(x,h1440)
plot(x,h2880)
plot(x,h4320)
plot(x,h5760)
plot(x,h7200)
plot(x,h8640)
plot(x,h10080)
plot(x,h11520)
plot(x,h12960)
plot(x,h14400)
hold off
Qzemni=[q;q1;q2;q3;q4;q5;q6;q7;q8;q9];
filename = 'discharge.xlsx';
sheet = 1;
xlRange = 'A1';
xlswrite(filename,Qzemni,sheet,xlRange)