作者Ecampus (7.7)
看板Physics
標題Re: [疑問] 跑擴散方程式 到了終點..卻不會衰退
時間Mon Oct 16 00:47:10 2017
: http://imgur.com/TxeHaFC
:
: 這是一維的汙染物濃度擴散方程式
: u是速度 c是濃度
: E是擴散係數(延散係數)
: P是衰減係數
:
:
: 我用台灣某條河川 做實地測驗、與跑程式模擬兩個結果
: 實地測驗得到的曲線
: 就是高斯分布那一種型態(從0→到濃度最高點→到終點濃度衰退為0)
:
:
: 跑程式的話
: 依照各參數輸入 得到的結果卻是(0→濃度最高點→到終點仍然濃度最高點 沒衰退衰減過)
:
: → wohtp: 有圖嗎？你的初始條件和邊界條件是什麼？ 10/10 10:36
: → wohtp: 終點是時間上的終點嗎？多遠？ 10/10 10:39
您好 初始條件 我是假設開始的釋放點(0公尺處) 濃度為100ppb
邊界條件 我假設終點站是N 又C(N-1)=C(N+1) 所以濃度會一直無法下降
後來我改成C(N)=0
(距離是 0M~1600M)
(時間設定是0分鐘~160分鐘)
但是 也是到後面快終點的時候 濃度才會降下來
我想得到的圖形是類似高斯分布這樣的
https://imgur.com/w11jVaI
但我實際得到的結果是
https://imgur.com/vZSSWeh
=.=
我想做的是 時變的系統 濃度隨時間變化的
不過現在卡關 我知道時不變系統怎麼做....不知道時變系統怎麼做
求高手幫忙解惑 薄酬2000P
CC(I)=CC(I)+E*DT*(CC(I+1)-2*CC(I)+CC(I-1))/DX/DX-&
U*DT*(CC(I+1)-CC(I-1))/2./DX-&
DE*DT*CC(I)
CC(N)=0
我的時不變系統是這樣寫
如果改成時變..我還在想要如何改
: → Vulpix: 聽起來除了沒衰退有點詭異外，一段時間後濃度最高點位置不 10/10 20:20
: → Vulpix: 變，這件事不是什麼大問題啊。當然視乎你的IC、BC、t。 10/10 20:22
: → wohtp: 我本來想說是邊界條件沒設對，擴散到最旁邊又被擋回來之類 10/11 10:21
: → wohtp: 但是發現他還有用手放進去衰減 10/11 10:22
: → wohtp: 最後沒有變成零就是數值做錯了，沒話好說 10/11 10:23
: → saltlake: 確定參數輸入正確? 跑的程式誰寫的? 10/11 11:17
--
※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 125.230.93.84
※ 文章網址: https://webptt.com/m.aspx?n=bbs/Physics/M.1508086032.A.2D8.html
※ 編輯: Ecampus (125.230.93.84), 10/16/2017 00:57:05
1F：推 wohtp: 這你的方程式有解析解啦 10/16 02:38
2F：→ wohtp: 你要不要乾脆放棄數值算了 10/16 02:39
3F：→ Ecampus: 解析解我有 只是我要做的是數值解跟解析解的對比 10/16 03:10
4F：→ wohtp: 我覺得你首先連方程式是什麼意思好像都不清楚啊... 10/16 10:50
5F：→ wohtp: 能夠魔改成沒有時間就已經很神奇了 10/16 10:51
6F：→ wohtp: 然後你覺得拿掉時間以後得到的解是什麼意思？ 10/16 10:52
7F：→ wohtp: 我是不覺得你魔改過的方程式有什麼意義啦 10/16 10:55
8F：→ wohtp: 然後你的邊界條件又是什麼意思？ 10/16 10:56
9F：→ wohtp: 時變就把 dc/dt --> c(t+1)-c(t) = ... 10/16 11:40
10F：→ wohtp: 右邊那一堆都只跟 c(t) 有關 10/16 11:40
11F：→ wohtp: 給定 c(t) 你可以算 c(t+1) - c(t)，就可以找到 c(t+1) 10/16 11:41
12F：→ wohtp: 不過你最好想清楚你的邊界條件要什麼 10/16 11:42
13F：→ wohtp: 把邊邊直接釘在零實在很可疑 10/16 11:42
14F：→ Ecampus: 我試一下大大說的時變寫法看看 唉 當初沒好好學基本程式= 10/16 17:10
15F：→ Ecampus: =.=導致現在的苦果 10/16 17:10
17F：推 bluecadence: 用python寫了個簡單的模擬程式 畫出來的圖長這樣 10/18 03:54
19F：→ bluecadence: 橫軸是空間位置，縱軸是濃度。黑色曲線是起始濃度分 10/18 03:55
20F：→ bluecadence: 部(假設是高斯分佈)，然後其他紅色漸漸變淡的曲線是 10/18 03:57
21F：→ bluecadence: 隨時間改變的濃度空間分佈變化。 10/18 03:57
24F：推 sputtering: 樓上Good Job 10/18 07:55
25F：→ wohtp: 問一下b大邊界條件是什麼？ 10/18 10:26
26F：→ bluecadence: 回w大 還沒有放任何邊界條件 10/18 11:40
27F：→ bluecadence: 用的是forwrd Euler scheme 10/18 12:02
28F：→ bluecadence: forward 10/18 12:03
30F：→ wohtp: 不可能不放邊界條件吧？不然邊界附近的空間導數怎麼定義？ 10/18 14:37
31F：→ Vulpix: 大概可以當作Neumann BC吧？ 10/18 15:30
32F：→ Vulpix: 看那幾張圖，這個條件感覺比較符合code。 10/18 15:31
33F：→ Vulpix: 不過邊界的空間導數也可以用one-side的寫法。 10/18 15:39
34F：→ wohtp: 濃度到了邊界可不可以漏出去／怎麼漏出去可是很重要的 10/18 15:57
35F：→ wohtp: 要是隨便放個periodic bc或是會反射的，那整體衰減其實就跟 10/18 15:58
36F：→ wohtp: 擴散完全無關，不是嗎？ 10/18 15:59
37F：→ wohtp: 我想問的是這個 10/18 15:59
38F：推 bluecadence: 感謝w大的指點。昨天半夜睡不著無聊爬起來亂寫一通:p 10/18 17:19
39F：→ bluecadence: 我的確沒考量到邊界條件，也沒去處理這問題 10/18 17:20
40F：→ bluecadence: 我剛剛仔細重新研究一下 我用的方法是Crank-Nicolson 10/18 17:22
41F：→ bluecadence: implicit method 。過程中我沒有指定邊界值，但因為 10/18 17:25
42F：→ bluecadence: 解聯立方程式的矩陣的第一列和最後一列是有初始給定 10/18 17:27
43F：→ bluecadence: 特定值的,所以依舊算出來空間邊界的濃度，只是這個 10/18 17:28
44F：→ bluecadence: 算出來的值有沒有意義，的確有很大的問題。 10/18 17:30
45F：→ bluecadence: 第一次試著用數值方法解PDE，感謝指點 10/18 17:32
46F：→ bluecadence: 濃度似乎沒有"漏"出去 因為面積的積分隨時間變化，一 10/18 17:35
47F：→ bluecadence: 值保持相同 (面積~1) 10/18 17:36
49F：→ bluecadence: 看來邊界條件也很費心思。我想問在一個開放的場域， 10/18 19:06
50F：→ bluecadence: 原則上無限遠處濃度才等於0，那該如何設定合理的邊界 10/18 19:08
51F：→ bluecadence: 條件 10/18 19:08
52F：→ wohtp: 我也沒做過啊。會做的人也不必問這問題了？ 10/18 20:37
55F：→ bluecadence: 研究了一下，在數值方法中，邊界條件的給定 1. 當然 10/19 12:15
56F：→ bluecadence: 如果有實際物理系統的量值測最好 2. 根據物理系統給 10/19 12:16
57F：→ bluecadence: 與合理假設值 3. 假設線性，用"格子"內的點做外插 10/19 12:17
58F：→ bluecadence: 4. 用無限遠處的邊界值，回推估計。每種方法都必須檢 10/19 12:19
59F：→ bluecadence: 討誤差以及有效範圍。要得到接近現實的模擬，得費些 10/19 12:20
60F：→ bluecadence: 功夫就是。 10/19 12:21
61F：→ wohtp: 我想到的： 10/19 13:16
62F：→ wohtp: 1. 系統取得夠大，邊界上濃度一直不高 --> 直接釘在零 10/19 13:17
63F：→ wohtp: 2. 假設初始濃度都在原點附近，邊界上直接用 exp(-a x^2) 10/19 13:20
64F：→ wohtp: 的形式去 fit，然後把該有的衰減用手插進去 10/19 13:21
65F：→ sputtering: 可以問一下原PO那兩條方程的資料來源為何 10/19 15:20
66F：→ Ecampus: ??第1張綠色背景的圖 是1D濃度擴散方程式的公式Y 10/19 15:57
67F：→ Ecampus: 高斯分布那張圖 是TAYLOR G.I 在1953年提出的瞬時點源解 10/19 15:58
68F：→ Ecampus: TAYLOR的瞬時點源解 被稱為解析解 10/19 15:58
70F：→ sputtering: 又上了一課 10/19 18:36
不好意思 我的表達能力太差 讓大家看不懂我想問的東西...XD
其實我想問的 是下面這問題：
===============================
我們設定的最初濃度
在不同時間的時候 濃度應該也會隨時間而衰減吧~~~
(譬如 我在0m處的起點 丟了一把濃度10的粉末
過了1分鐘以後 0m處的濃度自然也會有衰減...
所以1分鐘後的0m處之濃度 也不會是10了)
請問這個問題 大家是如何處理的
https://imgur.com/TxeHaFC
譬如這個一維濃度擴散方程式而言
是用迭代處理嗎
顆顆 ..我打算直接用各種不同的T 手動帶進去
※ 編輯: Ecampus (125.230.93.173), 10/19/2017 19:08:02
71F：→ bluecadence: ㄜ..我完全誤解你的問題 你問的問題是 Taylor 10/19 20:18
72F：→ bluecadence: dispersion.. 所以請忽略我的圖 :p 10/19 20:18
73F：→ bluecadence: 你還多了一項 shear flow 的影響 10/19 20:25
74F：→ bluecadence: 還有多一項 external source 10/19 20:43
75F：→ bluecadence: 應該是把我給你的程式裡面的A和B矩陣的元素重新寫過 10/19 21:43
76F：→ bluecadence: (把shear flow 和external source考慮進去) 就可以了 10/19 21:43
77F：→ bluecadence: 但這要花一些時間推導公式。我考慮的是解整個PDE， 10/19 21:46
78F：→ bluecadence: 而不光是只想知道擴散物質的濃度在原點隨時間的變化 10/19 21:48
79F：→ bluecadence: 你稱作衰減係數的P 就是我說的 external source or 10/19 21:52
80F：→ bluecadence: sink 10/19 21:52
81F：→ wohtp: 把 d/dx 代換成 (d/dx - E/2u) 總之配完整平方搞掉一階項 10/19 23:11
82F：→ wohtp: 然後本來的解照樣用 10/19 23:12
83F：→ wohtp: 應該這樣就可以了吧 10/19 23:12
84F：→ Ecampus: !!!!!!!是的 立馬TRY一下 10/19 23:16
86F：→ bluecadence: 也就是你的 u > 0 的時候 10/20 01:17
87F：→ bluecadence: 在50的位置撒入一個gaussain分佈的擴散物，"河流"往 10/20 01:18
88F：→ bluecadence: 橫軸的右方流動 10/20 01:19
89F：→ bluecadence: (邊界條件必須更嚴格討論 不過大致上就這樣了) 10/20 01:20
90F：→ bluecadence: 要取空間任一處濃度隨時間變化 都沒問題的。 10/20 01:23
92F：→ bluecadence: 要扔 "delta function" 顯然在這數值計算會罩成困擾 10/20 01:55
93F：→ bluecadence: 圖形出來不合理 我就不放了 10/20 01:56
94F：→ bluecadence: 當把模擬空間加到更大例如1000，基本上不管是高斯或 10/20 01:58
95F：→ bluecadence: 是方形分佈，都能慢慢逼近"delta function" 10/20 01:59
96F：→ bluecadence: 這個東西玩二維應該會更好玩，diffusion 和 "河流"方 10/20 02:00
97F：→ bluecadence: 向的交互影響，還有流速快慢和diffusion constant的 10/20 02:01
98F：→ bluecadence: 大小變化 10/20 02:02
99F：→ bluecadence: 我這模擬狀況故意把河流流速放的非常慢，只是要展示 10/20 02:07
100F：→ bluecadence: 當 shear flow 出現時，對原本被動擴散的影響 10/20 02:08
101F：→ wohtp: 下游邊界的濃度似乎累積得蠻高的，所以可能不太準？ 10/20 15:41
102F：→ bluecadence: 純粹邊界條件問題 我上游設定0 下游用外插法 所以我 10/20 16:25
103F：→ bluecadence: 說邊界條件要更嚴格討論 不過我玩夠了 10/20 16:26
104F：→ bluecadence: 而且如果仔細看 下游累積到某個高點 已經開始下降 10/20 16:28
105F：→ bluecadence: 累積量和河流流速以及模擬時間長度會有關係 我肉眼 10/20 16:29
106F：→ bluecadence: 看不出來 暫時也不會想要嚴格探討誤差 好玩罷了 10/20 16:30
107F：→ bluecadence: 要讓邊界不累積很容易設定邊界0就好 但是會比較有意 10/20 16:32
108F：→ bluecadence: 義嗎 這我是不知道 10/20 16:32
109F：→ bluecadence: 我是沒把所有參數/條件全都清楚po上來 所以好像也沒 10/20 16:39
110F：→ bluecadence: 甚麼價值討論。原 po 如果想做類似的模擬，就把我 10/20 16:43
111F：→ bluecadence: 給的程式自己改一下矩陣A和B的元素，把 shear flow 10/20 16:44
112F：→ bluecadence: 項影響後造成的改變放進去就好了。source/sink 在 10/20 16:46
113F：→ bluecadence: 原po目前的狀況，大概還不必考慮 10/20 16:46
115F：→ bluecadence: 當初只是想玩玩不同的邊界處理方式 沒想太多 10/20 17:00
116F：→ Ecampus: 感謝大家 我來研究一下 感恩Q_Q 10/20 17:47
117F：→ Ecampus: 感謝大家 我來研究一下 感恩Q_Q 10/20 17:48