第一章仿真 -3p13y-2(3ky)3+g (1.72) 动量守恒方程的边界条件为:进水口和出水口处流体垂直流速为0,在进水 冂或出水凵断面处流体水平流速各点相同。在沉淀池水面,流体垂直流速为0, 水平流速有对称性。在池壁与水的界面,与池壁接触的流体在与池壁相切方向 的流速和垂直方向的流速均为0,通过池壁的净质量通量为0。根据计算流体力 学,池壁对流体的剪切张力x0为 r6=C2kk△a/nEX (1.73) 式中,C是经验常数,约为0.09;k是湍流动能;κ是常数,对于平滑池壁为 0.4187;E也是常数,对于平滑池壁为9793;X是无量纲数,为11.63。由此 上式变成: ro=0.048m△ (1.74) 动量方程中压力项的含义是流场中某一点的压力与该点流体静压力之差 当流体静止不动时,压力项p处处为0;此时若流体密度恒定而浮力很小,则动 量方程中的重力项可以忽略不计。但若流体密度不是处处相同,且浮力项须考 虑的话,此时重力项须写成(P-p)8,其中p是流体参考密度,而压力项p则 须定义为流场中某一点的压力与该点流体参考静压力Pg△y之差。 5.湍流动能方程 根据流体力学,不可压缩流体湍流的平均动能方程为 a(映)t=-V·[uk]+V·1[p1+(pr1)·k+(p12)(△:△)- (1.75) 式中,是湍流动能,H1是层流粘度,r是湍流粘度流体粘度μ=1+pr,Hr 0.09k21e,e是湍流动能耗散速率。方程(1.75)左侧是流场中单位体积内流 体湍流动能变化的速率。右侧第一项是因对流引起的单位体积湍流动能的变 化;第二项是因扩散引起的单位体积流体湍流动能的变化,扩散系数与粘度有 关,σ是常数,一般为1.0;第三项是单位体积流体内湍流动能的产生项,Δ:Δ在 二维系统中等于4[(n,1x)2+(ov,/3y)2]+2(a,1x)+(bv,|ly)]2;第四 项是单位体积流体湍流动能的耗散项。 由于湍流动能是标量,二维空间湍流动能的守恒方程为 ,t+引(Pm,p一[1+(p1,)](1x)x+pk-【m1+(m7/) (31y)1y=pr2(au2|x)2+(au,|ay)2]+[(at,1x) +(au.lay)]:-pe (1.76) 若平均动能产生项及层流粘度项可忽略,则方程(1.75,1.76)分别变成 a(ak )/at=-v. ouk]+V[(urlok)Vk]-pe (1.77)
第一节模型的建立 27 7,|t+Dmk-(m1a)(ak1x)]1x+,k-(1.)(k(y)]1y= 湍流动能方程的边界条件是:在沉淀池入水口,湍流动能为: 式中,Ⅰ为平均湍流强度,数值为0.1-0.3;在沉淀池池壁处,湍流动能的产生 与耗散达到平衡,湍流动能为: 在沉淀池废水液面,ky=0;在沉淀池出水口,ak|ax=0。 6.湍流涡流耗散 根据流体力学,流场中的湍流动能耗散速率方程为 a(c)13t=-V·[pe]+V·[1.+(pr/a)v·e +c1(ek)(p/2)(△:△)-c2E2/k (1.81) 式中,;、C2是经验常数,为1.44,c2为1.92。方程(181)左侧是流场中单位 体积内流体湍流动能耗散的变化速率。右侧第一项是因对流引起的单位体积湍 流动能耗散的变化速率;第二项是因扩散引起的单位体积流体湍流动能耗散的 变化速率,扩散系数与粘度有关,a,是常数,一般为1.22;第三项是单位体积流 体内湍流动能的产生对湍流动能耗散速率的影响;第四项是因湍流涡流的解体 与耗散引起的单位体积流体湍流动能的耗散的变化速率。 出于湍流动能耗散速率是标量,二维空间湍流动能耗散速率方程为 o(xe)t+4=-[1+(pr1)](e1x)1x+Pm,-[m +(Hr/a.)](ae3y)13y=c1(e/k)yr2[(0u,x)2+(,1by)2] +[(bulx)+(1(y)]2}-c22 (1.82) 若忽略平均动能产生项及层流粘度项,则方程(181,1.82)分别变成 a(e)3t=-V·pwe]+V·[(Hrl)Ve]-c2E/k (1.83) a( pe)aL +alpu, e-(plo,(ae/ar)Jax+aou, e-(urAa, )(dE/ay)]/dy 力k 湍流动能耗散速率方程的边界条件是:在沉淀池人水口,湍流动能为 k=3(Iv2)2/2 (1.85) 式中,为平均湍流强度,数值为0.1~0.3;在沉淀池池壁处,湍流动能的产生 与耗散达到平衡,湍流动能为 k=stolid 在沉淀池水,ky=0;在沉淀池出水口,ak1x=0,经验常数C=0.09。 7.思浮固体质量守恒 若悬浮固体质量分数为s,则流体密度为:
第一章仿 p=p, /1-s[1-(p, p)iF (1.87) 式中pn是水的密度、即参考密度,,是悬浮固体的密度。悬浮固体的浓度为 1.88 悬浮固体质量守恒方程为: a(a5)/t=-V·[s]+V·[1+(pr{o、)]vs-V·[ps] (I.89) 方程(1.89)左侧为悬浮固体的浓度变化速率。右侧第一项是因对流引起的 单位体积流体内的固体质量变化速率;第二项是因涡流扩散引起的固体浓度变 化速率,σ为常数,般为1;第三项是因沉降引起的固体浓度变化速率,v是固 体相对于液体的沉降速率向量。 在二维空间,方程(1.89)可以写成 (ps )/at+ai pu S-LAL +(urlo,)(as)ax)1/ax+aipu, s-[uy +(r{a,)](s/3y)}ay=a(mu,5)y (1.90) 若忽略层流粘度项,方程(1.90)可变成: a(ps)/at+ aLou s-(u as(as/ax)1/ax +aLou, s-(urha(ds/ay)]/ay=a(ov,s)/ay (1.91) 固体悬浮颗粒的沉降速率与悬浮物浓度有关,计算沉降速率v,的方程为: p-rep (1.92) 式中,v。是最大理论沉降速率,r是拥挤沉降参数,rp是低浓度缓慢沉降组分 参数。P=P.-Pm=p,- Fns p,Pm是出水可达到的最小悬浮固体浓度,F 是进水悬浮固体中不可沉降部分的分数,pn是进水中固体的浓度。 悬浮固体质量守恒方程的边界条件是;在池壁和池水表面,固体净通量为 0。在池底,固体在接触池底后即假定被除去,此时形成固体沉降的净通量。 8.基本微分方程 比较流体的质量守恒方程(1.43)动量守恒方程(148)、湍流动能守恒方程 (1.75)湍流动能耗散速率方程(181)及悬浮固体质量守恒方程(1.89),可知 若用基本变量来代替上述方程中的相应状态变量,可以写出一个基本微分方 程来表示上述5种方程的性质。该基本微分方程为 a(p中)t=-V·[中]+V·[rV中+S 式中,中是单位质量流体中的基本量(标量或矢量),u是流体流速向量,ρ为流 体密度,t为时间。方程(1.93)中右侧第一项为对流项,对流通量为pk中;第 项为扩散项,扩散通量与梯度∮成正比,r为基本扩散系数,具体内容由卓确 定;S为其它项或称源项。 对于质量守恒方程(143),基本微分方程(1.93)在形式上与其完全一致。 对于动量守恒方程(1.71,1.72),基本微分方程(1.93)中的源项为方程
第二节模型的分析 29 (1.71)右侧的p/3x-2(0k|3x)/3,或(1.72)中的-0p/3y-2(3ly)3+ 对于湍流动能守恒方程(1.75),基本微分方程(1.93)中的源项为方程 (1.75)右侧的-c。 对于湍流动能耗散速率方程(1.81),基本微分方程(1.93)中的源项为方程 (1.81)右侧的-c22/k。 对于悬浮固体守恒方程(1.90),基本微分方程(1.93)中的源项为方程 (1.90)右侧的(v5)by。 方程(1.93)即为本例所求的描述沉淀池二维流场中流体流速分布和固体颗 粒浓度分布的基本数学模型。 第二节模型的分析 仿真的第二步,是要对模型进行分析。所谓对模型进行分析,就是在不同的 边界条件或参数设定下对模型求解,并从求得的解中获得所研究对象或过程的 动态性质。 具有集总参数特征的过程,其数学模型一般为微分方程(组),可用四阶龙格 库塔法求解;具有分布参数特征的过程,其数学模型一般为偏微分方程(组), 可用有限差分法求解。龙格一库塔法是在单一方向上(如时间)的差分法;有限 差分法是多个方向上(如时间和空间)的差分法。 四阶龙格一库塔法 若有微分方程dydt=f(t,y),已知t=tn时,y=yn,时问t迭代计算的步 长为h,则在t=tn时 y+1=y+(k1+2k2+2k3+k4)/6 式中,k=h·f(tn,yn) k2=h·f(t+h12,yn+k12) k3=h·f(,+h/2,yn+k2/2) k,=h·f(tn+h,yn+k3) 由以上方程可知,只要已知初始时刻的y值即可通过迭代计算,算出任意 时刻y的数值。计算原理可从图18(a)、(b)、(c)、(d)中看出。 在图1.8(a)中,根据方程dy/dt=f(t,y),将tn、y,值代入/(,y)即可得到 dydt值也即该方程曲线在tn、y,点的斜率。由图可见,斜率为f(ln,yn)的
第一章仿真 直线与y=yn、t=tn+h的直线相交形成三角形。三角形底边长为h,斜边斜 率为∫(tn,yn),因此三角形另一直角边的长度为h·f(tn,yn),即k1 在图1.8(b)中,根据方程dy/dt=f(t,y),将t+h2、y,+k12值代入 f(4,y)即可得到又个dydt值,也即该方程在tn+h12、yn+k1/2点的斜率 (注意:表示该斜率的虚线不与方程曲线相切)。通过tn、y,点作虚线的平行线, 与!=tn+h的直线相交后形成三角形。三角形底边长为h,斜边斜率为 f(t+h12,yn+k1/2),因此三角形另一直角边的长度为h·f(tn+h2,yn+k1 2),即k,。 图1.8(c)、(d)中的情况与1.8(b)中相似,读者可自行推导。 对于含多个变量的微分方程组,上述计算方法仍然有效。例如,若有一阶微 分方程组 /dt=g(t, x, y) 已知在tn时,x=xn,y=yn,时间步长为h,则在t=tn+h时 In+=I,+(k,+2k,+2k3+k4)/6 y=yn+(l1+2l2+2l3+l4)6 式中,k1=hf(tn,x,,yn) k2=hf(t,+h/2,r,+k/2, y,+1112) k3=hf(Lm+h/2,I+k2/2, y,+12/2) k4=hf(t,+ l3) t,r hg(tn+h2,x。+k1/2 l3=hg(tn + h/2, I,+k2/2, y+ 12/2) L=hg(t, +h,xn+k3, y+l3) 在用四阶龙格一库塔法求解微分方程组模型时,初值的选取分重要。原 因是,微分方程组的收敛点不一定是惟一的,而龙格一库塔法是数值方法,只要 本次迭代计算结果与上次计算结果之差小于某一个较小的数(以某种方式表 示),即承认本次迭代结果为方程组的解。因此,若初值选取不当,则有可能在计 算结束时,得到的解并不符合实际。 正确的初值选取在很大程度上要依靠对研究对象的深刻了解。为了避免差 错,可以进行多组不同初值的选取与试算,并比较各次计算结果。在计算完成 后,计算结果还必须与工程实际进行比照