|
楼主 |
发表于 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
|
|