[求助]充填过程温度场的求解问题

kuangtq

现要计算充填过程的温度场，为薄壁型腔，考虑厚度方向的热传导和流动方向的热对流。不知哪有问题，算出的结果很大（上千K，设定的入口温度为513K，模壁温度为313K），显然不合理。麻烦有兴趣的大大研究一下，问题在哪里？后面调用的求解程序为徐士良的常用算法。
!*****************************************************!
!  !
!  求解温度场  !
!  !
!*****************************************************!
SUBROUTINE SOLVE_T()
USE DATA_MODULE
DIMENSION JS(NNODE*NZ)
CALL GLOBAL_T_STIFFNESS()! 形成总刚
CALL IMPLEMENT_T_BOUNDARY()! 施加边界条件
CALL GGJE(GLB_T_STIF,NNODE*NZ,GLB_T_R,L,JS)! 调用全选主元高斯－约当消去法求解系数矩阵为大型稀疏矩阵的方程组
DO I=1, NNODE
DO J=1, NZ
T(I, J)=GLB_T_R((J-1)*NNODE+I)! 获取温度值
IF(T(I, J)<313 .OR. T(I, J)>513)THEN! 出错信息
WRITE(6, "('THE THERMAL IS OUT OF THE RANGE: 313-513')")
WRITE(6, "('T=', E12.6, 5X, 'NODE:', I4, 5X, 'LAYER:', I4, 5X, 'TIME_STEP:', I4)") T(I, J), I, J, NUM_TIMESTEP
STOP
ENDIF
ENDDO
ENDDO
END

!***************************************!
!  形成温度总刚 !
!***************************************!
SUBROUTINE GLOBAL_T_STIFFNESS()
USE DATA_MODULE
DIMENSION IJM(3)
INTEGER E_I! 共用节点I的单元标志
DOUBLE PRECISION XC, YC, U_XC, U_YC, VISCOS_C, M, TEMP, EL_UP_CV, ZJ
INTEGER L, N1, ROW
DOUBLE PRECISION DENSITY, HEAT_TRANS_RATE, CP
! 计算节点的上风控制体积为对流项的计算作准备
UP_CV=0
DO I=1, NNODE
DO J=1, NELEM! 遍历共用节点I的所有单元
E_I=0
DO K=1, 3
IF(EL_NODE(K, J)==I) E_I=1 ! 单元J共用节点I
ENDDO
IF(E_I==1)THEN! 计算共用节点I的单元J对I节点控制体积的贡献
XC=0
YC=0
U_XC=0
U_YC=0
DO K=1, 3! 求单元形心处的相关值
N=EL_NODE(K, J)
XC=XC+ND_XY(1, N)/3.0
YC=YC+ND_XY(2, N)/3.0
U_XC=U_XC+U_X(N)/3.0
U_YC=U_YC+U_Y(N)/3.0
ENDDO
DOT=U_XC*(ND_XY(1, I)-XC)+U_YC*(ND_XY(2, I)-YC)! 注意：熔体前沿气体区，速度场为0，因而其节点上风控制体积为0
IF(DOT>0) UP_CV(I)=UP_CV(I)+EL_BCAREA(7, J)/3.0! 计算节点的上风控制体积
ENDIF
ENDDO
ENDDO
! 计算系数矩阵
GLB_T_STIF=0! 置零
GLB_T_R=0
DO N=1, NNODE! 逐节点
DENSITY=F(1, N)*DENSITY1+(1-F(1, N))*DENSITY0! 加权计算物性参数
HEAT_TRANS_RATE=F(1, N)*HEAT_TRANS_RATE1+(1-F(1, N))*HEAT_TRANS_RATE0
CP=F(1, N)*CP1+(1-F(1, N))*CP0
M=DELTA_T/(DENSITY*CP)
DO L=1, NELEM! 遍历共用节点N的所有单元
E_I=0
DO K=1, 3
IF (EL_NODE(K, L)==N) E_I=1
ENDDO
IF(E_I==1)THEN
! 热传导项系数  涉及到相邻两层
DO K1=1, 3
N1=EL_NODE(K1, L)
TEMP=M/ND_CV(N)*(EL_BCAREA(7, L)/3.0)*HEAT_TRANS_RATE/(3.0*DELTA_Z*DELTA_Z)
DO JJ=1, NZ! 逐层
J=JJ-1
ROW=J*NNODE+N! J*NNODE+N   N节点在第LAYER层的差分节点的总体编号, ROW方程行号
IF(JJ==1 .OR. JJ==NZ)THEN! 贴模壁层
IF(JJ==1)THEN! J=0 无 J-1
GLB_T_STIF(ROW, (J+1)*NNODE+N1)=GLB_T_STIF(ROW, (J+1)*NNODE+N1)-2*TEMP
ELSE! J=NZ-1 无 J+1
GLB_T_STIF(ROW, (J-1)*NNODE+N1)=GLB_T_STIF(ROW, (J-1)*NNODE+N1)-2*TEMP
ENDIF
GLB_T_STIF(ROW, J*NNODE+N1)=GLB_T_STIF(ROW, J*NNODE+N1)+2*(1+DELTA_Z*HEAT_TRANSF_W/HEAT_TRANS_RATE)*TEMP
GLB_T_R(ROW)=GLB_T_R(ROW)+2*DELTA_Z*HEAT_TRANSF_W*T_W*TEMP/HEAT_TRANS_RATE! 只有贴模壁层有此项
ELSE! 内部层
GLB_T_STIF(ROW, (J+1)*NNODE+N1)=GLB_T_STIF(ROW, (J+1)*NNODE+N1)-TEMP
GLB_T_STIF(ROW, J*NNODE+N1)=GLB_T_STIF(ROW, J*NNODE+N1)+2*TEMP
GLB_T_STIF(ROW, (J-1)*NNODE+N1)=GLB_T_STIF(ROW, (J-1)*NNODE+N1)-TEMP
ENDIF
ENDDO
ENDDO
! 计算节点的上风单元控制体积贡献，为计算热对流项而准备
XC=0
YC=0
U_XC=0
U_YC=0
VISCOS_C=0
DO K2=1, 3
N1=EL_NODE(K2, L)
XC=XC+ND_XY(1, N1)/3.0
YC=YC+ND_XY(2, N1)/3.0
U_XC=U_XC+U_X(N1)/3.0
U_YC=U_YC+U_Y(N1)/3.0
VISCOS_C=VISCOS_C+VISCOS(N1)/3.0
ENDDO
DOT=U_XC*(ND_XY(1, N)-XC)+U_YC*(ND_XY(2, N)-YC)
IF(DOT>0) THEN
EL_UP_CV=EL_BCAREA(7, L)/3.0
ELSE
EL_UP_CV=0
ENDIF
DO JJ=1, NZ! 逐层
J=JJ-1
ROW=J*NNODE+N! J*NNODE+N   N节点在第LAYER层的差分节点的总体编号, ROW方程行号
! 热对流项系数 仅考虑来自节点上游单元的贡献
ZJ=(J-0.5)*DELTA_Z-HALF_THICKNESS! J层的Z坐标值(中性面Z坐标为0)
IF(UP_CV(N)/=0)THEN! 熔体区
TEMP=DELTA_T/(UP_CV(N))*EL_UP_CV*((ZJ*ZJ-HALF_THICKNESS*HALF_THICKNESS)/(2*VISCOS_C))/(2*EL_BCAREA(7, L))
ELSE! 熔体前沿及气体区，认为无热传导
TEMP=0
ENDIF
DO K3=1, 3
N1=EL_NODE(K3, L)
GLB_T_STIF(ROW, J*NNODE+N1)=GLB_T_STIF(ROW, J*NNODE+N1)+TEMP*(-DP_DX_C(L)*EL_BCAREA(K3, L)+(-DP_DY_C(L)*EL_BCAREA(K3+3, L)))
ENDDO
! 瞬态项系数为1
GLB_T_STIF(ROW, ROW)=1
! 右端项
GLB_T_R(ROW)=GLB_T_R(ROW)+T(N, JJ)! +M*STREAR_HEAT  STREAR_HEAT为剪切热
ENDDO
ENDIF
ENDDO
ENDDO
END
!***************************************!
!施加温度边界条件 !
!***************************************!
SUBROUTINE IMPLEMENT_T_BOUNDARY()
USE DATA_MODULE
DO J=1, NZ
N=(J-1)*NNODE+NODE_ENTRY! 入口处各差分节点的总编号
GLB_T_STIF(N, N)=GLB_T_STIF(N, N)*1E20! 入口节点 熔体入口温度，对角元素乘大数
GLB_T_R(N)=GLB_T_STIF(N, N)*T_ENTRY! 右端项
ENDDO
DO I=1, NNODE
IF(F(1, I)==0) THEN! 熔体前沿气体区  温度为模温
DO J=1, NZ
N=(J-1)*NNODE+I
GLB_T_STIF(N, N)=GLB_T_STIF(N, N)*1E20
GLB_T_R(N)=GLB_T_STIF(N, N)*T_W
ENDDO
ENDIF
ENDDO
DO I=1, NZ*NNODE! 检查总刚是否正定
IF(GLB_T_STIF(I, I)<=0)THEN
WRITE(6, "(/10X, 'MAIN ELEMENT OF GLB_T_STIF IS NONPOSITIVE. N=', I4, 5X, E12.6)") I, GLB_T_STIF(I, I)
STOP
ENDIF
ENDDO
END
!***************************************!
!全选主元高斯－约当消去法 !
!***************************************!
SUBROUTINE GGJE(A,N,B,L,JS)! 引自《Fortran常用算法程序集-徐士良》AGGJE.FOR
DIMENSION A(N, N), B(N), JS(N)! 全选主元高斯－约当消去法求解系数矩阵为大型稀疏矩阵的方程组
DOUBLE PRECISION A, B
DOUBLE PRECISION TEMP

L=1
DO K=1,N
IF(K==48)THEN
KK=1
ENDIF
D=0.0
DO I=K,N
DO J=K,N
IF (ABS(A(I,J)).GT.D) THEN
D=ABS(A(I,J))
JS(K)=J
IS=I
END IF
ENDDO
ENDDO
IF (D+1.0.EQ.1.0) THEN
WRITE(*,20)
L=0
RETURN
END IF
20FORMAT(1X,'  FAIL  ')
DO J=K, N
TEMP=A(K, J)
A(K, J)=A(IS, J)
A(IS, J)=TEMP
ENDDO
TEMP=B(K)
B(K)=B(IS)
B(IS)=TEMP
DO I=1, N
TEMP=A(I,K)
A(I, K)=A(I,JS(K))
A(I, JS(K))=TEMP
ENDDO
TEMP=A(K,K)
DO J=K+1,N
IF (A(K,J).NE.0.0) A(K,J)=A(K,J)/TEMP
ENDDO
B(K)=B(K)/TEMP
DO J=K+1,N
IF (A(K,J).NE.0.0) THEN
DO I=1,N
IF ((I.NE.K).AND.(A(I,K).NE.0.0)) THEN
A(I,J)=A(I,J)-A(I,K)*A(K,J)
END IF
ENDDO
END IF
ENDDO
DO I=1,N
IF ((I.NE.K).AND.(A(I,K).NE.0.0)) THEN
B(I)=B(I)-A(I,K)*B(K)
END IF
ENDDO
ENDDO
DO K=N,1,-1
IF (K.NE.JS(K)) THEN
TEMP=B(K)
B(K)=B(JS(K))
B(JS(K))=TEMP
END IF
ENDDO
RETURN
END