找回密码
 注册
查看: 3581|回复: 5

代数应力模型

[复制链接]
发表于 2003-8-6 16:26:06 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?注册

x
偶自编程序计算分层流,希望用代数应力模型,但是根据模型的实现代码。
代码本人经过无数遍确认,保证无误。不过没有特殊处理。
但是就是的不到想要的结果。而且计算的主应力值之和不等于 2*k.
请有经验的同道们救火!
发表于 2003-8-7 10:11:25 | 显示全部楼层

代数应力模型

俺也准备用代数应力模型,但还没开始写程序。
 楼主| 发表于 2003-8-7 14:08:54 | 显示全部楼层

代数应力模型

我是用W.Rodi提出的方程。程序代码很简单,就是只按书中的关系编写的程序算出来结果有问题,等晚上我把程序贴上来,大家讨论讨论。
我估计是由于数值计算中因连续性要求不能完全满足和梯度估算问题产生的,所以可能根据可实现性,如许瓦兹不等式等,需要作一些限制。
 楼主| 发表于 2003-8-7 20:07:15 | 显示全部楼层

代数应力模型

  SUBROUTINE COMASM1( )
INTEGER :: IC
REAL*8 :: A(6,6),B(6),X(6),GRADU(3),GRADV(3),GRADW(3)
REAL*8 :: T1,T11,T12,D1,D2,D3,D4,D5,D6, BETAG1, BETAG2, THMC, TTS
DO IC=1,NCELLS
  TTS = TKE(iC)/(EDR(iC)+SMALL)
  IF (TTS.LT.0) THEN
TTS = 0
  ELSE IF(TTS.GT.100)THEN
TTS=100
  ENDIF
GRADU = DUDXI(IC,
GRADV = DVDXI(IC,
GRADW = DWDXI(IC,
!IF(IC.EQ.5342)THEN
!write(*,*)ic
!ENDIF
IF(SYSBY(6).NE.0)THEN
BETAG1 = BETA(THM(IC))/den(ic)
!BETAG1 = BETA(THM(IC))/(0.999842594 + BETA(THMENV)*THMENV)
ENDIF
IF(SYSBY(7).NE.0)THEN
IF(SYSBY(6).NE.0)THEN
  THMC = THM(IC)
ELSE
  THMC = THMENV
END IF
BETAG2 = ETA(THMC,PLT(IC))/den(ic)!(0.999842594 + BETA(THMENV)*THMENV)
ENDIF
  T1 = -TTS/C1
  T11 = 2*(1-C2)
  T12 = 2*C2/3.0
  D1 = T1 * ( T11+T12 ) * GRADU(1)
  D2 = T1 * ( T11*GRADU(2)+T12*(GRADU(2)+GRADV(1)) )
  D3 = T1 * ( T11*GRADU(3)+T12*(GRADU(3)+GRADW(1)) )
  D4 = T1 * T12 * GRADV(2)
  D5 = T1 * T12 * ( GRADV(3) + GRADW(2) )
  D6 = T1 * T12 * GRADW(3)
  A(1,1) = 1-D1
  A(1,2) =  -D2  
  A(1,3) =  -D3   
  A(1,4) =  -D4   
  A(1,5) =  -D5   
  A(1,6) =  -D6
  B(1) = T1*2*(1-C1)/3*EDR(IC)
!  IF(SYSBY(9).NE.0)THEN
IF(SYSBY(6).NE.0)THEN
  B(1)=B(1)+T1*(2.*C3/3)*BETAG1*GRAVITY*WTA(IC)
ENDIF
IF(SYSBY(7).NE.0)THEN
  B(1)=B(1)+T1*(2.*C3/3)*BETAG2*GRAVITY*WSA(IC)
ENDIF
!  END IF
  T1 = -(1-C2)*TTS/C1
  D1 = T1 *  GRADV(1)
  D2 = T1 * ( GRADU(1)+GRADV(2) )
  D3 = T1 * GRADV(3)
  D4 = T1 * GRADU(2)
  D5 = T1 * GRADU(3)
  D6 = 0
  A(2,1) =  -D1
  A(2,2) =  1-D2  
  A(2,3) =  -D3   
  A(2,4) =  -D4   
  A(2,5) =  -D5   
  A(2,6) =  -D6
  B(2) = 0
  T1 = -(1-C2)*TTS/C1
  D1 = T1 * GRADW(1)
  D2 = T1 * GRADW(2)
  D3 = T1 * ( GRADU(1)+GRADW(3) )
  D4 = 0
  D5 = T1 * GRADU(2)
  D6 = T1 * GRADU(3)
  A(3,1) =  -D1
  A(3,2) =  -D2  
  A(3,3) =  1-D3   
  A(3,4) =  -D4   
  A(3,5) =  -D5   
  A(3,6) =  -D6
  B(3) = 0
!  IF(SYSBY(9).NE.0)THEN
IF(SYSBY(6).NE.0)THEN
  B(3) = B(3) - (1-C3)*TTS/C1*BETAG1*GRAVITY*UTA(IC)
ENDIF
IF(SYSBY(7).NE.0)THEN
  B(3) = B(3) - (1-C3)*TTS/C1*BETAG2*GRAVITY*USA(IC)
ENDIF
!  END IF
  T1 = -TTS/C1
  T11 = 2*(1-C2)
  T12 = 2*C2/3.0
  D1 = T1 * T12 * GRADU(1)
  D2 = T1 * ( T11*GRADV(1) + T12*(GRADU(2)+GRADV(1)) )
  D3 = T1 * T12 * ( GRADU(3)+GRADW(1) )
  D4 = T1 * ( T11+T12 )*GRADV(2)
  D5 = T1 * ( T11*GRADV(3) + T12*(GRADV(3)+GRADW(2)) )
  D6 = T1 * T12 * GRADW(3)
  A(4,1) =  -D1
  A(4,2) =  -D2  
  A(4,3) =  -D3   
  A(4,4) =  1-D4   
  A(4,5) =  -D5   
  A(4,6) =  -D6
  B(4) = T1*2*(1-C1)*EDR(IC)/3
!  IF(SYSBY(9).NE.0)THEN
IF(SYSBY(6).NE.0)THEN
  B(4) = B(4) + T1*2*C3/3*BETAG1*GRAVITY*WTA(IC)
ENDIF
IF(SYSBY(7).NE.0)THEN
  B(4) = B(4) + T1*2*C3/3*BETAG2*GRAVITY*WSA(IC)
ENDIF
!  END IF
  T1 = -(1-C2)*TTS/C1
  D1 = 0
  D2 = T1 * GRADW(1)
  D3 = T1 * GRADV(1)
  D4 = T1 * GRADW(2)
  D5 = T1 * ( GRADV(2)+GRADW(3) )
  D6 = T1 * GRADV(3)
  A(5,1) =  -D1
  A(5,2) =  -D2  
  A(5,3) =  -D3   
  A(5,4) =  -D4   
  A(5,5) =  1-D5   
  A(5,6) =  -D6
  B(5) = 0
!  IF(SYSBY(9).NE.0)THEN
IF(SYSBY(6).NE.0)THEN
  B(5) = B(5) - TTS/C1*(1-C3)*BETAG1*GRAVITY*VTA(IC)
ENDIF
IF(SYSBY(7).NE.0)THEN
  B(5) = B(5) - TTS/C1*(1-C3)*BETAG2*GRAVITY*VSA(IC)
ENDIF
!  END IF

  T1 = -TTS/C1
  T11 = 2*(1-C2)
  T12 = 2*C2/3.0
  D1 = T1 * T12 * GRADU(1)
  D2 = T1 * T12 * ( GRADU(2)+GRADV(1) )
  D3 = T1 * ( T11*GRADW(1) + T12*(GRADU(3)+GRADW(1)) )
  D4 = T1 * T12 * GRADV(2)
  D5 = T1 * ( T11*GRADW(2) + T12*(GRADV(3)+GRADW(2)) )
  D6 = T1 * ( T11+T12 ) * GRADW(3)
  A(6,1) =  -D1
  A(6,2) =  -D2  
  A(6,3) =  -D3   
  A(6,4) =  -D4   
  A(6,5) =  -D5   
  A(6,6) =  1-D6
  B(6) = 2*T1*(1-C1)*EDR(IC)/3
!  IF(SYSBY(9).NE.0)THEN
IF(SYSBY(6).NE.0)THEN
  B(6) = B(6) + 2*T1*(1-2.*C3/3)*BETAG1*GRAVITY*WTA(IC)
ENDIF
IF(SYSBY(7).NE.0)THEN
  B(6) = B(6) + 2*T1*(1-2.*C3/3)*BETAG2*GRAVITY*WSA(IC)
ENDIF
!  END IF

  CALL DLSARG (6, A, 6, B, 1, X)
  UUA(IC) = X(1)
  UVA(IC) = X(2)
  UWA(IC) = X(3)
  VVA(IC) = X(4)
  VWA(IC) = X(5)
  WWA(IC) = X(6)
END DO
  END SUBROUTINE
 楼主| 发表于 2003-8-7 20:12:09 | 显示全部楼层

代数应力模型

这是计算程序,希望大家发表意见。
SYSBY(I) 为变量计算控制参数,=1---计算该变量, =0 ---N
       I=6,温度        =7,为浓度
发表于 2003-11-10 16:39:13 | 显示全部楼层

代数应力模型

可压不可压?
您需要登录后才可以回帖 登录 | 注册

本版积分规则

快速回复 返回顶部 返回列表