帮忙看看为什么分析解和LBM计算的剪切应力不同呢?

healthykxy

剪切应力计算问题
LBM模拟参数：沿X 方向的2D poiseuille流，入口和出口压强边界条件rou_in=1.03，rou_out=1.00，上下反弹边界，计算区域长Lx =60，宽Ly =20，（所有量都是晶格单位），tau=0.9，在LBM中计算剪切应力采用的公式为：
σxy =– (1-1/2τ) ∑_(i=0)^8((f_i-f_i^eq ) e_ix*e_iy )
分析解：
y0=dble(Ly-2)/2.d0+1.5d0
gradp=(rou_in-rou_out)*c_squ/dble(Lx)
u_max=gradp/2.d0/((2.d0*tau-1.d0)/6.d0)*(dble(Ly-2)/2.d0)**2.d0
r=abs(y0-iy)
stress=u_max*2.d0*r/((dble(Ly-2)/2.d0)**2.d0)*(2.d0*tau-1.d0)/6.d0
在中间取一横截面，比较两种方法计算剪切应力结果如图