[求助]能量方程的源项？不知道这样理解对不对？请看看

zlf81

    以下是在看书和程序的时候碰到的问题。
    能量方程为：
    d(ρCpΦ)/dt+d(ρCpuΦ)/dx+d(ρCpvΦ)/dy=k(d^2Φ/dx^2)+k(d^2Φ/dy^2)+S  （1）
    而simplec算法求解的是一个通用的控制方程，
    d(ρΦ)/dt+d(ρuΦ)/dx+d(ρvΦ)/dy=Г(d^2Φ/dx^2)+Г(d^2Φ/dy^2)+S   （2）
    要求解的话，需要将式（1）转化为式（2），这样的话，是不是
    方程（2）中扩散系数Г＝k/Cp？？
    方程（2）中源项S应该等于方程（2）中的源项S除以比热Cp，即S=S/Cp。？？
    不知道这样的理解对不对？？望指点！！thanks！！

anyone

解温度基本就这样吧.
也可以解焓或内能.
如果 Φ=CpT, S就不用除Cp.

zlf81

    你的意思是指：
    Φ=CpT，
    则，扩散系数Г还是＝k/Cp，式（1）中的源项S保持不变。
    利用simplec算出Φ以后，
    再除以Cp，
    就得到了温度T。
    是这样理解把？
    谢谢anyone了！

zlf81

    不对啊，我用两种方法计算的结果相差很大，不知道为什么？
    我是这样做的：
    第一种：在void user_gamsor()这样设置
    { 。。。
for(j=1;j<=m1;j++)
  for(i=1;i<=l1;i++)
  gam(i,j)=thercon(i,j)/spe(i,j);//即Г＝k/Cp
         for (j=2;j<=m2;j+＋）
            for (i=2;i<=l2;i++)
{  con[j]=con[j]/spe(i,j);
   ap[j]=ap[j]/spe(i,j);}
      。。。
     }
     然后在user_output()中输出温度t(i,j)
   
     第二种方法
    在void user_gamsor()这样设置
    { 。。。
for(j=1;j<=m1;j++)
  for(i=1;i<=l1;i++)
  gam(i,j)=thercon(i,j)/spe(i,j);//即Г＝k/Cp，不变
         for (j=2;j<=m2;j+＋）
            for (i=2;i<=l2;i++)
{  con[j]=con[j];  //不除以spe(i,j）
   ap[j]=ap[j];}
      。。。
     }
     在user_output()中这样设置
     {for(i = 2; i <= l2; i ++)
for(j =2; j <= m2; j ++)
    t(i,j) = t(i,j)/ spe(i,j);  //spe(i,j)是各点Cp值。
     }
     
     这两种方法算出的温度相差很大，第一种最高温度有1100度，第二种最高温度才300度。更实际结果相比，第二种相差太大，当然第一种结果也有问题。这其中我用了附加源项法处理了第二和三类边界条件，不知道这个有没有影响？？
     急问为什么啊？？哪错了？？

anyone

边界条件处理会有关系.
到底具体是什么问题? 可以画个简图看看.

zlf81

    二维的情况，是模拟焊接溶池的温度场。
    具体问题就是在一个热源，热流密度已知为q，在长方形（金属）的上边移动，可以说是一个移动热源。热源以外的上边以及下边、左右边的边界条件，均为换热和辐射即h(ta-tw)+...。
    求得是不同时刻这个长方形（计算区域）的温度场。
    这是非稳态的。
    我在热源处采用的是第二边界条件，附加源项法，其他边界都是第三边界条件，也是采用附加源项法。计算出各点问题。
    我用上面的一种方法计算的温度过高，我想是由于金属熔化要吸热，我还没考虑熔化潜热的影响（还不会处理这个源项^_^），只是单纯的热传导。
    但是用第二种方法计算的温度300就有问题了。
    就是这样的。我再调试调试。不知道你做过没有？谢谢了

zlf81

   我重新调试了一下第二种方法。
   用一个新数组存储t(i,j)/Cp的值，即真正的温度T，这样就不影响程序计算t(i,j)即（温度×Cp的值）。
   另外在附温度初值上我把初始温度乘以了Cp，作为计算t(i,j)的初值。
   这样计算的结果，比第一种方法更合理些。又搞不懂第一种方法哪错了^_^。
   咳，随便该点程序，就产生这么大的差别。我采用第二种方法算了^_^。
   不知道你会不会带潜热的源项处理，比如S＝d(ρ*fl*Lm)/dt。fl为液相分数，Lm为潜热。是不是把它看成Sc，并对时间t离散就行了呢。
   问的有点贪心，呵呵，谢谢了！

anyone

你有没有用简单算例验证过你的程序?
比如运行一些知道解析解的一维问题(定常,非定常, 各种边界条件)?
如果没有, 强烈建议你先用一些简单算例验证你的程序.
潜热没做过, 感觉没那么简单, 你又引进了新变量fl, 要新的方程来封闭.

anyone

刚才没看清楚, 潜热方程有了. 不过FL是未知的, 原则上可以迭代求出, 但因为温度方程在溶解温度下出现间断, 处理不好, 就会产生震荡, 收敛不了,  

miaoji

我有一个疑问啊
如果直接采用Φ=CpT，那么扩散系数Г是不是应该不用除Cp的啊
希望能指点一下
谢谢

下面引用由zlf81在 2005/04/01 09:30am 发表的内容：
你的意思是指：
    Φ=CpT，
    则，扩散系数Г还是＝k/Cp，式（1）中的源项S保持不变。
    利用simplec算出Φ以后，
...

zlf81

我觉得还是要把。
d(ρCpΦ)/dt+d(ρCpuΦ)/dx+d(ρCpvΦ)/dy=k(d^2Φ/dx^2)+k(d^2Φ/dy^2)+S  （1）
转化一下：
d(ρ（CpΦ）)/dt+d(ρu（CpΦ）)/dx+d(ρv（CpΦ）)/dy=k/Cp(d^2（CpΦ)/dx^2)+k/Cp(d^2(CpΦ)/dy^2)+S  （1）‘
令Φ＝CpΦ；
d(ρΦ)/dt+d(ρuΦ)/dx+d(ρvΦ)/dy=k/Cp(d^2Φ/dx^2)+k/Cp(d^2Φ/dy^2)+S   （2）
－》Г＝k/Cp
我觉得是这样的，你认为呢？

miaoji

对哦，没有看清楚
不好意思
我还有个问题，源项中是不是可以是求解量的偏导之类的表达式。

maximus

当然可以，只是在离散时显式的离散就是了。

zlf81

可以。
我的控制方程的源项中也含有偏导，甚至更复杂一些。
好难解啊

miaoji

你的好难解是什么意思？
是指收敛的问题吗？