林翊展--養豬廢水處理設備的自動化設計--完整版

前言
前一陣子曾經協助老師進行一項雲林縣政府的計畫,該計畫內容為替雲林縣政府設計一座養豬廢水處理的示範場.在設計過程當中,有需要計算的地方,雖然用到的計算方法不難,但是希望利用matlab程式,只需要輸入所需的參數,即可以自動設計出符合標準的處理設備規格.

背景說明

養豬廢水處理

台灣的畜牧業以養豬為大宗,而且多是採企業化經營,規模大又集中,導致畜牧場所產生的廢水量相當驚人.目前台灣的總養豬頭數大約700萬頭,大部分都集中在台灣中南部,而一頭豬的糞尿廢水汙染力相當於人的六倍,可以折算成約4200萬人口.在過去,養豬廢水都是以直接排放的方式處理,對於台灣河川的破壞相當大,產生了許多”黑水河”.現在依環保署的規定,養豬規模超過20頭要強制建設廢水處理設備,而處理過的廢水也有品質的要求,如此一來,大大減少台灣河川所遭受的污染,而且在處理豬糞尿的過程中,還會產生許多具有附加價值的產物,如沼氣、有機堆肥等,達到保護環境與廢棄物再利用的雙重功效.

目前國內豬糞尿水的處理一般都採用三段式處理,過程主要包括固液分離、厭氣醱酵及好氣處理等步驟,其中再加入調整池、計量槽、初沉池、污泥回流及污泥處理等構成完整的養豬廢水處理系統,其流程如下圖:




而實際處理圖片則如下所示

1.固液分離




2.厭氣發酵




3.好氣發酵




經過三段式處理後的廢水即放流到河川之中

放流廢水




目的

程式設計的目標為輸入養豬頭數後,自動計算出廢水處理槽的規格參數,節省人力計算所需的時間.


名詞解釋
1.BOD：biochemical oxygen demand,生化需氧量.水中的雜質經由好氧菌分解所需要的氧氣量.由此數值可以看出廢水的生物可分解性,亦是廢水污染力的重要指標之一.BOD5則是指讓廢水暴露在氧氣充足的情況下,五天所耗費的氧氣量.

2.COD：chemical oxygen demand,化學需氧量.利用化學方法(強氧化劑),讓水中雜質氧化的耗氧量.廢水污染指標之一,且COD值必大於BOD值,而COD值通常被當作終極BOD值來作計算.

3.HRT:水力停留時間
HRT= reactor volume / effluent volume flux

4.BOD loading:反應槽對於進流廢水的負載,可以用BOD或COD表示,通常為經驗值,過多或過少都會對反應槽的效率產生影響.


設計原則
1.假設一頭肉豬100 kg,產生30L/day的廢水量
2.根據ASAE在1994公佈的STANDARDS 裡,每1000kg豬產生的BOD5約 3.1 kg/day.
3.反應槽的BOD loading假設為0.5 kg/m^3/day


程式內容

k=menu('請選擇計算方式','依廢水體積計算','依BOD負載計算');
pig=input('請輸入豬隻頭數 : ');
HRT=input('請輸入水力停留時間(hour) : ');
if k==1
    q=pig*30*0.001;
    V=q*(HRT/24);
    disp(['反應槽體積= ',num2str(V),' 立方公尺']);
elseif k==2
    V=pig*0.31/0.5*(HRT/24);
    disp(['反應槽體積= ',num2str(V),' 立方公尺']);
else
    break
end

執行結果
程式計算方式有兩種,第一是根據廢水量來計算,第二則是根據預估的反應槽loadind來做計算.兩種結果算出來後,取其較大值作為依據.

選擇圖示



根據廢水量計算



根據反應槽的loading來計算




根據以上的結果,養豬1000頭,水力停留時間6 hours的廢水處理槽,體積需要155立方公尺.


結果討論
1.兩種計算方式的數據會有不同,因為第一個根據廢水量的計算,只單就以物理性而言,槽體的大小是為了不讓水滿出來,但是第二項則是實驗的經驗值,參雜了細菌生長的好壞等生物性因素,所以兩者的參數不相同,必須取其大者.

2.水力停留時間6~8小時,是傳統反應槽的設計,如果有需要,例如loading改變,細菌生長狀況差異等,則會另行調整,以便達到處理效率.

3. 現行環保署畜牧業放流水標準
依據水污染防制法第七條第二項之規定,國內畜牧業適用之放流水標準,其水質項目及限制如下：
(1)非草食性動物(豬)：生化需氧量 < 80 mg/L；化學需氧量 < 600 mg/L；懸浮固體 < 150 mg/L。
(2)草食性動物(牛)：生化需氧量 < 80 mg/L；化學需氧量 < 450 mg/L；懸浮固體 < 150 mg/L。
(行政院環保署，2003)

4.以現實情況來說,廢水處理達不到標準通常受限於成本與土地大小,造成反應槽的大小不足,再加上生物性因素會比較難以掌握.

5.固液分離後的豬糞可以進行發酵做為有機堆肥,而厭氧反應過程中會產生沼氣提供發電,不過要先想辦法去除沼氣中的硫化氫,因為它會產生硫酸腐蝕發電機組.

6.此篇的程式與其他同學比較起來相對的簡單很多,因為生物反應槽大部分都是實驗參數,設計上比較難的地方在於參數的取用,用到計算上都只是簡單的計算,而現實上則隨時根據生物狀況來做調整,不過此篇可以提供同學們一點關於廢水處理的相關資訊.

林翊展--養豬廢水處理設備的自動化設計 v1

前言

前一陣子曾經協助老師進行一項雲林縣政府的計畫,該計畫內容為替雲林縣政府設計一座養豬廢水處理的示範場.在設計過程當中,有需要計算的地方,雖然用到的計算方法不難,但是希望利用matlab程式,只需要輸入所需的參數,即可以自動設計出符合標準的處理設備規格.

背景說明

養豬廢水處理

台灣的畜牧業以養豬為大宗，而且多是採企業化經營，規模大又集中，導致畜牧場所產生的廢水量相當驚人。目前台灣的總養豬頭數大約700萬頭，大部分都集中在台灣中南部，而一頭豬的糞尿廢水汙染力相當於人的六倍，可以折算成約4200萬人口。在過去，養豬廢水都是以直接排放的方式處理，對於台灣河川的破壞相當大，產生了許多”黑水河”。現在依環保署的規定，養豬規模超過20頭要強制建設廢水處理設備，而處理過的廢水也有品質的要求，如此一來，大大減少台灣河川所遭受的污染，而且在處理豬糞尿的過程中，還會產生許多具有附加價值的產物，如沼氣、有機堆肥等，達到保護環境與廢棄物再利用的雙重功效。

目前國內豬糞尿水的處理一般都採用三段式處理，過程主要包括固液分離、厭氣醱酵及好氣處理等步驟，其中再加入調整池、計量槽、初沉池、污泥回流及污泥處理等構成完整的養豬廢水處理系統，其流程如下圖:




而實際處理圖片則如下所示

1.固液分離




2.厭氣發酵




3.好氣發酵




經過三段式處理後的廢水即放流到河川之中

放流廢水




目的

程式設計的目標為輸入原廢水的化學性質，及其他必需化學參數後，自動計算出廢水處理槽的規格參數，節省人力計算所需的時間。

林翊展:廢水處理生化需氧量曲線

廢水處理工程裡,生化需氧量(BOD)是一項重要的參數,可以了解該廢水對於環境影響的程度.在BOD試驗當中,BOD值與反應時間為一條對數曲線,y=x*(1-exp(-k*t)),y為BOD值,x為終極BOD值,k為反應常數,t為時間.

u=input('請輸入終極BOD值 : ');
d=input('請輸入五天BOD值 : ');
k=-1/5*log(1-d/u);
t=0:30;
y=u*(1-exp(-k*t));
plot(t,y);
xlabel('reaction time (day)');
ylabel('BOD value (ppm)');
grid on;

程式會先要求輸入五天BOD值及終極BOD值,便可以畫出反應時間與BOD的曲線圖.
下圖為終極BOD=150,五天BOD=123的曲線圖

林翊展--水分蒸發潛熱

在濕氣特性的章節中,根據老師所提供的公式(Brooker, 1967)

1.當 273.16K ≦ 亁球溫度(Tdb) ≦338.72K 時

hfg = 2,502,535.259 – 2,385.76424 (Tdb- 273.16)

2.當 338.72K ≦ 亁球溫度(Tdb) ≦533.16K 時

hfg = (7,329,155,978,000 – 15,995,964.08 Tdb2)^1/2

3.當 255.38K ≦ 亁球溫度(Tdb) ≦273.16K 時

hsg =-2,839,683.144 – 212.56384(Tdb-255.38)


但是根據老師所給的程式,似乎無法跑出結果,所以另行撰寫了一個程式,可以直接輸入攝氏溫度,得到該溫度下相對應的水分蒸發潛熱.

程式碼如下:

####################################################

function [h]=hfg(Tc)
% function [h]=hfg(Tc)
% find the latent heat, kJ/kg, Tc dry-bulb temp in C, in column matrix
T=Tc+273.15;
if T>=273.16 & T<338.72
   h=2502535.259-2385.76424.*(T-273.16);
elseif T>=338.72 & T<533.16            
   h=(7329155978000-15995964.08.*T.^2).^1/2;       
elseif T>=255.38 & T<273.16
   h=-2839683.144-212.56384.*(T-255.38);
else
   disp('overrange')
end

####################################################

此functuion命名為"hfg"
執行結果如下:

指令執行結果一

如果以矩陣輸入,結果如下

指令執行結果二

再加入兩條指令,畫成曲線圖(攝氏1~10度)

>>plot(hfg([1:10]'));
>>xlabel('Temperatures, C');
>>ylabel('Evaporation heat, J/kg');

曲線圖

Problem #02 林翊展

第一題:
題旨:將量測的華式溫度轉換成攝氏後,找出體溫大於38的人,以矩陣D表示.當大於38顯示為1,小於則為0.
分析:將B矩陣中所有的元素進行轉換後,與38比較大小顯示於矩陣D.再利用尋找非零元素的方法,顯示其位置而容易尋找.
程式流程與分析:

>> clear
>> B=[98 100 102 104 98.4 98.2 98.5 101 102 99.5];
>> C=(B-32).*5/9; %華式溫標轉換成攝氏溫標的轉換公式
>> D=C>38

D =

0 0 1 1 0 0 0 1 1 0

>> [E]=find(D) %尋找非零元素並顯示其位址

E =

3 4 8 9

討論:雖然題目的要求只有顯示D矩陣,但是如果當被測者數量眾多(資料量龐大),海關人員難以尋找到底是哪一位的體溫超過38度,所以後面加了一個command,使得最後顯示為超過38度的人的序號,如果是二維矩陣則為位址.


第二題:
題旨:矩陣A與B,利用MATLAB程式進行比較,將結果顯示於D.
分析:當進行矩陣的運算時,要先注意其維數是否符合進行運算的條件.如果為矩陣與常數運算,則無此顧慮,但運算元前面要加"."
程式流程與分析:

>> clear
>> A=[-10 8 6 4 -5 20];B=[2 8 5 10 -6 3];
>> D=(A>B)

D =

0 0 1 0 1 1

>> D=(A>5)

D =

0 1 1 0 0 1 %比大小時,當符合則顯示1,不符合則顯示0

>> D=A+B

D =

-8 16 11 14 -11 23

討論:此題的兩個矩陣維數相同,所以能輕易進行運算.


第三題:
題旨:兩個矩陣x與y,進行加減乘除運算
分析:跟第二題一樣,矩陣進行運算須注意其維度,因為此為兩個列矩陣,無法進行矩陣乘除,只能讓相同位址的元素進行運算,此時運算元間要加".",常數則不用,而三角則直接進行元素的運算.

程式流程與分析:

>> clear
>> x=[10 20 30];y=[1 4 6];
>> A=3*x+y

A =

31 64 96

>> B=5*y./x

B =

0.5000 1.0000 1.0000

>> C=4*x.*(2*y)

C =

80 640 1440

>> D=sin(x).*cos(y)


D =

-0.2939 -0.5967 -0.9487 %三角函數以逕度量進行運算

>> E=5*x.*sin(2*y)

E =

45.4649 98.9358 -80.4859

討論:當沒有加上"."號時,程式會告知無法計算.三角函數則以逕度量進行運算,如果要以角度進行運算,指令為sind與cosd.


第四題:
題旨:建立M-file,先在檔案中將程式打好,然後在matlab主程式中呼叫即可.
分析:建立M-file的路徑"file->new->M-file",程式打好後存檔,在matlab主程式中,直接輸入M-file的檔名即可執行.

程式流程與分析:

1.M-file:
R=[10 30 200 400];
series=sum(R) %串聯的總電阻為全部相加
parallel=1/sum(1./R) %並聯的總電阻為全部的倒數相加的倒數

(檔名:resistor)

2.matlab主程式

>> clear
>> resistor

series =

640


parallel =

7.1006

討論:M-file可以很方便的用來debug,尤其當程式很大的時候.此程式的另一個方法可以使用函數檔.如果需要常常用到相同的計算式時,可以利用M-file製成函數檔,只須呼叫並將變數帶入即可.


第五題:
題旨:一個1~24的矩陣A,將它變成3x8,6x4,2x12的矩陣.
分析:利用resharp的指令來改變其行列.
程式流程與分析:

>> clear
>> A=1:24

A =

Columns 1 through 13

1 2 3 4 5 6 7 8 9 10 11 12 13

Columns 14 through 24

14 15 16 17 18 19 20 21 22 23 24

>> B=reshape(A,3,8)

B =

1 4 7 10 13 16 19 22
2 5 8 11 14 17 20 23
3 6 9 12 15 18 21 24

>> C=reshape(A,6,4)

C =

1 7 13 19
2 8 14 20
3 9 15 21
4 10 16 22
5 11 17 23
6 12 18 24

>> D=reshape(A,2,12)

D =

1 3 5 7 9 11 13 15 17 19 21 23
2 4 6 8 10 12 14 16 18 20 22 24

討論:建立一個元素為連續的矩陣時,可用":",可以省去麻煩.使用resharp指令時,要注意行數與列數的相乘是否符合,否則無法resharp.

Problem #01: 林翊展

題旨:matlab用於算數的基本操作

QUESTION 1:

>> 200*(20/60)/0.101972 %拉力(kg)*速度(m/s),換算成KW

ans =

653.7742

THE POWER IS 653.7742KW

QUESTION 2:

>> I=110/(1+20+200), %電流=電壓/電阻

I =

0.4977

>> P=110*I, %功率=電壓*電流

P =

54.7511

THE CURRENT IS 0.4977A , THE POWER IS 54.7511W

QUESTION 3:

>> m=650;n=428;o=282;s=(m+n+o)/2;
>> A=(s*(s-m)*(s-n)*(s-o))^(1/2) %海龍公式

A =

4.5233e+004

the area of triangle is 4.5233e+004 平方公分

QUESTION 4:

>> hog=(52000+84300-25000)/600

hog =

185.5000

it can raise 185 hogs.

QUESTION 5:

>> L=4;H=0.9;Q=3.33*[L-0.2*H]^(3/2)

Q =

24.8622

the flow rate is 24.8622 立方呎/秒

討論:
很容易忘記加乘號(*),matlab會以紅字提示格式錯誤的輸入要求,並標出錯誤點.

r95631028林翊展

學生姓名:林翊展
blog帳號:r95631028
e-mail:r95631028@ntu.edu.tw