代数应力模型

ljin

偶自编程序计算分层流，希望用代数应力模型，但是根据模型的实现代码。
代码本人经过无数遍确认，保证无误。不过没有特殊处理。
但是就是的不到想要的结果。而且计算的主应力值之和不等于 2*k.
请有经验的同道们救火！

maximus

俺也准备用代数应力模型，但还没开始写程序。

ljin

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

ljin

  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

ljin

这是计算程序，希望大家发表意见。
SYSBY（I） 为变量计算控制参数，=1---计算该变量， =0 ---N
       I=6，温度        =7，为浓度

xyyr

可压不可压？