Dtransu の流速 その2

飽和時の計算について、マニュアルには何も書かれていませんでした。


仕方ないので、ソースを見てみました。

ここですね。


C     DETERMINE RELATIVE HYDRAULIC CONDUCTIVITY
C
        LTB = MTBL3D( MAT )
        DO 61 K=1,NN
          I = LM(K)
C
c-------- confined problem ----

          if( kan2.eq.0 ) then

            CR(I) = 1.0
            TXX(K) = 100.0D0
            go to 61
          end if
c
          IF( P(I).LT.0.0 ) THEN
          CALL INTERP(XYP( 1,NMTB(1,2,LTB)+1 ),TX,P(I),NMTB(2,2,LTB),1)
          TXX(K) = TX
          IF( CR(I).LE.0.0 )
     2    CALL INTERP(XYK( 1,NMTB(1,1,LTB)+1 ),TX,CR(I),NMTB(2,1,LTB),0)
          ELSE
          CR(I) = 1.0
          TXX(K) = 100.0D0
          END IF
C
   61   CONTINUE



      SUBROUTINE SETCR( KODE,P,P1,CR,T,CS,POR,XYK,NK,XYP,NP,XYC,NC )
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      COMMON / CTPRM2/ KAN1,KAN2
C
      DIMENSION XYK(2,NK),XYP(2,NP),XYC(2,NC)
C
      IF( CR.GE.0.0D0 ) RETURN
      IF( KAN2.EQ.0 ) GO TO 2000
      IF( P.LT.0.0D0 )  GO TO 1000
C
C---- SATURATED ---------------------------
C
 2000 CR = 1.0D0
      T  = POR
      CS = 0.0D0
C
      RETURN
C
C---- UNSATURATED -------------------------
C
 1000 IF( KODE.LT.1 ) CM = P
      IF( KODE.GE.1 ) CM = (P+P1)*0.5D0
      CALL INTERP( XYP,T,CM,NP,1 )
      CALL INTERP( XYK,T,CR,NK,0 )
      CALL INTERP( XYC,T,CS,NC,0 )
      T = T*POR*0.01D0
      CS= CS*POR
C
      RETURN
C
C
      END




コメントがないので誤っているかもしれませんが、飽和のフラグを立てると、比透水係数1、体積含水率 = 有効間隙率で計算しているのでしょう。単純に、この条件で負圧となれば不飽和時よりも上昇流や水平方向の流れが発生しやすいということでしょうね。


結局、変換ツールのソースを、圧力水頭がマイナス時は、Vi ⇔ { 0, 0, -0.1} として吐き出すように修正しました。不飽和の部分では鉛直下降流が起こっている、という、よくある簡略化です。

そうすると、TECPLOT の streamline が、思い描いたように綺麗になりました。


うーん、いいのか？