作者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/cn.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