From 46f69b726ac4100f9716c9b337e7e1c62488b840 Mon Sep 17 00:00:00 2001 From: mcgratta Date: Tue, 16 Apr 2024 17:06:45 -0400 Subject: [PATCH] FDS Source: Continue streamlining wall.f90 --- Source/wall.f90 | 536 ++++++++++++++++++------------------------------ 1 file changed, 203 insertions(+), 333 deletions(-) diff --git a/Source/wall.f90 b/Source/wall.f90 index 47b931a04a6..b20ea51f74c 100644 --- a/Source/wall.f90 +++ b/Source/wall.f90 @@ -27,55 +27,11 @@ MODULE WALL_ROUTINES SUBROUTINE WALL_BC(T,DT,NM) USE COMP_FUNCTIONS, ONLY: CURRENT_TIME +USE CC_SCALARS, ONLY : CFACE_THERMAL_GASVARS REAL(EB) :: TNOW REAL(EB), INTENT(IN) :: T,DT INTEGER, INTENT(IN) :: NM LOGICAL :: CALL_HT_1D - -IF (LEVEL_SET_MODE==1) RETURN ! No need for boundary conditions if the simulation is uncoupled fire spread only - -TNOW=CURRENT_TIME() - -CALL POINT_TO_MESH(NM) - -! Compute the temperature TMP_F at all boundary cells, including PYROLYSIS and 1-D heat transfer - -CALL THERMAL_BC(T,NM,CALL_HT_1D) - -! Compute rho*D at WALL cells - -CALL DIFFUSIVITY_BC - -! Special boundary routines -IF (DEPOSITION .AND. .NOT.INITIALIZATION_PHASE .AND.CORRECTOR .AND. .NOT.SOLID_PHASE_ONLY) CALL DEPOSITION_BC(DT) -IF (HVAC_SOLVE .AND. .NOT.INITIALIZATION_PHASE) CALL HVAC_BC - -! Compute the species mass fractions, ZZ_F, at all boundary cells - -CALL SPECIES_BC(T,DT,NM,CALL_HT_1D) - -! Compute the density, RHO_F, at WALL cells only - -CALL DENSITY_BC - -T_USED(6)=T_USED(6)+CURRENT_TIME()-TNOW -END SUBROUTINE WALL_BC - - -!> \brief Thermal boundary conditions for all boundaries. -!> \param T Time (s) -!> \param NM Mesh number -!> \param CALL_HT_1D Indicator if 1-D conduction solver is called -!> \details One dimensional heat transfer and pyrolysis is done in PYROLYSIS. -!> Note also that gas phase values are assigned here to be used for all subsequent BCs. -!> \callgraph - -SUBROUTINE THERMAL_BC(T,NM,CALL_HT_1D) - -USE CC_SCALARS, ONLY : CFACE_THERMAL_GASVARS -REAL(EB), INTENT(IN) :: T -INTEGER, INTENT(IN) :: NM -LOGICAL, INTENT(OUT) :: CALL_HT_1D REAL(EB) :: DT_BC INTEGER :: IW,IP,ICF,ITW TYPE(WALL_TYPE), POINTER :: WC @@ -85,6 +41,13 @@ SUBROUTINE THERMAL_BC(T,NM,CALL_HT_1D) TYPE(LAGRANGIAN_PARTICLE_CLASS_TYPE), POINTER :: LPC TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1 +TYPE(BOUNDARY_PROP2_TYPE), POINTER :: B2 + +IF (LEVEL_SET_MODE==1) RETURN ! No need for boundary conditions if the simulation is uncoupled fire spread only + +TNOW=CURRENT_TIME() + +CALL POINT_TO_MESH(NM) IF (PREDICTOR) THEN UU => US @@ -105,8 +68,8 @@ SUBROUTINE THERMAL_BC(T,NM,CALL_HT_1D) ENDIF ! For thermally-thick boundary conditions, set the flag, CALL_HT_1D, to call the subroutine SOLID_HEAT_TRANSFER. -! This routine is called every WALL_INCREMENT time steps. -! Also, if any 3D heat transfer calculations are done, set a new "sweep direction" (1, 2, or 3) +! This routine is called every WALL_INCREMENT time steps. Also, if any 3D heat transfer calculations are done, set a new "sweep +! direction" (1, 2, or 3) CALL_HT_1D = .FALSE. @@ -120,20 +83,46 @@ SUBROUTINE THERMAL_BC(T,NM,CALL_HT_1D) ENDIF ENDIF -! Loop through all wall cells and apply heat transfer BCs +! Sweep through all WALL cells and apply boundary conditions WALL_CELL_LOOP: DO IW=1,N_EXTERNAL_WALL_CELLS+N_INTERNAL_WALL_CELLS + WC=>WALL(IW) IF (WC%BOUNDARY_TYPE==NULL_BOUNDARY) CYCLE WALL_CELL_LOOP SF => SURFACE(WC%SURF_INDEX) BC => BOUNDARY_COORD(WC%BC_INDEX) B1 => BOUNDARY_PROP1(WC%B1_INDEX) + B2 => BOUNDARY_PROP2(WC%B2_INDEX) + CALL NEAR_SURFACE_GAS_VARIABLES(T,SF,BC,B1,WALL_INDEX=IW) + IF (.NOT.SF%THERMAL_BC_INDEX==THERMALLY_THICK) THEN CALL SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX=IW) ELSEIF (CALL_HT_1D) THEN CALL SOLID_HEAT_TRANSFER(NM,T,SF%HT_DIM*DT_BC,WALL_INDEX=IW) ENDIF + + IF (N_TRACKED_SPECIES>1 .AND. WC%BOUNDARY_TYPE/=OPEN_BOUNDARY .AND. WC%BOUNDARY_TYPE/=INTERPOLATED_BOUNDARY) THEN + CALL CALCULATE_RHO_D_F(B1,BC,WALL_INDEX=IW) + 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) SURFACE(CFA%SURF_INDEX) BC => BOUNDARY_COORD(CFA%BC_INDEX) B1 => BOUNDARY_PROP1(CFA%B1_INDEX) + B2 => BOUNDARY_PROP2(CFA%B2_INDEX) CALL CFACE_THERMAL_GASVARS(ICF,SF,B1,T) IF (.NOT.SF%THERMAL_BC_INDEX==THERMALLY_THICK) THEN CALL SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,CFACE_INDEX=ICF) ELSEIF (CALL_HT_1D) THEN CALL SOLID_HEAT_TRANSFER(NM,T,DT_BC,CFACE_INDEX=ICF) ENDIF + + IF (N_TRACKED_SPECIES>1) CALL CALCULATE_RHO_D_F(B1,BC,CFACE_INDEX=ICF) + + 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) LAGRANGIAN_PARTICLE(IP) LPC => LAGRANGIAN_PARTICLE_CLASS(LP%CLASS_INDEX) IF (.NOT.LPC%SOLID_PARTICLE .AND. .NOT.LPC%MASSLESS_TARGET) CYCLE PARTICLE_LOOP SF => SURFACE(LPC%SURF_INDEX) BC => BOUNDARY_COORD(LP%BC_INDEX) B1 => BOUNDARY_PROP1(LP%B1_INDEX) + CALL NEAR_SURFACE_GAS_VARIABLES(T,SF,BC,B1,LP=LP,PARTICLE_INDEX=IP) + IF (.NOT.SF%THERMAL_BC_INDEX==THERMALLY_THICK) THEN CALL SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,PARTICLE_INDEX=IP) ELSEIF (CALL_HT_1D) THEN CALL SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX=IP) ENDIF + + IF (LPC%SOLID_PARTICLE) THEN + CALL CALCULATE_ZZ_F(T,DT,PARTICLE_INDEX=IP) + IF (CORRECTOR) CALL DEPOSIT_PARTICLE_MASS(NM,IP,LP,LPC,CALL_HT_1D) ! Add the particle off-gas to the gas phase mesh + ENDIF + ENDDO PARTICLE_LOOP + ENDIF -END SUBROUTINE THERMAL_BC +T_USED(6)=T_USED(6)+CURRENT_TIME()-TNOW +END SUBROUTINE WALL_BC !> \brief Set various near-surface variables to be used for boundary conditions @@ -825,142 +845,12 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I END SUBROUTINE SURFACE_HEAT_TRANSFER -!> \brief Transfer heat from one 3-D sweep dirction to the other two. -!> \param NM Mesh index -!> \details Sweep over all 3D wall cells and look for those that are not oriented in the current sweep direction. -!> Copy the updated temperatures of the sweep direction to the non-sweep directions. - -SUBROUTINE HT3D_TEMPERATURE_EXCHANGE(NM) - -INTEGER, INTENT(IN) :: NM -INTEGER :: NWP,I_NODE,I,II,IWA,NM2,IW,ITW -REAL(EB) :: TMP_1,TMP_NWP -TYPE(BOUNDARY_ONE_D_TYPE), POINTER :: ONE_D2 -TYPE(BOUNDARY_THR_D_TYPE), POINTER :: THR_D -TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC,BC2 -TYPE(WALL_TYPE), POINTER :: WC,WC2 -TYPE(THIN_WALL_TYPE), POINTER :: TW,TW2 -TYPE(MESH_TYPE), POINTER :: M -TYPE(SURFACE_TYPE), POINTER :: SF -TYPE(BOUNDARY_ONE_D_TYPE), POINTER :: ONE_D - -M => MESHES(NM) - -WALL_LOOP: DO IW=1,M%N_EXTERNAL_WALL_CELLS+M%N_INTERNAL_WALL_CELLS - - WC => M%WALL(IW) - SF => SURFACE(WC%SURF_INDEX) - IF (SF%HT_DIM==1 .OR. WC%BOUNDARY_TYPE/=SOLID_BOUNDARY) CYCLE WALL_LOOP - BC => M%BOUNDARY_COORD(WC%BC_INDEX) - IF (ABS(BC%IOR)==M%HT_3D_SWEEP_DIRECTION) CYCLE WALL_LOOP - - ONE_D => M%BOUNDARY_ONE_D(WC%OD_INDEX) - NWP = SUM(ONE_D%N_LAYER_CELLS(1:ONE_D%N_LAYERS)) - THR_D => M%BOUNDARY_THR_D(WC%TD_INDEX) - - TMP_1 = ONE_D%TMP(1) - TMP_NWP = ONE_D%TMP(NWP) - NODE_LOOP: DO I=1,NWP ! Nodes of the chosen wall cell. - IF (.NOT.THR_D%NODE(I)%HT3D) CYCLE NODE_LOOP - IF (THR_D%NODE(I)%ALTERNATE_WALL_COUNT==0) CYCLE NODE_LOOP - IF (.NOT.ANY(ABS(THR_D%NODE(I)%ALTERNATE_WALL_IOR(:))==M%HT_3D_SWEEP_DIRECTION)) CYCLE NODE_LOOP - WEIGHT_LOOP: DO II=1,THR_D%NODE(I)%ALTERNATE_WALL_COUNT - IWA = THR_D%NODE(I)%ALTERNATE_WALL_INDEX(II) - NM2 = THR_D%NODE(I)%ALTERNATE_WALL_MESH(II) - I_NODE = THR_D%NODE(I)%ALTERNATE_WALL_NODE(II) - IF (THR_D%NODE(I)%ALTERNATE_WALL_TYPE(II)==0) THEN - WC2 => MESHES(NM2)%WALL(IWA) - BC2 => MESHES(NM2)%BOUNDARY_COORD(WC2%BC_INDEX) - ONE_D2 => MESHES(NM2)%BOUNDARY_ONE_D(WC2%OD_INDEX) - ELSE - TW2 => MESHES(NM2)%THIN_WALL(IWA) - BC2 => MESHES(NM2)%BOUNDARY_COORD(TW2%BC_INDEX) - ONE_D2 => MESHES(NM2)%BOUNDARY_ONE_D(TW2%OD_INDEX) - ENDIF - IF (ABS(BC2%IOR)/=M%HT_3D_SWEEP_DIRECTION) CYCLE WEIGHT_LOOP - ONE_D%TMP(I) = (ONE_D%RHO_C_S(I)*ONE_D%TMP(I) + & - THR_D%NODE(I)%ALTERNATE_WALL_WEIGHT(II)*ONE_D2%RHO_C_S(I_NODE)*ONE_D2%DELTA_TMP(I_NODE))/ONE_D%RHO_C_S(I) - - ENDDO WEIGHT_LOOP - IF (THR_D%NODE(I)%MESH_NUMBER==NM) M%TMP(THR_D%NODE(I)%I,THR_D%NODE(I)%J,THR_D%NODE(I)%K) = ONE_D%TMP(I) - ENDDO NODE_LOOP - ONE_D%TMP(0) = ONE_D%TMP(0) + ONE_D%TMP(1) - TMP_1 - ONE_D%TMP(NWP+1) = ONE_D%TMP(NWP+1) + ONE_D%TMP(NWP) - TMP_NWP - -ENDDO WALL_LOOP - -THIN_WALL_LOOP: DO ITW=1,M%N_THIN_WALL_CELLS - - TW => M%THIN_WALL(ITW) - SF => SURFACE(TW%SURF_INDEX) - BC => M%BOUNDARY_COORD(TW%BC_INDEX) - IF (ABS(BC%IOR)==M%HT_3D_SWEEP_DIRECTION) CYCLE THIN_WALL_LOOP - ONE_D => M%BOUNDARY_ONE_D(TW%OD_INDEX) - NWP = SUM(ONE_D%N_LAYER_CELLS(1:SF%N_LAYERS)) - THR_D => M%BOUNDARY_THR_D(TW%TD_INDEX) - IF (.NOT.ALLOCATED(THR_D%NODE)) CYCLE THIN_WALL_LOOP - - NODE_LOOP_2: DO I=1,NWP ! Nodes of the chosen wall cell. - IF (THR_D%NODE(I)%ALTERNATE_WALL_COUNT==0) CYCLE NODE_LOOP_2 - IF (.NOT.ANY(ABS(THR_D%NODE(I)%ALTERNATE_WALL_IOR(:))==M%HT_3D_SWEEP_DIRECTION)) CYCLE NODE_LOOP_2 - WEIGHT_LOOP_2: DO II=1,THR_D%NODE(I)%ALTERNATE_WALL_COUNT - IWA = THR_D%NODE(I)%ALTERNATE_WALL_INDEX(II) - NM2 = THR_D%NODE(I)%ALTERNATE_WALL_MESH(II) - I_NODE = THR_D%NODE(I)%ALTERNATE_WALL_NODE(II) - IF (THR_D%NODE(I)%ALTERNATE_WALL_TYPE(II)==0) THEN - WC2 => MESHES(NM2)%WALL(IWA) - BC2 => MESHES(NM2)%BOUNDARY_COORD(WC2%BC_INDEX) - ONE_D2 => MESHES(NM2)%BOUNDARY_ONE_D(WC2%OD_INDEX) - ELSE - TW2 => MESHES(NM2)%THIN_WALL(IWA) - BC2 => MESHES(NM2)%BOUNDARY_COORD(TW2%BC_INDEX) - ONE_D2 => MESHES(NM2)%BOUNDARY_ONE_D(TW2%OD_INDEX) - ENDIF - IF (ABS(BC2%IOR)/=M%HT_3D_SWEEP_DIRECTION) CYCLE WEIGHT_LOOP_2 - ONE_D%TMP(I) = ONE_D%TMP(I) + THR_D%NODE(I)%ALTERNATE_WALL_WEIGHT(II)*ONE_D2%DELTA_TMP(I_NODE) - ENDDO WEIGHT_LOOP_2 - ENDDO NODE_LOOP_2 - -ENDDO THIN_WALL_LOOP - -END SUBROUTINE HT3D_TEMPERATURE_EXCHANGE - - -!> \brief Calculate the term RHO_D_F=RHO*D at the wall. - -SUBROUTINE DIFFUSIVITY_BC - -INTEGER :: IW,ICF -TYPE(CFACE_TYPE), POINTER :: CFA -TYPE(WALL_TYPE), POINTER :: WC -TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1 -TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC - -IF (N_TRACKED_SPECIES==1) RETURN - -! Loop over all WALL cells - -WALL_LOOP: DO IW=1,N_EXTERNAL_WALL_CELLS+N_INTERNAL_WALL_CELLS - WC=>WALL(IW) - IF (WC%BOUNDARY_TYPE==NULL_BOUNDARY .OR. & - WC%BOUNDARY_TYPE==OPEN_BOUNDARY .OR. & - WC%BOUNDARY_TYPE==INTERPOLATED_BOUNDARY) CYCLE WALL_LOOP - B1 => BOUNDARY_PROP1(WC%B1_INDEX) - BC => BOUNDARY_COORD(WC%BC_INDEX) - CALL CALCULATE_RHO_D_F(B1,BC,WALL_INDEX=IW) -ENDDO WALL_LOOP - -! Loop over all cut face cells - -CFACE_LOOP: DO ICF=INTERNAL_CFACE_CELLS_LB+1,INTERNAL_CFACE_CELLS_LB+N_INTERNAL_CFACE_CELLS - CFA => CFACE(ICF) - B1 => BOUNDARY_PROP1(CFA%B1_INDEX) - BC => BOUNDARY_COORD(CFA%BC_INDEX) - CALL CALCULATE_RHO_D_F(B1,BC,CFACE_INDEX=ICF) -ENDDO CFACE_LOOP - -END SUBROUTINE DIFFUSIVITY_BC - +!> \brief Calculate the diffusion coefficient at the boundary +!> \param B1 Pointer to boundary properties +!> \param BC Pointer to boundary coordinates +!> \param WALL_INDEX Optional WALL cell index +!> \param CFACE_INDEX Optional immersed boundary (CFACE) index +!> \details Calculate the diffusion coefficient, rho*D (kg/m/s) SUBROUTINE CALCULATE_RHO_D_F(B1,BC,WALL_INDEX,CFACE_INDEX) @@ -1000,52 +890,13 @@ SUBROUTINE CALCULATE_RHO_D_F(B1,BC,WALL_INDEX,CFACE_INDEX) END SUBROUTINE CALCULATE_RHO_D_F -!> \brief Compute the species mass fractions at the boundary, ZZ_F. -!> +!> \brief Calculate the species mass fractions, ZZ, at the boundary !> \param T Current time (s) !> \param DT Current time step (s) -!> \param NM Mesh number -!> \param CALL_HT_1D Indicates whether 1-D conduction solver is called - -SUBROUTINE SPECIES_BC(T,DT,NM,CALL_HT_1D) - -REAL(EB), INTENT(IN) :: T,DT -INTEGER, INTENT(IN) :: NM -LOGICAL, INTENT(IN) :: CALL_HT_1D -INTEGER :: IW,ICF,IP -TYPE(WALL_TYPE), POINTER :: WC -TYPE(LAGRANGIAN_PARTICLE_TYPE), POINTER :: LP -TYPE(LAGRANGIAN_PARTICLE_CLASS_TYPE), POINTER :: LPC - -! Add evaporating gases from solid particles to the mesh using a volumetric source term - -PARTICLE_LOOP: DO IP=1,NLP - LP => LAGRANGIAN_PARTICLE(IP) - LPC => LAGRANGIAN_PARTICLE_CLASS(LP%CLASS_INDEX) - IF (.NOT.LPC%SOLID_PARTICLE) CYCLE PARTICLE_LOOP - CALL CALCULATE_ZZ_F(T,DT,PARTICLE_INDEX=IP) - IF (PREDICTOR) CYCLE PARTICLE_LOOP ! Only do basic boundary conditions during the PREDICTOR stage of time step. - CALL DEPOSIT_PARTICLE_MASS(NM,IP,LP,LPC,CALL_HT_1D) ! Add the particle off-gas to the gas phase mesh -ENDDO PARTICLE_LOOP - -! Loop through the wall cells, apply mass boundary conditions - -WALL_CELL_LOOP: DO IW=1,N_EXTERNAL_WALL_CELLS+N_INTERNAL_WALL_CELLS - WC => WALL(IW) - IF (WC%BOUNDARY_TYPE==OPEN_BOUNDARY) CYCLE WALL_CELL_LOOP - IF (WC%BOUNDARY_TYPE==NULL_BOUNDARY) CYCLE WALL_CELL_LOOP - IF (WC%BOUNDARY_TYPE==INTERPOLATED_BOUNDARY) CYCLE WALL_CELL_LOOP - CALL CALCULATE_ZZ_F(T,DT,WALL_INDEX=IW) -ENDDO WALL_CELL_LOOP - -! Loop through the cut face cells, apply mass boundary conditions - -DO ICF=INTERNAL_CFACE_CELLS_LB+1,INTERNAL_CFACE_CELLS_LB+N_INTERNAL_CFACE_CELLS - CALL CALCULATE_ZZ_F(T,DT,CFACE_INDEX=ICF) -ENDDO - -END SUBROUTINE SPECIES_BC - +!> \param WALL_INDEX Optional WALL cell index +!> \param CFACE_INDEX Optional immersed boundary (CFACE) index +!> \param PARTICLE_INDEX Optional particle index +!> \details Calculate the diffusion coefficient, rho*D (kg/m/s) SUBROUTINE CALCULATE_ZZ_F(T,DT,WALL_INDEX,CFACE_INDEX,PARTICLE_INDEX) @@ -1445,6 +1296,12 @@ END SUBROUTINE CALCULATE_ZZ_F !> \brief Deposit particle off-gas onto mesh +!> \param NM Mesh number +!> \param IP Particle index +!> \param LP Pointer to Lagrangian Particle derived type variable +!> \param LPC Pointer to Lagrangian Particle Class +!> \param CALL_HT_1D Logical indicating if the 1-D heat transfer routine was called +!> \details Deposit the particle off-gas onto the mesh SUBROUTINE DEPOSIT_PARTICLE_MASS(NM,IP,LP,LPC,CALL_HT_1D) @@ -1571,38 +1428,6 @@ SUBROUTINE DEPOSIT_PARTICLE_MASS(NM,IP,LP,LPC,CALL_HT_1D) END SUBROUTINE DEPOSIT_PARTICLE_MASS -!> \brief Compute density at wall from wall temperatures and mass fractions - -SUBROUTINE DENSITY_BC - -INTEGER :: IW,ICF -TYPE(WALL_TYPE), POINTER :: WC -TYPE(CFACE_TYPE), POINTER :: CFA -TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1 -TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC - -! Loop over all wall cells - -WALL_CELL_LOOP: DO IW=1,N_EXTERNAL_WALL_CELLS+N_INTERNAL_WALL_CELLS - WC => WALL(IW) - IF (WC%BOUNDARY_TYPE==NULL_BOUNDARY .OR. WC%BOUNDARY_TYPE==INTERPOLATED_BOUNDARY) CYCLE WALL_CELL_LOOP - B1 => BOUNDARY_PROP1(WC%B1_INDEX) - BC => BOUNDARY_COORD(WC%BC_INDEX) - CALL CALCULATE_RHO_F(BC,B1,WALL_INDEX=IW) -ENDDO WALL_CELL_LOOP - -! Loop over all cut face cells - -CFACE_LOOP: DO ICF=INTERNAL_CFACE_CELLS_LB+1,INTERNAL_CFACE_CELLS_LB+N_INTERNAL_CFACE_CELLS - CFA => CFACE(ICF) - B1 => BOUNDARY_PROP1(CFA%B1_INDEX) - BC => BOUNDARY_COORD(CFA%BC_INDEX) - CALL CALCULATE_RHO_F(BC,B1,CFACE_INDEX=ICF) -ENDDO CFACE_LOOP - -END SUBROUTINE DENSITY_BC - - !> \brief Compute density, RHO_F, at non-iterpolated boundaries !> \param BC Boundary Coordinates derived type !> \param B1 Boundary Properties derived type @@ -1649,44 +1474,13 @@ SUBROUTINE CALCULATE_RHO_F(BC,B1,WALL_INDEX,CFACE_INDEX) END SUBROUTINE CALCULATE_RHO_F -SUBROUTINE DEPOSITION_BC(DT) - -REAL(EB), INTENT(IN) :: DT -INTEGER:: IW,ICF -TYPE(SURFACE_TYPE), POINTER :: SF -TYPE(WALL_TYPE), POINTER :: WC -TYPE(CFACE_TYPE), POINTER :: CFA -TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC -TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1 -TYPE(BOUNDARY_PROP2_TYPE), POINTER :: B2 - -WALL_CELL_LOOP: DO IW=1,N_EXTERNAL_WALL_CELLS+N_INTERNAL_WALL_CELLS - WC=>WALL(IW) - B1=>BOUNDARY_PROP1(WC%B1_INDEX) - ! No deposition if the boundary isn't solid or has a specified flow - IF (WC%BOUNDARY_TYPE/=SOLID_BOUNDARY .OR. B1%NODE_INDEX /= 0) CYCLE WALL_CELL_LOOP - SF=>SURFACE(WC%SURF_INDEX) - IF (ABS(SF%VEL)>TWO_EPSILON_EB .OR. ANY(ABS(SF%MASS_FLUX)>TWO_EPSILON_EB) .OR. ABS(SF%VOLUME_FLOW)>TWO_EPSILON_EB) & - CYCLE WALL_CELL_LOOP - BC=>BOUNDARY_COORD(WC%BC_INDEX) - B2=>BOUNDARY_PROP2(WC%B2_INDEX) - CALL CALC_DEPOSITION(DT,BC,B1,B2,WALL_INDEX=IW) -ENDDO WALL_CELL_LOOP - -CFACE_LOOP: DO ICF=INTERNAL_CFACE_CELLS_LB+1,INTERNAL_CFACE_CELLS_LB+N_INTERNAL_CFACE_CELLS - CFA=>CFACE(ICF) - IF (CFA%BOUNDARY_TYPE/=SOLID_BOUNDARY) CYCLE CFACE_LOOP - SF=>SURFACE(CFA%SURF_INDEX) - IF (ABS(SF%VEL)>TWO_EPSILON_EB .OR. ANY(ABS(SF%MASS_FLUX)>TWO_EPSILON_EB) .OR. ABS(SF%VOLUME_FLOW)>TWO_EPSILON_EB) & - CYCLE CFACE_LOOP - BC=>BOUNDARY_COORD(CFA%BC_INDEX) - B1=>BOUNDARY_PROP1(CFA%B1_INDEX) - B2=>BOUNDARY_PROP2(CFA%B2_INDEX) - CALL CALC_DEPOSITION(DT,BC,B1,B2,CFACE_INDEX=ICF) -ENDDO CFACE_LOOP - -END SUBROUTINE DEPOSITION_BC - +!> \brief Calculate aerosol deposition onto a solid surface +!> \param DT Current time step +!> \param BC Pointer to Boundary Coordinate derived type variable +!> \param B1 Pointer to Boundary Property derived type variable +!> \param B2 Pointer to Boundary Property derived type variable +!> \param WALL_INDEX Optional WALL index +!> \param CFACE_INDEX Optional CFACE index SUBROUTINE CALC_DEPOSITION(DT,BC,B1,B2,WALL_INDEX,CFACE_INDEX) @@ -1774,39 +1568,10 @@ SUBROUTINE CALC_DEPOSITION(DT,BC,B1,B2,WALL_INDEX,CFACE_INDEX) END SUBROUTINE CALC_DEPOSITION -!> \brief Compute density at wall from wall temperatures and mass fractions - -SUBROUTINE HVAC_BC - -INTEGER :: IW,ICF -TYPE(CFACE_TYPE), POINTER :: CFA -TYPE(WALL_TYPE), POINTER :: WC -TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1 -TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC -TYPE(SURFACE_TYPE), POINTER :: SF - -! Loop over all internal and external wall cells - -WALL_CELL_LOOP: DO IW=1,N_EXTERNAL_WALL_CELLS+N_INTERNAL_WALL_CELLS - WC => WALL(IW) - B1 => BOUNDARY_PROP1(WC%B1_INDEX) - IF (B1%NODE_INDEX == 0) CYCLE WALL_CELL_LOOP - BC => BOUNDARY_COORD(WC%BC_INDEX) - SF => SURFACE(WC%SURF_INDEX) - CALL CALC_HVAC_BC(BC,B1,SF) -ENDDO WALL_CELL_LOOP - -CFACE_LOOP: DO ICF=INTERNAL_CFACE_CELLS_LB+1,INTERNAL_CFACE_CELLS_LB+N_INTERNAL_CFACE_CELLS - CFA => CFACE(ICF) - B1 => BOUNDARY_PROP1(CFA%B1_INDEX) - IF (B1%NODE_INDEX == 0) CYCLE CFACE_LOOP - BC => BOUNDARY_COORD(CFA%BC_INDEX) - SF => SURFACE(CFA%SURF_INDEX) - CALL CALC_HVAC_BC(BC,B1,SF) -ENDDO CFACE_LOOP - -END SUBROUTINE HVAC_BC - +!> \brief Calculate HVAC boundary conditions +!> \param BC Pointer to Boundary Coordinate derived type variable +!> \param B1 Pointer to Boundary Property derived type variable +!> \param SF Pointer to Surface Propertes derived type variable SUBROUTINE CALC_HVAC_BC(BC,B1,SF) @@ -3631,6 +3396,111 @@ REAL(EB) FUNCTION HEAT_TRANSFER_COEFFICIENT(NMX,DELTA_N_TMP,H_FIXED,SFX,WALL_IND END FUNCTION HEAT_TRANSFER_COEFFICIENT +!> \brief Transfer heat from one 3-D sweep dirction to the other two. +!> \param NM Mesh index +!> \details Sweep over all 3D wall cells and look for those that are not oriented in the current sweep direction. +!> Copy the updated temperatures of the sweep direction to the non-sweep directions. + +SUBROUTINE HT3D_TEMPERATURE_EXCHANGE(NM) + +USE COMP_FUNCTIONS, ONLY: CURRENT_TIME +INTEGER, INTENT(IN) :: NM +INTEGER :: NWP,I_NODE,I,II,IWA,NM2,IW,ITW +REAL(EB) :: TMP_1,TMP_NWP,TNOW +TYPE(BOUNDARY_ONE_D_TYPE), POINTER :: ONE_D2 +TYPE(BOUNDARY_THR_D_TYPE), POINTER :: THR_D +TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC,BC2 +TYPE(WALL_TYPE), POINTER :: WC,WC2 +TYPE(THIN_WALL_TYPE), POINTER :: TW,TW2 +TYPE(MESH_TYPE), POINTER :: M +TYPE(SURFACE_TYPE), POINTER :: SF +TYPE(BOUNDARY_ONE_D_TYPE), POINTER :: ONE_D + +TNOW = CURRENT_TIME() + +M => MESHES(NM) + +WALL_LOOP: DO IW=1,M%N_EXTERNAL_WALL_CELLS+M%N_INTERNAL_WALL_CELLS + + WC => M%WALL(IW) + SF => SURFACE(WC%SURF_INDEX) + IF (SF%HT_DIM==1 .OR. WC%BOUNDARY_TYPE/=SOLID_BOUNDARY) CYCLE WALL_LOOP + BC => M%BOUNDARY_COORD(WC%BC_INDEX) + IF (ABS(BC%IOR)==M%HT_3D_SWEEP_DIRECTION) CYCLE WALL_LOOP + + ONE_D => M%BOUNDARY_ONE_D(WC%OD_INDEX) + NWP = SUM(ONE_D%N_LAYER_CELLS(1:ONE_D%N_LAYERS)) + THR_D => M%BOUNDARY_THR_D(WC%TD_INDEX) + + TMP_1 = ONE_D%TMP(1) + TMP_NWP = ONE_D%TMP(NWP) + NODE_LOOP: DO I=1,NWP ! Nodes of the chosen wall cell. + IF (.NOT.THR_D%NODE(I)%HT3D) CYCLE NODE_LOOP + IF (THR_D%NODE(I)%ALTERNATE_WALL_COUNT==0) CYCLE NODE_LOOP + IF (.NOT.ANY(ABS(THR_D%NODE(I)%ALTERNATE_WALL_IOR(:))==M%HT_3D_SWEEP_DIRECTION)) CYCLE NODE_LOOP + WEIGHT_LOOP: DO II=1,THR_D%NODE(I)%ALTERNATE_WALL_COUNT + IWA = THR_D%NODE(I)%ALTERNATE_WALL_INDEX(II) + NM2 = THR_D%NODE(I)%ALTERNATE_WALL_MESH(II) + I_NODE = THR_D%NODE(I)%ALTERNATE_WALL_NODE(II) + IF (THR_D%NODE(I)%ALTERNATE_WALL_TYPE(II)==0) THEN + WC2 => MESHES(NM2)%WALL(IWA) + BC2 => MESHES(NM2)%BOUNDARY_COORD(WC2%BC_INDEX) + ONE_D2 => MESHES(NM2)%BOUNDARY_ONE_D(WC2%OD_INDEX) + ELSE + TW2 => MESHES(NM2)%THIN_WALL(IWA) + BC2 => MESHES(NM2)%BOUNDARY_COORD(TW2%BC_INDEX) + ONE_D2 => MESHES(NM2)%BOUNDARY_ONE_D(TW2%OD_INDEX) + ENDIF + IF (ABS(BC2%IOR)/=M%HT_3D_SWEEP_DIRECTION) CYCLE WEIGHT_LOOP + ONE_D%TMP(I) = (ONE_D%RHO_C_S(I)*ONE_D%TMP(I) + & + THR_D%NODE(I)%ALTERNATE_WALL_WEIGHT(II)*ONE_D2%RHO_C_S(I_NODE)*ONE_D2%DELTA_TMP(I_NODE))/ONE_D%RHO_C_S(I) + + ENDDO WEIGHT_LOOP + IF (THR_D%NODE(I)%MESH_NUMBER==NM) M%TMP(THR_D%NODE(I)%I,THR_D%NODE(I)%J,THR_D%NODE(I)%K) = ONE_D%TMP(I) + ENDDO NODE_LOOP + ONE_D%TMP(0) = ONE_D%TMP(0) + ONE_D%TMP(1) - TMP_1 + ONE_D%TMP(NWP+1) = ONE_D%TMP(NWP+1) + ONE_D%TMP(NWP) - TMP_NWP + +ENDDO WALL_LOOP + +THIN_WALL_LOOP: DO ITW=1,M%N_THIN_WALL_CELLS + + TW => M%THIN_WALL(ITW) + SF => SURFACE(TW%SURF_INDEX) + BC => M%BOUNDARY_COORD(TW%BC_INDEX) + IF (ABS(BC%IOR)==M%HT_3D_SWEEP_DIRECTION) CYCLE THIN_WALL_LOOP + ONE_D => M%BOUNDARY_ONE_D(TW%OD_INDEX) + NWP = SUM(ONE_D%N_LAYER_CELLS(1:SF%N_LAYERS)) + THR_D => M%BOUNDARY_THR_D(TW%TD_INDEX) + IF (.NOT.ALLOCATED(THR_D%NODE)) CYCLE THIN_WALL_LOOP + + NODE_LOOP_2: DO I=1,NWP ! Nodes of the chosen wall cell. + IF (THR_D%NODE(I)%ALTERNATE_WALL_COUNT==0) CYCLE NODE_LOOP_2 + IF (.NOT.ANY(ABS(THR_D%NODE(I)%ALTERNATE_WALL_IOR(:))==M%HT_3D_SWEEP_DIRECTION)) CYCLE NODE_LOOP_2 + WEIGHT_LOOP_2: DO II=1,THR_D%NODE(I)%ALTERNATE_WALL_COUNT + IWA = THR_D%NODE(I)%ALTERNATE_WALL_INDEX(II) + NM2 = THR_D%NODE(I)%ALTERNATE_WALL_MESH(II) + I_NODE = THR_D%NODE(I)%ALTERNATE_WALL_NODE(II) + IF (THR_D%NODE(I)%ALTERNATE_WALL_TYPE(II)==0) THEN + WC2 => MESHES(NM2)%WALL(IWA) + BC2 => MESHES(NM2)%BOUNDARY_COORD(WC2%BC_INDEX) + ONE_D2 => MESHES(NM2)%BOUNDARY_ONE_D(WC2%OD_INDEX) + ELSE + TW2 => MESHES(NM2)%THIN_WALL(IWA) + BC2 => MESHES(NM2)%BOUNDARY_COORD(TW2%BC_INDEX) + ONE_D2 => MESHES(NM2)%BOUNDARY_ONE_D(TW2%OD_INDEX) + ENDIF + IF (ABS(BC2%IOR)/=M%HT_3D_SWEEP_DIRECTION) CYCLE WEIGHT_LOOP_2 + ONE_D%TMP(I) = ONE_D%TMP(I) + THR_D%NODE(I)%ALTERNATE_WALL_WEIGHT(II)*ONE_D2%DELTA_TMP(I_NODE) + ENDDO WEIGHT_LOOP_2 + ENDDO NODE_LOOP_2 + +ENDDO THIN_WALL_LOOP + +T_USED(6)=T_USED(6)+CURRENT_TIME()-TNOW +END SUBROUTINE HT3D_TEMPERATURE_EXCHANGE + + !> \brief Perform a numerical TGA (thermo-gravimetric analysis) at the start of the simulation SUBROUTINE TGA_ANALYSIS