diff --git a/Source/ccib.f90 b/Source/ccib.f90 index 751786d8d04..02067cc5704 100644 --- a/Source/ccib.f90 +++ b/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 --------------------------------- diff --git a/Source/wall.f90 b/Source/wall.f90 index e54a85035d9..ecec1d3a324 100644 --- a/Source/wall.f90 +++ b/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)