Merge pull request #13316 from marcosvanella/master

FDS Source : Issue #13235, move wall model call to CFACE loop in wall…

Source/ccib.f90

@@ -92,7 +92,6 @@ MODULE CC_SCALARS
 TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC
 TYPE(EXTERNAL_WALL_TYPE),  POINTER :: EWC
 TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1
-TYPE(BOUNDARY_PROP2_TYPE), POINTER :: B2
 
 ! Baroclininc torque linking variables.
 REAL(EB), SAVE, ALLOCATABLE, DIMENSION(:) :: F_LINK,A_LINK
@@ -5798,7 +5797,6 @@ END SUBROUTINE CC_MATCH_VELOCITY
 SUBROUTINE CC_COMPUTE_VISCOSITY(DT,NM)
 
 USE TURBULENCE, ONLY: WALE_VISCOSITY
-USE TURBULENCE, ONLY : WALL_MODEL
 
 REAL(EB), INTENT(IN):: DT
 INTEGER, INTENT(IN) :: NM
@@ -5809,11 +5807,7 @@ SUBROUTINE CC_COMPUTE_VISCOSITY(DT,NM)
 REAL(EB), POINTER, DIMENSION(:,:,:,:) :: ZZP=>NULL()
 REAL(EB) :: NU_EDDY,DELTA,A_IJ(3,3),DUDX,DUDY,DUDZ,DVDX,DVDY,DVDZ,DWDX,DWDY,DWDZ
 
-REAL(EB) :: NVEC(IAXIS:KAXIS), TVEC1(IAXIS:KAXIS), TVEC2(IAXIS:KAXIS), VEL_WALL(IAXIS:KAXIS), VEL_CELL(IAXIS:KAXIS), &
-            DN, RHO_WALL, MU_WALL, SLIP_FACTOR, TT(IAXIS:KAXIS), SS(IAXIS:KAXIS), &
-            U_NORM, U_ORTH, U_STRM, U_CELL, V_CELL, W_CELL, U_RELA(IAXIS:KAXIS), TNOW
-INTEGER  :: ICF,IND1,IND2
-TYPE(SURFACE_TYPE), POINTER :: SF
+REAL(EB) :: TNOW
 
 ! Dummy assignments:
 TNOW = DT
@@ -5871,45 +5865,6 @@ SUBROUTINE CC_COMPUTE_VISCOSITY(DT,NM)
 
 ENDIF LES_IF
 
-! Now compute U_TAU and Y_PLUS on CFACES:
-DO ICF=INTERNAL_CFACE_CELLS_LB+1,INTERNAL_CFACE_CELLS_LB+N_INTERNAL_CFACE_CELLS
-
-   CFA => CFACE(ICF)
-   BC => BOUNDARY_COORD(CFA%BC_INDEX)
-   B1 => BOUNDARY_PROP1(CFA%B1_INDEX)
-   B2 => BOUNDARY_PROP2(CFA%B2_INDEX)
-   SF => SURFACE(CFA%SURF_INDEX)
-
-   ! Surface Velocity:
-   NVEC = BC%NVEC
-
-   ! right now VEL_T not defined for CFACEs
-   TVEC1=(/ 0._EB,0._EB,0._EB/)
-   TVEC2=(/ 0._EB,0._EB,0._EB/)
-   ! velocity vector of the surface
-   VEL_WALL = -B1%U_NORMAL*NVEC + SF%VEL_T(1)*TVEC1 + SF%VEL_T(2)*TVEC2
-
-   ! find cut-cell adjacent to CFACE
-   IND1 = CFA%CUT_FACE_IND1
-   IND2 = CFA%CUT_FACE_IND2
-   CALL GET_UVWGAS_CFACE(U_CELL,V_CELL,W_CELL,IND1,IND2)
-   ! velocity vector in the centroid of the gas (cut) cell
-   VEL_CELL = (/U_CELL,V_CELL,W_CELL/) ! (/1._EB,0._EB,0._EB/) ! test
-
-   CALL GET_MUDNS_CFACE(MU_WALL,IND1,IND2)
-
-   ! Gives local velocity components U_STRM , U_ORTH , U_NORM
-   ! in terms of unit vectors SS,TT,NN:
-   U_RELA(IAXIS:KAXIS) = VEL_CELL(IAXIS:KAXIS)-VEL_WALL(IAXIS:KAXIS)
-   CALL GET_LOCAL_VELOCITY(U_RELA,NVEC,TT,SS,U_NORM,U_ORTH,U_STRM)
-
-   DN  = 1._EB/B1%RDN
-   RHO_WALL = B1%RHO_F
-
-   CALL WALL_MODEL(SLIP_FACTOR,B2%U_TAU,B2%Y_PLUS,MU_WALL/RHO_WALL,SF%ROUGHNESS,0.5_EB*DN,U_STRM)
-
-ENDDO
-
 T_USED(14) = T_USED(14) + CURRENT_TIME() - TNOW
 IF (TIME_CC_IBM) &
    T_CC_USED(CC_COMPUTE_VISCOSITY_TIME_INDEX) = T_CC_USED(CC_COMPUTE_VISCOSITY_TIME_INDEX) + CURRENT_TIME() - TNOW
@@ -15267,46 +15222,46 @@ END SUBROUTINE GET_LINKED_FV
 ! END SUBROUTINE DEBUG_WAIT
 ! #endif /* defined(DEBUG_CC_INTERPOLATION) */
 
-! ---------------------------- GET_LOCAL_VELOCITY ------------------------------
+! ! ---------------------------- GET_LOCAL_VELOCITY ------------------------------
 
-SUBROUTINE GET_LOCAL_VELOCITY(U_RELA,NN,TT,SS,U_NORM,U_ORTH,U_STRM)
-USE MATH_FUNCTIONS, ONLY: CROSS_PRODUCT
+! SUBROUTINE GET_LOCAL_VELOCITY(U_RELA,NN,TT,SS,U_NORM,U_ORTH,U_STRM)
+! USE MATH_FUNCTIONS, ONLY: CROSS_PRODUCT
 
-REAL(EB), INTENT(IN) :: NN(IAXIS:KAXIS), U_RELA(IAXIS:KAXIS)
-REAL(EB), INTENT(OUT):: TT(IAXIS:KAXIS), SS(IAXIS:KAXIS), U_NORM, U_ORTH, U_STRM
+! REAL(EB), INTENT(IN) :: NN(IAXIS:KAXIS), U_RELA(IAXIS:KAXIS)
+! REAL(EB), INTENT(OUT):: TT(IAXIS:KAXIS), SS(IAXIS:KAXIS), U_NORM, U_ORTH, U_STRM
 
-! Local Variables:
-REAL(EB), DIMENSION(3), PARAMETER :: E1=(/1._EB,0._EB,0._EB/),E2=(/0._EB,1._EB,0._EB/),E3=(/0._EB,0._EB,1._EB/)
-REAL(EB), DIMENSION(3,3) :: C
-
-! find a vector TT in the tangent plane of the surface and orthogonal to U_VELO-U_SURF
-CALL CROSS_PRODUCT(TT,NN,U_RELA) ! TT = NN x U_RELA
-IF (ABS(NORM2(TT))<=TWO_EPSILON_EB) THEN
-   ! tangent vector is completely arbitrary, just perpendicular to NN
-   IF (ABS(NN(1))>=TWO_EPSILON_EB .OR.  ABS(NN(2))>=TWO_EPSILON_EB) TT = (/NN(2),-NN(1),0._EB/)
-   IF (ABS(NN(1))<=TWO_EPSILON_EB .AND. ABS(NN(2))<=TWO_EPSILON_EB) TT = (/NN(3),0._EB,-NN(1)/)
-ENDIF
-TT = TT/NORM2(TT) ! normalize to unit vector
-CALL CROSS_PRODUCT(SS,TT,NN) ! define the streamwise unit vector SS
-
-! directional cosines (see Pope, Eq. A.11)
-C(1,1) = DOT_PRODUCT(E1,SS)
-C(1,2) = DOT_PRODUCT(E1,TT)
-C(1,3) = DOT_PRODUCT(E1,NN)
-C(2,1) = DOT_PRODUCT(E2,SS)
-C(2,2) = DOT_PRODUCT(E2,TT)
-C(2,3) = DOT_PRODUCT(E2,NN)
-C(3,1) = DOT_PRODUCT(E3,SS)
-C(3,2) = DOT_PRODUCT(E3,TT)
-C(3,3) = DOT_PRODUCT(E3,NN)
-
-! transform velocity (see Pope, Eq. A.17)
-U_STRM = C(1,1)*U_RELA(1) + C(2,1)*U_RELA(2) + C(3,1)*U_RELA(3)
-U_ORTH = C(1,2)*U_RELA(1) + C(2,2)*U_RELA(2) + C(3,2)*U_RELA(3)
-U_NORM = C(1,3)*U_RELA(1) + C(2,3)*U_RELA(2) + C(3,3)*U_RELA(3)
+! ! Local Variables:
+! REAL(EB), DIMENSION(3), PARAMETER :: E1=(/1._EB,0._EB,0._EB/),E2=(/0._EB,1._EB,0._EB/),E3=(/0._EB,0._EB,1._EB/)
+! REAL(EB), DIMENSION(3,3) :: C
 
-RETURN
-END SUBROUTINE GET_LOCAL_VELOCITY
+! ! find a vector TT in the tangent plane of the surface and orthogonal to U_VELO-U_SURF
+! CALL CROSS_PRODUCT(TT,NN,U_RELA) ! TT = NN x U_RELA
+! IF (ABS(NORM2(TT))<=TWO_EPSILON_EB) THEN
+!    ! tangent vector is completely arbitrary, just perpendicular to NN
+!    IF (ABS(NN(1))>=TWO_EPSILON_EB .OR.  ABS(NN(2))>=TWO_EPSILON_EB) TT = (/NN(2),-NN(1),0._EB/)
+!    IF (ABS(NN(1))<=TWO_EPSILON_EB .AND. ABS(NN(2))<=TWO_EPSILON_EB) TT = (/NN(3),0._EB,-NN(1)/)
+! ENDIF
+! TT = TT/NORM2(TT) ! normalize to unit vector
+! CALL CROSS_PRODUCT(SS,TT,NN) ! define the streamwise unit vector SS
+
+! ! directional cosines (see Pope, Eq. A.11)
+! C(1,1) = DOT_PRODUCT(E1,SS)
+! C(1,2) = DOT_PRODUCT(E1,TT)
+! C(1,3) = DOT_PRODUCT(E1,NN)
+! C(2,1) = DOT_PRODUCT(E2,SS)
+! C(2,2) = DOT_PRODUCT(E2,TT)
+! C(2,3) = DOT_PRODUCT(E2,NN)
+! C(3,1) = DOT_PRODUCT(E3,SS)
+! C(3,2) = DOT_PRODUCT(E3,TT)
+! C(3,3) = DOT_PRODUCT(E3,NN)
+
+! ! transform velocity (see Pope, Eq. A.17)
+! U_STRM = C(1,1)*U_RELA(1) + C(2,1)*U_RELA(2) + C(3,1)*U_RELA(3)
+! U_ORTH = C(1,2)*U_RELA(1) + C(2,2)*U_RELA(2) + C(3,2)*U_RELA(3)
+! U_NORM = C(1,3)*U_RELA(1) + C(2,3)*U_RELA(2) + C(3,3)*U_RELA(3)
+
+! RETURN
+! END SUBROUTINE GET_LOCAL_VELOCITY
 
 
 ! ------------------------------- CC_NO_FLUX ---------------------------------

Source/wall.f90

@@ -166,8 +166,14 @@ SUBROUTINE WALL_BC(T,DT,NM)
 
    IF (N_TRACKED_SPECIES>1) CALL CALCULATE_RHO_D_F(B1,BC,CFACE_INDEX=ICF)
 
+   IF (CFA%BOUNDARY_TYPE==SOLID_BOUNDARY .AND. &
+      (ANY(SPECIES_MIXTURE%CONDENSATION_SMIX_INDEX>0) .OR. DEPOSITION .OR. OUTPUT_WALL_QUANTITIES)) THEN
+      CALL WALL_MODEL(SLIP_COEF,B2%U_TAU,B2%Y_PLUS,MU_DNS(BC%IIG,BC%JJG,BC%KKG)/RHO(BC%IIG,BC%JJG,BC%KKG),SF%ROUGHNESS,&
+                      0.5_EB/B1%RDN,B1%U_TANG)
+   ENDIF
+
    IF (DEPOSITION .AND. .NOT.INITIALIZATION_PHASE .AND. CORRECTOR .AND. .NOT.SOLID_PHASE_ONLY) THEN
-      IF (WC%BOUNDARY_TYPE==SOLID_BOUNDARY .AND. B1%NODE_INDEX==0 .AND. ABS(SF%VEL)<TWO_EPSILON_EB .AND. &
+      IF (CFA%BOUNDARY_TYPE==SOLID_BOUNDARY .AND. B1%NODE_INDEX==0 .AND. ABS(SF%VEL)<TWO_EPSILON_EB .AND. &
           ANY(ABS(SF%MASS_FLUX)<TWO_EPSILON_EB) .AND. ABS(SF%VOLUME_FLOW)<TWO_EPSILON_EB) THEN
          CALL CALC_DEPOSITION(DT,BC,B1,B2,CFACE_INDEX=ICF)
       ENDIF
@@ -2158,7 +2164,7 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX,
    ENDIF PYROLYSIS_PREDICTED_IF
 
    ! Add internal heat source specified by user
-   
+
    DO I=1,NWP
       Q_S(I) = Q_S(I) + ONE_D%HEAT_SOURCE(LAYER_INDEX(I))*EVALUATE_RAMP(T-T_BEGIN,ONE_D%RAMP_IHS_INDEX(LAYER_INDEX(I)))
    ENDDO