污水排放问题

这次出的题目是一个污水扩散的一个模型, 需要计算在下游3km处的一个工厂的排放污染量最大排多少才能使得在下游10km处的一个监测站监测的数据不超过10mg/L.

这次的模型相对复杂, 最后我们结果也不是太满意, 做了一个2维扩散的模型, 我是编程的, 对模型部分就不多做讨论了, 主要进行两个方面的讨论, 工厂是连续投放污染物还是一次性投放污染物, 对于一次投放的计算代码如下.

i = 1;
start = 1.5;
en = 3;
len = size(start:0.01:en,2);
temp = 0.005;
flag = 0;
resf = zeros(1,len);
res = zeros(1,len);
for ux = start:0.01:en
    M = 0.01;
    while flag < 10
        M = M + temp;
        time = 7000/ux;
        x = targetx(7000,time-1:0.01:time+1,ux,M);
        flag = max(x);
    end
    resf(i) = flag;
    flag = 0;
    i
    res(i) = M;
    i = i+1;
end
figure(1);
plot(start:0.01:en,res,'r');
axis tight;title('不同流速下最大一次排放量');
xlabel('流速/(m/s)'); ylabel('所能承受最大一次性排放量');
figure(2)
plot(start:0.01:en,resf-10,'r');
axis tight;title('所求最大排放量超出10mg/L的数量');
xlabel('流速/(m/s)'); ylabel('超出量/(mg/L)');
figure(3)
plot(start:0.01:en,7000./(start:0.01:en),'r');
residual = resf-10;
axis tight;title('所求最大排放量在7000米处到达污染峰值时间');
xlabel('流速/(m/s)'); ylabel('时间/t');
save('resultx.mat','res','residual')

由于是求一个不等式, 解析解难求, 所以是数值求解, 具体求法是对参数进行逐步递增, 直到超出不等式要求, 然后保存结果. 由于题目中没有给速度, 所以自己给了1.5m/s-3.0m/s的速度以0.01m为间隔求解了150组结果



对于连续的情况, 新增了一个参数, 工厂每秒钟排放的升数, 因此有两个维度.这里同样速度在1.5m/s-3.0m/s,不过以0.1m/s为间隔, 对于排放量, 在50L-1500L之间以50L为间隔去了30组排放量, 也就有了480组结果.代码如下

start = 1.5;
en = 3;
du = 0.1;
tstart = 50;
ten = 1500;
dc = 50;
sequ = [start:du:en];
len = size(sequ,2);
seqc = tstart:dc:ten;
tlen = size(seqc,2);
temp = 0.1;
flag = 0;
resf = zeros(len,tlen);
res = zeros(len,tlen);
resc = zeros(len,tlen);
i = 1;
j = 1;
global ux cap;
for ux = sequ
    time = 7000/ux-1;
    j = 1;
    for cap = seqc
        t = 0.001;
        flag = 0;
        seq = (time/3):1:(time*3);
        x = targetx(7000,seq,ux,cap);
        seqtime = seq(find(x ~= 0));
        starttime = seqtime(1);
        endtime = seqtime(end);
        length = endtime-starttime;
        while flag <10
            flag = quadl(@conti3D,starttime,starttime+t);
            t = t+0.01;
            if t >= length
                flag = quadl(@conti3D,starttime,endtime);
                break;
            end
        end
        resc(i,j) = cap;
        res(i,j) = t;
        resf(i,j) = flag;
        flag;
        j = j+1
    end
    i = i+1
end
time = res;
residual = resf-10;
totalcap = resc.*res;
[x,y] = meshgrid(seqc,sequ);
velocity = y;
vm = x./5;
vm = fliplr(vm);
figure(1);
totalcap = totalcap./5;
surf(vm,velocity,totalcap);colorbar;
axis tight;xlabel('工厂每秒排放量/(L/s)'); ylabel('流速/(m/s)');zlabel('最大总排放量/L');
title('不同流速与每秒排放量下最大总排放量');
figure(2)
surf(vm,velocity,residual);colorbar;
axis tight;xlabel('工厂每秒排放量/(L/s)'); ylabel('流速/(m/s)');zlabel('最大排放量超出浓度/(mg/L)');
title('所求最大排放量在7km处最高浓度超出10mg/L的数量');
result = zeros(480,5);
for i = 1:16
    for j = 1:30
        result(30*(i-1)+j,1) = vm(i,j);
        result(30*(i-1)+j,2) = velocity(i,j);
        result(30*(i-1)+j,3) = time(i,j);
        result(30*(i-1)+j,4) = residual(i,j);
        result(30*(i-1)+j,5) = totalcap(i,j);
    end
end
save('resultxconti','vm','velocity','time','residual','totalcap','result');

结果如下

最后结果还好吧, 不过今晚答辩没有去, 就在这里总结一点吧, 关于第二题的我还没有总结, 下次再写吧, 现在都12点50了, 先睡觉明天考科四了, 昨天晚上10点回去刷到了1点多把1000道题刷完了= =, 不过今天没有怎么做了, 准备明天早上在模考几次, 就差不多了, 驾照这件事总算完结了, 从上个暑假弄到现在也是醉醉的…

后话

搞个数模确实有点累, 今天晚上坐车回了家, 准备明天考科四, 正好下午把这三天的任务做完了, 然后5点大致和队友们把论文弄差不多了我就直接坐车去车站了, 结果5点半到了车站才发现身份证没有带, 当时我还想了一下没有身份证能不能取票, 然后才想起没带身份证有票也不能坐车啊, 马丹.

赶紧又叫了辆车回去拿了身份证又坐车回来, 还好车子晚点10分钟了7点22才到, 我是7点快10分才跑到检票口那, 差点没敢到= =而且火车票只要12.5块, 确实便宜…