diff --git a/Build/makefile b/Build/makefile index 927807b34c8..9d3b6e6f515 100644 --- a/Build/makefile +++ b/Build/makefile @@ -59,6 +59,8 @@ FFLAGSMKL_CUSTOM = LFLAGSMKL_CUSTOM = FFLAGS_SUNDIALS = LFLAGS_SUNDIALS = +FFLAGS_HYPRE = +LFLAGS_HYPRE = ifdef MKLROOT # This assumes the MKL library is available. ifeq ($(shell echo "check_quotes"),"check_quotes") # windows @@ -109,6 +111,11 @@ ifdef SUNDIALS_HOME # This assumes the SUNDIALS library is available. LFLAGS_SUNDIALS_WIN = ${SUNDIALS_HOME}/lib/sundials_fcvode_mod.lib ${SUNDIALS_HOME}/lib/sundials_fnvecserial_mod.lib ${SUNDIALS_HOME}/lib/sundials_cvode.lib /link /NODEFAULTLIB:MSVCRTD /NODEFAULTLIB:libcmtd.lib endif +ifdef HYPRE_HOME # This assumes the HYPRE library is available. + FFLAGS_HYPRE = -DWITH_HYPRE -I${HYPRE_HOME}/include + LFLAGS_HYPRE = -L${HYPRE_HOME}/lib -lHYPRE -lm +endif + obj_mpi = prec.o cons.o chem.o prop.o devc.o type.o data.o mesh.o func.o gsmv.o smvv.o rcal.o turb.o soot.o \ pois.o geom.o ccib.o radi.o part.o vege.o ctrl.o hvac.o mass.o \ wall.o fire.o velo.o pres.o init.o dump.o read.o divg.o main.o @@ -337,6 +344,14 @@ ompi_gnu_linux_db : obj = fds_ompi_gnu_linux_db ompi_gnu_linux_db : setup $(obj_mpi) $(FCOMPL) $(FFLAGS) $(FOPENMPFLAGS) -o $(obj) $(obj_mpi) $(LFLAGSMKL) +ompi_gnu_linux_dv : FFLAGS = -m64 -O1 -fbacktrace -std=f2018 -frecursive -ffpe-summary=none -fall-intrinsics $(GITINFOGNU) $(FFLAGSMKL_GNU_OPENMPI) $(GFORTRAN_OPTIONS) +ompi_gnu_linux_dv : LFLAGSMKL = $(LFLAGSMKL_GNU_OPENMPI) +ompi_gnu_linux_dv : FCOMPL = mpifort +ompi_gnu_linux_dv : FOPENMPFLAGS = -fopenmp +ompi_gnu_linux_dv : obj = fds_ompi_gnu_linux_dv +ompi_gnu_linux_dv : setup $(obj_mpi) + $(FCOMPL) $(FFLAGS) $(FOPENMPFLAGS) -o $(obj) $(obj_mpi) $(LFLAGSMKL) + ompi_gnu_osx : FFLAGS = -m64 -O2 -std=f2018 -frecursive -ffpe-summary=none -fall-intrinsics $(GITINFOGNU) $(FFLAGSMKL_GNU_CUSTOM) $(GFORTRAN_OPTIONS) $(FFLAGS_SUNDIALS) ompi_gnu_osx : LFLAGSMKL = $(LFLAGSMKL_GNU_CUSTOM) $(CLT_VERSION) $(LFLAGS_SUNDIALS_OSX) ompi_gnu_osx : FCOMPL = mpifort diff --git a/Build/ompi_gnu_linux_dv/make_fds.sh b/Build/ompi_gnu_linux_dv/make_fds.sh new file mode 100755 index 00000000000..885bf58554f --- /dev/null +++ b/Build/ompi_gnu_linux_dv/make_fds.sh @@ -0,0 +1,6 @@ +#!/bin/bash +dir=`pwd` +target=${dir##*/} + +echo Building $target +make -j4 VPATH="../../Source" -f ../makefile $target diff --git a/Manuals/FDS_User_Guide/FDS_User_Guide.tex b/Manuals/FDS_User_Guide/FDS_User_Guide.tex index 2641f978e7b..1844760a055 100644 --- a/Manuals/FDS_User_Guide/FDS_User_Guide.tex +++ b/Manuals/FDS_User_Guide/FDS_User_Guide.tex @@ -3066,11 +3066,12 @@ \subsection{Liquid Fuels} \subsubsection{Evaporation of a Pure Liquid} \label{methanol_evaporation} -An example of liquid evaporation is given by the sample case found in the {\ct Pyrolysis} folder called {\ct methanol\_evaporation.fds}. A 1~m by 1~m pan filled with methanol at $T_\infty=20$~$^\circ$C is exposed to a uniform heat flux, $\dot{q}''=20$~\unit{kW/m^2}. The boiling temperature of methanol is $T_{\rm b}=64.65$~$^\circ$C, its specific heat, $c=2.48$~kJ/(kg$\cdot$K), and heat of vaporization, $h_{\rm v}=1099$~kJ/kg. The evaporation rate of a burning liquid in steady state is approximately +An example of liquid evaporation is given by the sample case found in the {\ct Pyrolysis} folder called {\ct methanol\_evaporation.fds}. A 1~m by 1~m pan filled with methanol at $T_\infty=20$~$^\circ$C is exposed to a uniform heat flux, $\dot{q}''=20$~\unit{kW/m^2}. The boiling temperature of methanol is $T_{\rm b}=64.65$~$^\circ$C, its specific heat, $c=2.48$~kJ/(kg$\cdot$K), and heat of vaporization, $h_{\rm v}=1099$~kJ/kg. At steady state, the heat balance at the pool surface is \be - \dot{m}'' \approx \frac{\dot{q}''}{h_{\rm g}} \quad ; \quad h_{\rm g} = c (T_{\rm b}-T_\infty) + h_{\rm v} + \dot{q}''_{total} - \dot{q}''_{c}= \dot{m}'' h_{\rm v}(T_s) \ee -In this example, the methanol evaporates in an oxygen-depleted atmosphere and no burning occurs. The left hand plot in Fig.~\ref{methanol_evaporation_plot} displays the computed evaporation rate, $\dot{m}''$, versus the ideal, $\dot{q}''/h_{\rm g}$. The former approaches the latter as all of the absorbed energy is used to evaporate the liquid. The right hand plot shows the computed liquid surface temperature versus the liquid boiling temperature. + +where $ \dot{q}''_{c}$ is the heat being conducted away from the surface. If $ \dot{q}''_{c}$ is made zero, which can be done by using {\ct BACKING='INSULATED'} and a high thermal conductivity, then the pool surface temperature,$T_s$, will approach the boiling temperature,$T_b$. In this example, the methanol evaporates in an oxygen-depleted atmosphere and no burning occurs. The left hand plot in Fig.~\ref{methanol_evaporation_plot} displays the computed evaporation rate, $\dot{m}''$, versus the ideal, $\dot{q}''_{}total}/h_{\rm v}(T_b)$. The former approaches the latter as all of the absorbed energy is used to evaporate the liquid. The right hand plot shows the computed liquid surface temperature versus the liquid boiling temperature. \begin{figure}[!ht] \includegraphics[width=3.2in]{SCRIPT_FIGURES/methanol_evaporation_mdot} \includegraphics[width=3.2in]{SCRIPT_FIGURES/methanol_evaporation_temp} diff --git a/Source/cons.f90 b/Source/cons.f90 index 1b62ffe18e0..9bb91be4862 100644 --- a/Source/cons.f90 +++ b/Source/cons.f90 @@ -211,6 +211,7 @@ MODULE GLOBAL_CONSTANTS LOGICAL :: CHECK_VN=.TRUE. !< Check the Von Neumann number LOGICAL :: CHECK_FO=.FALSE. !< Check the solid phase Fourier number LOGICAL :: SOLID_PARTICLES=.FALSE. !< Indicates the existence of solid particles +LOGICAL :: ORIENTED_PARTICLES=.FALSE. !< Indicates the existence of particles with a specified orientation LOGICAL :: HVAC=.FALSE. !< Perform an HVAC calculation LOGICAL :: BAROCLINIC=.TRUE. !< Include the baroclinic terms in the momentum equation LOGICAL :: GRAVITATIONAL_DEPOSITION=.TRUE. !< Allow aerosol gravitational deposition @@ -248,6 +249,7 @@ MODULE GLOBAL_CONSTANTS LOGICAL :: DUCT_HT_INSERTED=.FALSE. LOGICAL :: HVAC_QFAN=.FALSE. LOGICAL :: USE_ATMOSPHERIC_INTERPOLATION=.FALSE. +LOGICAL :: NEAR_WALL_PARTICLE_INTERPOLATION=.FALSE. LOGICAL :: POSITIVE_ERROR_TEST=.FALSE. LOGICAL :: OBST_SHAPE_AREA_ADJUST=.FALSE. LOGICAL :: STORE_SPECIES_FLUX=.FALSE. @@ -725,8 +727,8 @@ MODULE GLOBAL_CONSTANTS REAL(EB), POINTER, DIMENSION(:,:) :: ORIENTATION_VECTOR !< Global array of orientation vectors INTEGER, ALLOCATABLE, DIMENSION(:) :: NEAREST_RADIATION_ANGLE !< Index of the rad angle most opposite the given ORIENTATION_VECTOR -REAL(EB), POINTER, DIMENSION(:) :: ORIENTATION_VIEW_ANGLE !< View angle of the given ORIENTATION_VECTOR -REAL(EB), ALLOCATABLE, DIMENSION(:) :: VIEW_ANGLE_AREA !< View angle area ORIENTATION_VECTOR +REAL(EB), POINTER, DIMENSION(:) :: COS_HALF_VIEW_ANGLE !< View angle of the given ORIENTATION_VECTOR +REAL(EB), ALLOCATABLE, DIMENSION(:) :: VIEW_ANGLE_FACTOR !< View angle area ORIENTATION_VECTOR INTEGER :: N_ORIENTATION_VECTOR !< Number of ORIENTATION_VECTORs INTEGER :: TGA_MESH_INDEX=HUGE(INTEGER_ONE) !< Mesh for the special TGA calculation diff --git a/Source/dump.f90 b/Source/dump.f90 index a9eb3d58965..ab1a03c0999 100644 --- a/Source/dump.f90 +++ b/Source/dump.f90 @@ -3011,7 +3011,10 @@ SUBROUTINE INITIALIZE_DIAGNOSTIC_FILE(DT) WRITE(LU_OUTPUT,'(A,A,A,F8.2)')' ',SPECIES_MIXTURE(NS)%ID,': ',ML%NU_GAS(NS,1) ENDDO WRITE(LU_OUTPUT,'(A,F8.2)') ' Boiling temperature (C): ',ML%TMP_BOIL-TMPM - WRITE(LU_OUTPUT,'(A,ES10.3)')' H_R (kJ/kg) : ',ML%H_R(1,NINT(TMPA))/1000._EB + ITMP = NINT(TMPA) + WRITE(LU_OUTPUT,'(A,I4,A,ES10.3)') ' H_R (kJ/kg) TMPA, ',ITMP,' K: ',ML%H_R(1,ITMP)/1000._EB + ITMP = NINT(ML%TMP_REF(1)) + WRITE(LU_OUTPUT,'(A,I4,A,ES10.3)') ' H_R (kJ/kg) TMP_REF, ',ITMP,' K: ',ML%H_R(1,ITMP)/1000._EB ENDIF ENDDO MATL_LOOP diff --git a/Source/func.f90 b/Source/func.f90 index 1a779b2403a..3a1d783b50f 100644 --- a/Source/func.f90 +++ b/Source/func.f90 @@ -1477,7 +1477,6 @@ SUBROUTINE PACK_PARTICLE(NM,OS,LP,LPC_INDEX,RC,IC,LC,UNPACK_IT,COUNT_ONLY) IC=IC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%INTEGERS(IC),LP%TAG,UNPACK_IT) IC=IC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%INTEGERS(IC),LP%CLASS_INDEX,UNPACK_IT) -IC=IC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%INTEGERS(IC),LP%INITIALIZATION_INDEX,UNPACK_IT) IC=IC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%INTEGERS(IC),LP%ORIENTATION_INDEX,UNPACK_IT) IC=IC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%INTEGERS(IC),LP%WALL_INDEX,UNPACK_IT) IC=IC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%INTEGERS(IC),LP%DUCT_INDEX,UNPACK_IT) @@ -3497,100 +3496,6 @@ SUBROUTINE UPDATE_HISTOGRAM(NBINS,LIMITS,COUNTS,VAL,WEIGHT) COUNTS(IND)=COUNTS(IND)+WEIGHT END SUBROUTINE UPDATE_HISTOGRAM - -!> \brief Linearly interpolate the a mesh quantity onto a point -!> \param X The interpolated value of the 3D array -!> \param A The 3D array of values -!> \param I The lower x index of the array -!> \param J The lower y index of the array -!> \param K The lower z index of the array -!> \param P Fraction of the distance from the lower to upper x coordinate -!> \param R Fraction of the distance from the lower to upper y coordinate -!> \param S Fraction of the distance from the lower to upper z coordinate - -SUBROUTINE MESH_TO_PARTICLE(X,A,I,J,K,P,R,S) - -REAL(EB), INTENT(IN), DIMENSION(0:,0:,0:) :: A -INTEGER, INTENT(IN) :: I,J,K -REAL(EB), INTENT(IN) :: P,R,S -REAL(EB), INTENT(OUT) :: X -REAL(EB) :: PP,RR,SS - -PP = 1._EB-P -RR = 1._EB-R -SS = 1._EB-S -X = ((PP*A(I,J,K) +P*A(I+1,J,K) )*RR+(PP*A(I,J+1,K) +P*A(I+1,J+1,K) )*R)*SS + & - ((PP*A(I,J,K+1)+P*A(I+1,J,K+1))*RR+(PP*A(I,J+1,K+1)+P*A(I+1,J+1,K+1))*R)*S - -END SUBROUTINE MESH_TO_PARTICLE - - -!> \brief Linearly interpolate the value at a point onto the mesh -!> \param X The interpolated value of the 3D array -!> \param A The 3D array of values -!> \param I The lower x index of the array -!> \param J The lower y index of the array -!> \param K The lower z index of the array -!> \param P Fraction of the distance from the lower to upper x coordinate -!> \param R Fraction of the distance from the lower to upper y coordinate -!> \param S Fraction of the distance from the lower to upper z coordinate - -SUBROUTINE PARTICLE_TO_MESH(X,A,I,J,K,P,R,S) - -REAL(EB), INTENT(INOUT), DIMENSION(0:,0:,0:) :: A -INTEGER, INTENT(IN) :: I,J,K -REAL(EB), INTENT(IN) :: P,R,S -REAL(EB), INTENT(IN) :: X -REAL(EB) :: PP,RR,SS - -PP = 1._EB-P -RR = 1._EB-R -SS = 1._EB-S -A(I ,J ,K ) = A(I ,J ,K ) - X*PP*RR*SS -A(I+1,J ,K ) = A(I+1,J ,K ) - X*P *RR*SS -A(I ,J+1,K ) = A(I ,J+1,K ) - X*PP*R *SS -A(I+1,J+1,K ) = A(I+1,J+1,K ) - X*P *R *SS -A(I ,J ,K+1) = A(I ,J ,K+1) - X*PP*RR*S -A(I+1,J ,K+1) = A(I+1,J ,K+1) - X*P *RR*S -A(I ,J+1,K+1) = A(I ,J+1,K+1) - X*PP*R *S -A(I+1,J+1,K+1) = A(I+1,J+1,K+1) - X*P *R *S - -END SUBROUTINE PARTICLE_TO_MESH - - -!> \brief Trilinear interpolation https://paulbourke.net/miscellaneous/interpolation/ -!> \param V The interpolated value of the 3D array -!> \param A The 3D array of box corner values -!> \param X Fractional distance from the lower to upper x coordinate -!> \param Y Fractional distance from the lower to upper y coordinate -!> \param Z Fractional distance from the lower to upper z coordinate - -SUBROUTINE TRILIN_INTERP(V,A,X,Y,Z) - -REAL(EB), INTENT(IN), DIMENSION(0:1,0:1,0:1) :: A -REAL(EB), INTENT(IN) :: X,Y,Z -REAL(EB), INTENT(OUT) :: V -REAL(EB), DIMENSION(0:1,0:1,0:1) :: WGT -REAL(EB) :: XX,YY,ZZ - -XX = 1._EB-X -YY = 1._EB-Y -ZZ = 1._EB-Z - -WGT(0,0,0) = XX * YY * ZZ -WGT(1,0,0) = X * YY * ZZ -WGT(0,1,0) = XX * Y * ZZ -WGT(0,0,1) = XX * YY * Z -WGT(1,0,1) = X * YY * Z -WGT(0,1,1) = XX * Y * Z -WGT(1,1,0) = X * Y * ZZ -WGT(1,1,1) = X * Y * Z - -V = SUM(A*WGT) - -END SUBROUTINE TRILIN_INTERP - - !> \brief Calculate the value of polynomial function. !> \param N Number of coefficients in the polynomial !> \param TEMP The independent variable diff --git a/Source/main.f90 b/Source/main.f90 index b0c429037a9..c0458c16e8f 100644 --- a/Source/main.f90 +++ b/Source/main.f90 @@ -52,10 +52,9 @@ PROGRAM FDS INTEGER :: LO10,NM,IZERO,ANG_INC_COUNTER REAL(EB) :: T,DT,TNOW REAL :: CPUTIME -REAL(EB), ALLOCATABLE, DIMENSION(:) :: TC_GLB,TC_LOC,DT_NEW,TI_LOC,TI_GLB, & - DSUM_ALL,PSUM_ALL,USUM_ALL,DSUM_ALL_LOCAL,PSUM_ALL_LOCAL,USUM_ALL_LOCAL +REAL(EB), ALLOCATABLE, DIMENSION(:) :: TC_GLB,TC_LOC,DT_NEW,TI_LOC,TI_GLB,DSUM_ALL,PSUM_ALL,USUM_ALL REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: TC2_GLB,TC2_LOC -LOGICAL, ALLOCATABLE, DIMENSION(:,:) :: CONNECTED_ZONES_GLOBAL,CONNECTED_ZONES_LOCAL +LOGICAL, ALLOCATABLE, DIMENSION(:,:) :: CONNECTED_ZONES_ALL LOGICAL, ALLOCATABLE, DIMENSION(:) :: STATE_GLB,STATE_LOC INTEGER :: ITER TYPE (MESH_TYPE), POINTER :: M,M4 @@ -519,7 +518,7 @@ PROGRAM FDS ! If there are zones and HVAC pass PSUM -IF (HVAC_SOLVE .AND. N_ZONE>0) CALL EXCHANGE_DIVERGENCE_INFO +IF (HVAC_SOLVE .AND. N_ZONE>0 .AND. .NOT.SOLID_PHASE_ONLY) CALL EXCHANGE_DIVERGENCE_INFO ! Make an initial dump of global output quantities @@ -687,7 +686,7 @@ PROGRAM FDS ! If there are pressure ZONEs, exchange integrated quantities mesh to mesh for use in the divergence calculation - IF (N_ZONE>0) CALL EXCHANGE_DIVERGENCE_INFO + IF (N_ZONE>0 .AND. .NOT.SOLID_PHASE_ONLY) CALL EXCHANGE_DIVERGENCE_INFO ! Update global pressure matrices after zone connections @@ -886,7 +885,7 @@ PROGRAM FDS ! Exchange global pressure zone information - IF (N_ZONE>0) CALL EXCHANGE_DIVERGENCE_INFO + IF (N_ZONE>0 .AND. .NOT.SOLID_PHASE_ONLY) CALL EXCHANGE_DIVERGENCE_INFO ! Update global pressure matrices after zone connections @@ -1390,11 +1389,7 @@ SUBROUTINE MPI_INITIALIZATION_CHORES(TASK_NUMBER) ALLOCATE(DSUM_ALL(N_ZONE),STAT=IZERO) ALLOCATE(PSUM_ALL(N_ZONE),STAT=IZERO) ALLOCATE(USUM_ALL(N_ZONE),STAT=IZERO) - ALLOCATE(CONNECTED_ZONES_GLOBAL(0:N_ZONE,0:N_ZONE),STAT=IZERO) - ALLOCATE(DSUM_ALL_LOCAL(N_ZONE),STAT=IZERO) - ALLOCATE(PSUM_ALL_LOCAL(N_ZONE),STAT=IZERO) - ALLOCATE(USUM_ALL_LOCAL(N_ZONE),STAT=IZERO) - ALLOCATE(CONNECTED_ZONES_LOCAL(0:N_ZONE,0:N_ZONE),STAT=IZERO) + ALLOCATE(CONNECTED_ZONES_ALL(0:N_ZONE,0:N_ZONE),STAT=IZERO) ENDIF ALLOCATE(CONNECTED_ZONES(0:N_ZONE,0:N_ZONE,NMESHES),STAT=IZERO) @@ -1818,32 +1813,27 @@ SUBROUTINE EXCHANGE_DIVERGENCE_INFO TNOW = CURRENT_TIME() -CONNECTED_ZONES_LOCAL = .FALSE. +CONNECTED_ZONES_ALL = .FALSE. DO IPZ=1,N_ZONE - DSUM_ALL_LOCAL(IPZ) = 0._EB - PSUM_ALL_LOCAL(IPZ) = 0._EB - USUM_ALL_LOCAL(IPZ) = 0._EB + DSUM_ALL(IPZ) = 0._EB + PSUM_ALL(IPZ) = 0._EB + USUM_ALL(IPZ) = 0._EB DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX - DSUM_ALL_LOCAL(IPZ) = DSUM_ALL_LOCAL(IPZ) + DSUM(IPZ,NM) - PSUM_ALL_LOCAL(IPZ) = PSUM_ALL_LOCAL(IPZ) + PSUM(IPZ,NM) - USUM_ALL_LOCAL(IPZ) = USUM_ALL_LOCAL(IPZ) + USUM(IPZ,NM) + DSUM_ALL(IPZ) = DSUM_ALL(IPZ) + DSUM(IPZ,NM) + PSUM_ALL(IPZ) = PSUM_ALL(IPZ) + PSUM(IPZ,NM) + USUM_ALL(IPZ) = USUM_ALL(IPZ) + USUM(IPZ,NM) DO IOPZ=0,N_ZONE - IF (CONNECTED_ZONES(IPZ,IOPZ,NM)) CONNECTED_ZONES_LOCAL(IPZ,IOPZ) = .TRUE. + IF (CONNECTED_ZONES(IPZ,IOPZ,NM)) CONNECTED_ZONES_ALL(IPZ,IOPZ) = .TRUE. ENDDO ENDDO ENDDO IF (N_MPI_PROCESSES>1) THEN - CALL MPI_ALLREDUCE(DSUM_ALL_LOCAL(1),DSUM_ALL(1),N_ZONE,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IERR) - CALL MPI_ALLREDUCE(PSUM_ALL_LOCAL(1),PSUM_ALL(1),N_ZONE,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IERR) - CALL MPI_ALLREDUCE(USUM_ALL_LOCAL(1),USUM_ALL(1),N_ZONE,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IERR) - CALL MPI_ALLREDUCE(CONNECTED_ZONES_LOCAL(0,0),CONNECTED_ZONES_GLOBAL(0,0),(N_ZONE+1)**2,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,IERR) -ELSE - DSUM_ALL = DSUM_ALL_LOCAL - PSUM_ALL = PSUM_ALL_LOCAL - USUM_ALL = USUM_ALL_LOCAL - CONNECTED_ZONES_GLOBAL = CONNECTED_ZONES_LOCAL + CALL MPI_ALLREDUCE(MPI_IN_PLACE,DSUM_ALL(1),N_ZONE,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IERR) + CALL MPI_ALLREDUCE(MPI_IN_PLACE,PSUM_ALL(1),N_ZONE,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IERR) + CALL MPI_ALLREDUCE(MPI_IN_PLACE,USUM_ALL(1),N_ZONE,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IERR) + CALL MPI_ALLREDUCE(MPI_IN_PLACE,CONNECTED_ZONES_ALL(0,0),(N_ZONE+1)**2,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,IERR) ENDIF DO IPZ=1,N_ZONE @@ -1851,7 +1841,7 @@ SUBROUTINE EXCHANGE_DIVERGENCE_INFO DSUM(IPZ,NM) = DSUM_ALL(IPZ) PSUM(IPZ,NM) = PSUM_ALL(IPZ) USUM(IPZ,NM) = USUM_ALL(IPZ) - CONNECTED_ZONES(IPZ,:,NM) = CONNECTED_ZONES_GLOBAL(IPZ,:) + CONNECTED_ZONES(IPZ,:,NM) = CONNECTED_ZONES_ALL(IPZ,:) ENDDO ENDDO diff --git a/Source/part.f90 b/Source/part.f90 index 7f2779f4ec0..9bf57a0621f 100644 --- a/Source/part.f90 +++ b/Source/part.f90 @@ -1215,7 +1215,7 @@ SUBROUTINE INSERT_VOLUMETRIC_PARTICLES LP%DX = DX(II) LP%DY = DY(JJ) LP%DZ = DZ(KK) - LP%INITIALIZATION_INDEX = INIT_INDEX + LP%INIT_INDEX = INIT_INDEX ! Initialize particle properties @@ -1399,7 +1399,6 @@ SUBROUTINE INSERT_VOLUMETRIC_PARTICLES LP%PWT = LP%PWT*PWT0 ENDIF IF (ANY(IN%PATH_RAMP_INDEX>0)) LP%PATH_PARTICLE=.TRUE. - LP%INIT_INDEX = INIT_INDEX ENDDO ENDIF @@ -2512,7 +2511,7 @@ END SUBROUTINE MOVE_ON_SOLID SUBROUTINE MOVE_IN_GAS USE PHYSICAL_FUNCTIONS, ONLY : DRAG, GET_VISCOSITY, SURFACE_DENSITY -USE MATH_FUNCTIONS, ONLY : EVALUATE_RAMP, RANDOM_CHOICE, BOX_MULLER, MESH_TO_PARTICLE, PARTICLE_TO_MESH +USE MATH_FUNCTIONS, ONLY : EVALUATE_RAMP, RANDOM_CHOICE, BOX_MULLER USE SOOT_ROUTINES, ONLY: DROPLET_SCRUBBING REAL(EB) :: UBAR,VBAR,WBAR,UREL,VREL,WREL,QREL,RHO_G,TMP_G,MU_FILM, & U_OLD,V_OLD,W_OLD,ZZ_GET(1:N_TRACKED_SPECIES),WAKE_VEL,DROP_VOL_FRAC,RE_WAKE,& @@ -2521,10 +2520,9 @@ SUBROUTINE MOVE_IN_GAS GX_LOC,GY_LOC,GZ_LOC,DRAG_MAX(3)=0._EB,K_SGS,U_P,KN,M_DOT,& EMBER_DENSITY,EMBER_VOLUME=0._EB,ACCEL_X,ACCEL_Y,ACCEL_Z,& LP_FORCE,FACE_VOLS(2,2,2),VEL_G_INT(3),VOL_WGT(2,2,2),& - EMBER_PACKING_RATIO,LOCAL_PACKING_RATIO,LPC_GEOM_FACTOR,& - X_WGT,Y_WGT,Z_WGT,X_WGT2,Y_WGT2,Z_WGT2 -REAL(EB) :: WGT(2,2,2,3) -REAL(EB), POINTER, DIMENSION(:,:,:) :: FV_D=>NULL(),VEL_G=>NULL() + EMBER_PACKING_RATIO,LOCAL_PACKING_RATIO,LPC_GEOM_FACTOR +REAL(EB) :: WGT(2,2,2,3),VEL_G(2,2,2) +REAL(EB), POINTER, DIMENSION(:,:,:) :: FV_D=>NULL() REAL(EB), SAVE :: BETA INTEGER :: IIX,JJY,KKZ,IL,JL,KL,AXIS,N_LPC2 LOGICAL :: STUCK=.FALSE. @@ -2556,46 +2554,29 @@ SUBROUTINE MOVE_IN_GAS ENDIF ENDIF -IF (ICC>0) THEN - WGT=0._EB - DO AXIS=IAXIS,KAXIS - IL = IIX; JL = JJY; KL = KKZ - IF (AXIS==IAXIS) THEN - VEL_G => U - IL = FLOOR(XI) - ELSEIF (AXIS==JAXIS) THEN - VEL_G => V - JL = FLOOR(YJ) - ELSEIF (AXIS==KAXIS) THEN - VEL_G => W - KL = FLOOR(ZK) - ENDIF +WGT=0._EB +DO AXIS=IAXIS,KAXIS + IL = IIX; JL = JJY; KL = KKZ + IF (AXIS==IAXIS) THEN + IL = FLOOR(XI) + VEL_G = U(IL:IL+1,JL:JL+1,KL:KL+1) + ELSEIF (AXIS==JAXIS) THEN + JL = FLOOR(YJ) + VEL_G = V(IL:IL+1,JL:JL+1,KL:KL+1) + ELSEIF (AXIS==KAXIS) THEN + KL = FLOOR(ZK) + VEL_G = W(IL:IL+1,JL:JL+1,KL:KL+1) + ENDIF + IF (ICC>0) THEN CALL GET_FACE_IDW(AXIS,IL,JL,KL,BC%X,BC%Y,BC%Z,WGT(:,:,:,AXIS)) - VEL_G_INT(AXIS) = SUM(VEL_G(IL:IL+1,JL:JL+1,KL:KL+1)*WGT(:,:,:,AXIS)) - ENDDO - UBAR = VEL_G_INT(IAXIS) - VBAR = VEL_G_INT(JAXIS) - WBAR = VEL_G_INT(KAXIS) -ELSE - X_WGT = XI+.5_EB-IIX - Y_WGT = YJ+.5_EB-JJY - Z_WGT = ZK+.5_EB-KKZ - X_WGT2 = XI-FLOOR(XI) - Y_WGT2 = YJ-FLOOR(YJ) - Z_WGT2 = ZK-FLOOR(ZK) - - IF (X_WGT>=0.5_EB .AND. CELL(IC_OLD)%WALL_INDEX(-1)>0) X_WGT = 1._EB - IF (X_WGT< 0.5_EB .AND. CELL(IC_OLD)%WALL_INDEX( 1)>0) X_WGT = 0._EB - IF (Y_WGT>=0.5_EB .AND. CELL(IC_OLD)%WALL_INDEX(-2)>0) Y_WGT = 1._EB - IF (Y_WGT< 0.5_EB .AND. CELL(IC_OLD)%WALL_INDEX( 2)>0) Y_WGT = 0._EB - IF (Z_WGT>=0.5_EB .AND. CELL(IC_OLD)%WALL_INDEX(-3)>0) Z_WGT = 1._EB - IF (Z_WGT< 0.5_EB .AND. CELL(IC_OLD)%WALL_INDEX( 3)>0) Z_WGT = 0._EB - - CALL MESH_TO_PARTICLE(UBAR,U,IIG_OLD-1,JJY,KKZ,X_WGT2,Y_WGT,Z_WGT) - CALL MESH_TO_PARTICLE(VBAR,V,IIX,JJG_OLD-1,KKZ,X_WGT,Y_WGT2,Z_WGT) - CALL MESH_TO_PARTICLE(WBAR,W,IIX,JJY,KKG_OLD-1,X_WGT,Y_WGT,Z_WGT2) -ENDIF - + ELSE + CALL GET_FACE_TLW(AXIS,IL,JL,KL,BC%X,BC%Y,BC%Z,WGT(:,:,:,AXIS),VEL_G) + ENDIF + VEL_G_INT(AXIS) = SUM(VEL_G*WGT(:,:,:,AXIS)) +ENDDO +UBAR = VEL_G_INT(IAXIS) +VBAR = VEL_G_INT(JAXIS) +WBAR = VEL_G_INT(KAXIS) ! If the particle has a path, just follow the path and return @@ -2962,32 +2943,26 @@ SUBROUTINE MOVE_IN_GAS LP%ACCEL_Y = LP%ACCEL_Y + ACCEL_Y LP%ACCEL_Z = LP%ACCEL_Z + ACCEL_Z -IF (ICC>0) THEN - DO AXIS=IAXIS,KAXIS - IL = IIX; JL = JJY; KL = KKZ - IF (AXIS == IAXIS) THEN - LP_FORCE = ACCEL_X/LP%RVC - IL = FLOOR(XI) - FV_D => FVX_D - ELSEIF (AXIS == JAXIS) THEN - LP_FORCE = ACCEL_Y/LP%RVC - JL = FLOOR(YJ) - FV_D => FVY_D - ELSEIF (AXIS == KAXIS) THEN - LP_FORCE = ACCEL_Z/LP%RVC - KL = FLOOR(ZK) - FV_D => FVZ_D - ENDIF - CALL GET_FACE_VOLUMES(AXIS,IL,JL,KL,FACE_VOLS) - VOL_WGT = FACE_VOLS*WGT(:,:,:,AXIS) - IF (ANY(VOL_WGT>TWO_EPSILON_EB)) VOL_WGT=WGT(:,:,:,AXIS)/SUM(VOL_WGT) - FV_D(IL:IL+1, JL:JL+1, KL:KL+1) = FV_D(IL:IL+1, JL:JL+1, KL:KL+1) - LP_FORCE*VOL_WGT - ENDDO -ELSE - CALL PARTICLE_TO_MESH(ACCEL_X,FVX_D,IIG_OLD-1,JJY,KKZ,X_WGT2,Y_WGT,Z_WGT) - CALL PARTICLE_TO_MESH(ACCEL_Y,FVY_D,IIX,JJG_OLD-1,KKZ,X_WGT,Y_WGT2,Z_WGT) - CALL PARTICLE_TO_MESH(ACCEL_Z,FVZ_D,IIX,JJY,KKG_OLD-1,X_WGT,Y_WGT,Z_WGT2) -ENDIF +DO AXIS=IAXIS,KAXIS + IL = IIX; JL = JJY; KL = KKZ + IF (AXIS == IAXIS) THEN + LP_FORCE = ACCEL_X/LP%RVC + IL = FLOOR(XI) + FV_D => FVX_D + ELSEIF (AXIS == JAXIS) THEN + LP_FORCE = ACCEL_Y/LP%RVC + JL = FLOOR(YJ) + FV_D => FVY_D + ELSEIF (AXIS == KAXIS) THEN + LP_FORCE = ACCEL_Z/LP%RVC + KL = FLOOR(ZK) + FV_D => FVZ_D + ENDIF + CALL GET_FACE_VOLUMES(AXIS,IL,JL,KL,FACE_VOLS) + VOL_WGT = FACE_VOLS*WGT(:,:,:,AXIS) + IF (ANY(VOL_WGT>TWO_EPSILON_EB)) VOL_WGT=WGT(:,:,:,AXIS)/SUM(VOL_WGT) + FV_D(IL:IL+1, JL:JL+1, KL:KL+1) = FV_D(IL:IL+1, JL:JL+1, KL:KL+1) - LP_FORCE*VOL_WGT +ENDDO ! store C_DRAG for output @@ -3219,25 +3194,21 @@ SUBROUTINE GET_FACE_IDW(AXIS,I,J,K,P_X,P_Y,P_Z,IDW) ELSE XYZ_INT = (/X_F(II),Y_F(JJ),Z_F(KK)/) ENDIF - IF (DIST>TWO_EPSILON_EB) THEN - IDW(II-I+1,JJ-J+1,KK-K+1) = 0._EB + DIST = NORM2((/P_X,P_Y,P_Z/)-XYZ_INT) + ! Special case where location is directly on face + IF (DIST0) D_WGT = 0._EB - IF (CC_IBM .AND. D_WGT>0._EB) THEN - IF(FCVAR(II,JJ,KK,CC_FGSC,AXIS)==CC_SOLID) D_WGT = 0._EB - ENDIF - IDW(II-I+1,JJ-J+1,KK-K+1) = D_WGT + D_WGT = 1._EB/DIST**6._EB + ENDIF + ! face is solid + IF(CELL(CELL_INDEX(II,JJ,KK))%WALL_INDEX(AXIS)>0) D_WGT = 0._EB + IF (CC_IBM .AND. D_WGT>0._EB) THEN + IF(FCVAR(II,JJ,KK,CC_FGSC,AXIS)==CC_SOLID) D_WGT = 0._EB ENDIF + IDW(II-I+1,JJ-J+1,KK-K+1) = D_WGT ENDDO ENDDO ENDDO FACE_LOOP @@ -3264,6 +3235,126 @@ SUBROUTINE GET_FACE_IDW(AXIS,I,J,K,P_X,P_Y,P_Z,IDW) END SUBROUTINE GET_FACE_IDW +!> \brief Get Tri-Linear interpolation Weight (TLW) values for nearest gas faces +!> \param AXIS The axis of the face quantity +!> \param I The lower x index +!> \param J The lower y index +!> \param K The lower z index +!> \param P_X Sample point location in x +!> \param P_Y Sample point location in y +!> \param P_Z Sample point location in z +!> \param TLW 2x2x2 matrix of weights for cartesian faces +!> \param V 2x2x2 array of velocity box corner values + +SUBROUTINE GET_FACE_TLW(AXIS,I,J,K,P_X,P_Y,P_Z,TLW,V) + +REAL(EB), INTENT(IN) :: P_X,P_Y,P_Z +REAL(EB), INTENT(OUT) :: TLW(0:1,0:1,0:1) +INTEGER, INTENT(IN) :: AXIS,I,J,K +REAL(EB), INTENT(INOUT) :: V(2,2,2) +REAL(EB) :: P,PP,R,RR,S,SS +INTEGER :: IWC(-3:3) +REAL(EB), POINTER :: X_F(:),Y_F(:),Z_F(:) + +TLW=0._EB + +X_F=>XC +Y_F=>YC +Z_F=>ZC +SELECT CASE(AXIS) + CASE(IAXIS); X_F=>X + CASE(JAXIS); Y_F=>Y + CASE(KAXIS); Z_F=>Z +END SELECT + +P = (P_X-X_F(I))/(X_F(I+1)-X_F(I)) +R = (P_Y-Y_F(J))/(Y_F(J+1)-Y_F(J)) +S = (P_Z-Z_F(K))/(Z_F(K+1)-Z_F(K)) + +IWC = CELL(IC_OLD)%WALL_INDEX +IF (NEAR_WALL_PARTICLE_INTERPOLATION) THEN + + IF (AXIS/=IAXIS .AND. IIG_OLD> I .AND. IWC(-1)>0) THEN + IF (WALL(IWC(-1))%BOUNDARY_TYPE==SOLID_BOUNDARY) THEN + P=(P_X-X_F(I))/(BOUNDARY_COORD(WALL(IWC(-1))%BC_INDEX)%X-X_F(I)) + SELECT CASE(AXIS) + CASE(JAXIS); V(1,:,:)=SURFACE(WALL(IWC(-1))%SURF_INDEX)%VEL_T(1) + CASE(KAXIS); V(1,:,:)=SURFACE(WALL(IWC(-1))%SURF_INDEX)%VEL_T(2) + END SELECT + ENDIF + ENDIF + IF (AXIS/=IAXIS .AND. IIG_OLD==I .AND. IWC( 1)>0) THEN + IF (WALL(IWC( 1))%BOUNDARY_TYPE==SOLID_BOUNDARY) THEN + P=(P_X-X_F(I))/(BOUNDARY_COORD(WALL(IWC( 1))%BC_INDEX)%X-X_F(I)) + SELECT CASE(AXIS) + CASE(JAXIS); V(2,:,:)=SURFACE(WALL(IWC( 1))%SURF_INDEX)%VEL_T(1) + CASE(KAXIS); V(2,:,:)=SURFACE(WALL(IWC( 1))%SURF_INDEX)%VEL_T(2) + END SELECT + ENDIF + ENDIF + + IF (AXIS/=JAXIS .AND. JJG_OLD> J .AND. IWC(-2)>0) THEN + IF (WALL(IWC(-2))%BOUNDARY_TYPE==SOLID_BOUNDARY) THEN + R=(P_Y-Y_F(J))/(BOUNDARY_COORD(WALL(IWC(-2))%BC_INDEX)%Y-Y_F(J)) + SELECT CASE(AXIS) + CASE(IAXIS); V(:,1,:)=SURFACE(WALL(IWC(-2))%SURF_INDEX)%VEL_T(1) + CASE(KAXIS); V(:,1,:)=SURFACE(WALL(IWC(-2))%SURF_INDEX)%VEL_T(2) + END SELECT + ENDIF + ENDIF + IF (AXIS/=JAXIS .AND. JJG_OLD==J .AND. IWC( 2)>0) THEN + IF (WALL(IWC( 2))%BOUNDARY_TYPE==SOLID_BOUNDARY) THEN + R=(P_Y-Y_F(J))/(BOUNDARY_COORD(WALL(IWC( 2))%BC_INDEX)%Y-Y_F(J)) + SELECT CASE(AXIS) + CASE(IAXIS); V(:,2,:)=SURFACE(WALL(IWC( 2))%SURF_INDEX)%VEL_T(1) + CASE(KAXIS); V(:,2,:)=SURFACE(WALL(IWC( 2))%SURF_INDEX)%VEL_T(2) + END SELECT + ENDIF + ENDIF + + IF (AXIS/=KAXIS .AND. KKG_OLD> K .AND. IWC(-3)>0) THEN + IF (WALL(IWC(-3))%BOUNDARY_TYPE==SOLID_BOUNDARY) THEN + S=(P_Z-Z_F(K))/(BOUNDARY_COORD(WALL(IWC(-3))%BC_INDEX)%Z-Z_F(K)) + SELECT CASE(AXIS) + CASE(IAXIS); V(:,:,1)=SURFACE(WALL(IWC(-3))%SURF_INDEX)%VEL_T(1) + CASE(JAXIS); V(:,:,1)=SURFACE(WALL(IWC(-3))%SURF_INDEX)%VEL_T(2) + END SELECT + ENDIF + ENDIF + IF (AXIS/=KAXIS .AND. KKG_OLD==K .AND. IWC( 3)>0) THEN + IF (WALL(IWC( 3))%BOUNDARY_TYPE==SOLID_BOUNDARY) THEN + S=(P_Z-Z_F(K))/(BOUNDARY_COORD(WALL(IWC( 3))%BC_INDEX)%Z-Z_F(K)) + SELECT CASE(AXIS) + CASE(IAXIS); V(:,:,2)=SURFACE(WALL(IWC( 3))%SURF_INDEX)%VEL_T(1) + CASE(JAXIS); V(:,:,2)=SURFACE(WALL(IWC( 3))%SURF_INDEX)%VEL_T(2) + END SELECT + ENDIF + ENDIF + +ELSE + IF (AXIS/=IAXIS .AND. IIG_OLD> I .AND. IWC(-1)>0) P = 1._EB + IF (AXIS/=IAXIS .AND. IIG_OLD==I .AND. IWC( 1)>0) P = 0._EB + IF (AXIS/=JAXIS .AND. JJG_OLD> J .AND. IWC(-2)>0) R = 1._EB + IF (AXIS/=JAXIS .AND. JJG_OLD==J .AND. IWC( 2)>0) R = 0._EB + IF (AXIS/=KAXIS .AND. KKG_OLD> K .AND. IWC(-3)>0) S = 1._EB + IF (AXIS/=KAXIS .AND. KKG_OLD==K .AND. IWC( 3)>0) S = 0._EB +ENDIF + +PP = 1._EB-P +RR = 1._EB-R +SS = 1._EB-S + +TLW(0,0,0) = PP * RR * SS +TLW(1,0,0) = P * RR * SS +TLW(0,1,0) = PP * R * SS +TLW(0,0,1) = PP * RR * S +TLW(1,0,1) = P * RR * S +TLW(0,1,1) = PP * R * S +TLW(1,1,0) = P * R * SS +TLW(1,1,1) = P * R * S + +END SUBROUTINE GET_FACE_TLW + !> \brief Return face volumes for distribution of quantity onto faces !> \param AXIS The axis for the face centered quantity !> \param I The lower x index @@ -3282,9 +3373,9 @@ SUBROUTINE GET_FACE_VOLUMES(AXIS,I,J,K,FACE_VOLS) DX2 => DY DX3 => DZ SELECT CASE (AXIS) -CASE(IAXIS); DX1 => DXN -CASE(JAXIS); DX2 => DYN -CASE(KAXIS); DX3 => DZN + CASE(IAXIS); DX1 => DXN + CASE(JAXIS); DX2 => DYN + CASE(KAXIS); DX3 => DZN END SELECT DO KK=K,K+1 @@ -3466,7 +3557,7 @@ SUBROUTINE PARTICLE_MASS_ENERGY_TRANSFER(T,DT,NM) RETURN ELSE ALLOCATE(PART_WARNING(NLP)) - PART_WARNING(NLP)=0 + PART_WARNING=0 ENDIF ! Working arrays @@ -4064,7 +4155,16 @@ SUBROUTINE PARTICLE_MASS_ENERGY_TRANSFER(T,DT,NM) CALL GET_TEMPERATURE(TMP_G_NEW,H_NEW/M_GAS_NEW,ZZ_GET2) IF (TMP_G_NEW < 0._EB) THEN DT_SUBSTEP = DT_SUBSTEP * 0.5_EB - CYCLE TIME_ITERATION_LOOP + IF (DT_SUBSTEP <= 0.00001_EB*DT) THEN + DT_SUBSTEP = DT_SUBSTEP * 2.0_EB + TMP_G_NEW = 1._EB + IF (.NOT. BTEST(PART_WARNING(IP),3)) THEN + WRITE(LU_ERR,'(A,I0,A,I0,A,I0)') 'WARNING TMP_G_N < 0. Mesh: ',NM,'Particle: ',IP,' Tag: ',LP%TAG + PART_WARNING(IP) = IBSET(PART_WARNING(IP),3) + ENDIF + ELSE + CYCLE TIME_ITERATION_LOOP + ENDIF ENDIF ! Limit gas temperature change @@ -4073,9 +4173,9 @@ SUBROUTINE PARTICLE_MASS_ENERGY_TRANSFER(T,DT,NM) DT_SUBSTEP = DT_SUBSTEP * 0.5_EB IF (DT_SUBSTEP <= 0.00001_EB*DT) THEN DT_SUBSTEP = DT_SUBSTEP * 2.0_EB - IF (.NOT. BTEST(PART_WARNING(IP),3)) THEN - WRITE(LU_ERR,'(A,I0,A,I0,A,I0)') 'WARNING Delta TMP_G.Mesh: ',NM,'Particle: ',IP,' Tag: ',LP%TAG - PART_WARNING(IP) = IBSET(PART_WARNING(IP),3) + IF (.NOT. BTEST(PART_WARNING(IP),4)) THEN + WRITE(LU_ERR,'(A,I0,A,I0,A,I0)') 'WARNING Delta TMP_G. Mesh: ',NM,'Particle: ',IP,' Tag: ',LP%TAG + PART_WARNING(IP) = IBSET(PART_WARNING(IP),4) ENDIF ELSE CYCLE TIME_ITERATION_LOOP @@ -4089,9 +4189,9 @@ SUBROUTINE PARTICLE_MASS_ENERGY_TRANSFER(T,DT,NM) DT_SUBSTEP = DT_SUBSTEP * 0.5_EB IF (DT_SUBSTEP <= 0.00001_EB*DT) THEN DT_SUBSTEP = DT_SUBSTEP * 2.0_EB - IF (.NOT. BTEST(PART_WARNING(IP),4)) THEN - WRITE(LU_ERR,'(A,I0,A,I0,A,I0)') 'WARNING TMP_G_N < TMP_D_N.Mesh: ',NM,'Particle: ',IP,' Tag: ',LP%TAG - PART_WARNING(IP) = IBSET(PART_WARNING(IP),4) + IF (.NOT. BTEST(PART_WARNING(IP),5)) THEN + WRITE(LU_ERR,'(A,I0,A,I0,A,I0)') 'WARNING TMP_G_N < TMP_D_N. Mesh: ',NM,'Particle: ',IP,' Tag: ',LP%TAG + PART_WARNING(IP) = IBSET(PART_WARNING(IP),5) ENDIF ELSE CYCLE TIME_ITERATION_LOOP diff --git a/Source/radi.f90 b/Source/radi.f90 index de0dbef45be..686fef277ed 100644 --- a/Source/radi.f90 +++ b/Source/radi.f90 @@ -2788,8 +2788,8 @@ SUBROUTINE INIT_RADIATION USE RADCAL_CALC USE WSGG_ARRAYS REAL(EB) :: THETAUP,THETALOW,PHIUP,PHILOW,F_THETA,PLANCK_C2,KSI,LT,RCRHO,YY,YY2,BBF,AP0,AMEAN,RADIANCE,TRANSMISSIVITY,X_N2,& - THETA,PHI -INTEGER :: N,I,J,K,IPC,IZERO,NN,NI,II,JJ,IIM,JJM,IBND,NS,NS2,NRA,NSB,RADCAL_TEMP(16)=0,RCT_SKIP=-1,OR_IN,I1,I2,IO + THETA,PHI,DLO +INTEGER :: N,I,J,K,IPC,IZERO,NN,NI,II,JJ,IIM,JJM,IBND,NS,NS2,NRA,NSB,RADCAL_TEMP(16)=0,RCT_SKIP=-1,IO TYPE (LAGRANGIAN_PARTICLE_CLASS_TYPE), POINTER :: LPC REAL(EB), ALLOCATABLE, DIMENSION(:) :: COSINE_ARRAY TYPE (RAD_FILE_TYPE), POINTER :: RF @@ -3389,44 +3389,28 @@ SUBROUTINE INIT_RADIATION ENDDO ! Determine angle factors for Lagrangian particles with ORIENTATION +! COSINE_ARRAY holds the cosines of the angles formed by the orientation vector and the radiation directions. +! DLO is the integral of the orientation vector dotted with the directional solid angle of the radiation directions. +! VIEW_ANGLE_FACTOR is the reduction of the radiation due to a view angle less than 180, like a narrow field of view radiometer. IF (SOLID_PARTICLES) THEN - PARTICLE_CLASS_LOOP: DO IPC=1,N_LAGRANGIAN_CLASSES - LPC => LAGRANGIAN_PARTICLE_CLASS(IPC) - IF (LPC%N_ORIENTATION==0) CYCLE PARTICLE_CLASS_LOOP - I1 = LPC%ORIENTATION_INDEX - I2 = LPC%ORIENTATION_INDEX+LPC%N_ORIENTATION-1 - ALLOCATE(COSINE_ARRAY(I1:I2)) - ANGLE_LOOP: DO N=1,NRA - ORIENTATION_LOOP: DO OR_IN=I1,I2 - COSINE_ARRAY(OR_IN) = ORIENTATION_VECTOR(1,OR_IN)*DLX(N) + & - ORIENTATION_VECTOR(2,OR_IN)*DLY(N) + & - ORIENTATION_VECTOR(3,OR_IN)*DLZ(N) - ENDDO ORIENTATION_LOOP - ENDDO ANGLE_LOOP - DEALLOCATE(COSINE_ARRAY) - ENDDO PARTICLE_CLASS_LOOP - ALLOCATE(COSINE_ARRAY(1:NRA)) ALLOCATE(NEAREST_RADIATION_ANGLE(N_ORIENTATION_VECTOR)) - ALLOCATE(VIEW_ANGLE_AREA(N_ORIENTATION_VECTOR)) - VIEW_ANGLE_AREA = 0._EB + ALLOCATE(VIEW_ANGLE_FACTOR(N_ORIENTATION_VECTOR)) + VIEW_ANGLE_FACTOR = 0._EB DO IO=1,N_ORIENTATION_VECTOR + DLO = 0._EB DO N=1,NRA - COSINE_ARRAY(N) = ORIENTATION_VECTOR(1,IO)*DLX(N) + & - ORIENTATION_VECTOR(2,IO)*DLY(N) + & - ORIENTATION_VECTOR(3,IO)*DLZ(N) - IF (-(ORIENTATION_VECTOR(1,IO)*DLANG(1,N) + & - ORIENTATION_VECTOR(2,IO)*DLANG(2,N) + & - ORIENTATION_VECTOR(3,IO)*DLANG(3,N)) > ORIENTATION_VIEW_ANGLE(IO)) & - VIEW_ANGLE_AREA(IO) = VIEW_ANGLE_AREA(IO) - COSINE_ARRAY(N) + COSINE_ARRAY(N) = ORIENTATION_VECTOR(1,IO)*DLANG(1,N) + & + ORIENTATION_VECTOR(2,IO)*DLANG(2,N) + & + ORIENTATION_VECTOR(3,IO)*DLANG(3,N) + IF (-COSINE_ARRAY(N) > COS_HALF_VIEW_ANGLE(IO)) & + DLO = DLO - (ORIENTATION_VECTOR(1,IO)*DLX(N) + ORIENTATION_VECTOR(2,IO)*DLY(N) + ORIENTATION_VECTOR(3,IO)*DLZ(N)) ENDDO NEAREST_RADIATION_ANGLE(IO) = MINLOC(COSINE_ARRAY,DIM=1) - VIEW_ANGLE_AREA(IO) = PI/VIEW_ANGLE_AREA(IO) + VIEW_ANGLE_FACTOR(IO) = PI/DLO ENDDO - DEALLOCATE(COSINE_ARRAY) - ENDIF ! Allocate array needed by angle-specific RADF output files @@ -3483,7 +3467,7 @@ SUBROUTINE RADIATION_FVM USE PHYSICAL_FUNCTIONS, ONLY : GET_VOLUME_FRACTION, GET_MASS_FRACTION REAL(EB) :: RAP, AX, AXU, AXD, AY, AYU, AYD, AZ, AZU, AZD, VC, RU, RD, RP, AFD, & ILXU, ILYU, ILZU, QVAL, BBF, BBFA, NCSDROP, RSA_RAT,EFLUX,SOOT_MASS_FRACTION, & - AIU_SUM,A_SUM,VOL,VC1,AY1,AZ1,COS_DL,AILFU, & + AIU_SUM,A_SUM,VOL,VC1,AY1,AZ1,DLO,COS_DLO,AILFU, & RAD_Q_SUM_PARTIAL,KFST4_SUM_PARTIAL,ALPHA_CC INTEGER :: N,NN,IIG,JJG,KKG,I,J,K,IW,ICF,II,JJ,KK,IOR,IC,IWUP,IWDOWN, & @@ -4402,15 +4386,17 @@ SUBROUTINE RADIATION_FVM ENDDO OTHER_WALL_LOOP ENDDO INTERPOLATION_LOOP - ! Compute projected intensity on particles + ! Compute projected intensity on particles with a specified ORIENTATION - IF (SOLID_PARTICLES) THEN + IF (ORIENTED_PARTICLES) THEN PARTICLE_RADIATION_LOOP: DO IP=1,NLP LP => LAGRANGIAN_PARTICLE(IP) LPC => LAGRANGIAN_PARTICLE_CLASS(LP%CLASS_INDEX) + IF (LPC%N_ORIENTATION==0) CYCLE PARTICLE_RADIATION_LOOP BC => BOUNDARY_COORD(LP%BC_INDEX) - IF (LP%INITIALIZATION_INDEX > 0) THEN - IN => INITIALIZATION(LP%INITIALIZATION_INDEX) + TEMP_ORIENTATION(1:3) = ORIENTATION_VECTOR(1:3,LP%ORIENTATION_INDEX) + IF (LP%INIT_INDEX > 0) THEN + IN => INITIALIZATION(LP%INIT_INDEX) IF (ANY(IN%ORIENTATION_RAMP_INDEX > 0)) THEN TEMP_ORIENTATION(1) = EVALUATE_RAMP(T,IN%ORIENTATION_RAMP_INDEX(1)) TEMP_ORIENTATION(2) = EVALUATE_RAMP(T,IN%ORIENTATION_RAMP_INDEX(2)) @@ -4418,44 +4404,21 @@ SUBROUTINE RADIATION_FVM TEMP_ORIENTATION = TEMP_ORIENTATION / & (SQRT(TEMP_ORIENTATION(1)**2+TEMP_ORIENTATION(2)**2+TEMP_ORIENTATION(3)**2) & +TWO_EPSILON_EB) - COS_DL = -DOT_PRODUCT(TEMP_ORIENTATION(1:3),DLANG(1:3,N)) - IF (COS_DL>ORIENTATION_VIEW_ANGLE(LP%ORIENTATION_INDEX)) THEN - COS_DL = -(TEMP_ORIENTATION(1)*DLX(N) + & - TEMP_ORIENTATION(2)*DLY(N) + & - TEMP_ORIENTATION(3)*DLZ(N)) - BR => BOUNDARY_RADIA(LP%BR_INDEX) - IF (LPC%MASSLESS_TARGET) THEN - BR%BAND(IBND)%ILW(N) = COS_DL * IL(BC%IIG,BC%JJG,BC%KKG) * VIEW_ANGLE_AREA(LP%ORIENTATION_INDEX) - IF (N==NEAREST_RADIATION_ANGLE(LP%ORIENTATION_INDEX)) & - BR%IL(IBND) = IL(BC%IIG,BC%JJG,BC%KKG) - ELSE - ! IL_UP does not account for the absorption of radiation within the cell occupied by the particle - BR%BAND(IBND)%ILW(N) = COS_DL * IL_UP(BC%IIG,BC%JJG,BC%KKG) * VIEW_ANGLE_AREA(LP%ORIENTATION_INDEX) - ENDIF - ENDIF - CYCLE PARTICLE_RADIATION_LOOP ENDIF ENDIF - SELECT CASE(LPC%N_ORIENTATION) - CASE(0) - CYCLE PARTICLE_RADIATION_LOOP - CASE(1) - COS_DL = -DOT_PRODUCT(ORIENTATION_VECTOR(1:3,LP%ORIENTATION_INDEX),DLANG(1:3,N)) - IF (COS_DL>ORIENTATION_VIEW_ANGLE(LP%ORIENTATION_INDEX)) THEN - COS_DL = -(ORIENTATION_VECTOR(1,LP%ORIENTATION_INDEX)*DLX(N) + & - ORIENTATION_VECTOR(2,LP%ORIENTATION_INDEX)*DLY(N) + & - ORIENTATION_VECTOR(3,LP%ORIENTATION_INDEX)*DLZ(N)) - BR => BOUNDARY_RADIA(LP%BR_INDEX) - IF (LPC%MASSLESS_TARGET) THEN - BR%BAND(IBND)%ILW(N) = COS_DL * IL(BC%IIG,BC%JJG,BC%KKG) * VIEW_ANGLE_AREA(LP%ORIENTATION_INDEX) - IF (N==NEAREST_RADIATION_ANGLE(LP%ORIENTATION_INDEX)) & - BR%IL(IBND) = IL(BC%IIG,BC%JJG,BC%KKG) - ELSE - ! IL_UP does not account for the absorption of radiation within the cell occupied by the particle - BR%BAND(IBND)%ILW(N) = COS_DL * IL_UP(BC%IIG,BC%JJG,BC%KKG) * VIEW_ANGLE_AREA(LP%ORIENTATION_INDEX) - ENDIF - ENDIF - END SELECT + COS_DLO = -DOT_PRODUCT(TEMP_ORIENTATION(1:3),DLANG(1:3,N)) + IF (COS_DLO > COS_HALF_VIEW_ANGLE(LP%ORIENTATION_INDEX)) THEN + DLO = -(TEMP_ORIENTATION(1)*DLX(N) + TEMP_ORIENTATION(2)*DLY(N) + TEMP_ORIENTATION(3)*DLZ(N)) + BR => BOUNDARY_RADIA(LP%BR_INDEX) + IF (LPC%MASSLESS_TARGET) THEN + BR%BAND(IBND)%ILW(N) = DLO * IL(BC%IIG,BC%JJG,BC%KKG) * VIEW_ANGLE_FACTOR(LP%ORIENTATION_INDEX) + IF (N==NEAREST_RADIATION_ANGLE(LP%ORIENTATION_INDEX)) & + BR%IL(IBND) = IL(BC%IIG,BC%JJG,BC%KKG) + ELSE + ! IL_UP does not account for the absorption of radiation within the cell occupied by the particle + BR%BAND(IBND)%ILW(N) = DLO * IL_UP(BC%IIG,BC%JJG,BC%KKG) * VIEW_ANGLE_FACTOR(LP%ORIENTATION_INDEX) + ENDIF + ENDIF ENDDO PARTICLE_RADIATION_LOOP ENDIF diff --git a/Source/read.f90 b/Source/read.f90 index b92ad467758..c5d0e121502 100644 --- a/Source/read.f90 +++ b/Source/read.f90 @@ -97,9 +97,9 @@ SUBROUTINE READ_DATA(DT) N_ORIENTATION_VECTOR = 0 ALLOCATE(ORIENTATION_VECTOR(3,0:10)) -ALLOCATE(ORIENTATION_VIEW_ANGLE(0:10)) +ALLOCATE(COS_HALF_VIEW_ANGLE(0:10)) ORIENTATION_VECTOR(1:3,0) = (/0._EB,0._EB,-1._EB/) -ORIENTATION_VIEW_ANGLE = 0._EB +COS_HALF_VIEW_ANGLE = 0._EB ! Open the input file @@ -1736,7 +1736,7 @@ SUBROUTINE READ_MISC HUMIDITY,HVAC_LOCAL_PRESSURE,HVAC_MASS_TRANSPORT_CELL_L,HVAC_PRES_RELAX,HVAC_QFAN,IBLANK_SMV,I_MAX_TEMP,& LES_FILTER_TYPE,LEVEL_SET_ELLIPSE,LEVEL_SET_MODE,& MAXIMUM_VISIBILITY,MAX_LEAK_PATHS,MAX_RAMPS,& - MINIMUM_ZONE_VOLUME,MPI_TIMEOUT,NEIGHBOR_SEPARATION_DISTANCE,NORTH_BEARING,& + MINIMUM_ZONE_VOLUME,MPI_TIMEOUT,NEAR_WALL_PARTICLE_INTERPOLATION,NEIGHBOR_SEPARATION_DISTANCE,NORTH_BEARING,& NOISE,NOISE_VELOCITY,NO_PRESSURE_ZONES,NUCLEATION_SITES,ORIGIN_LAT,ORIGIN_LON,& OVERWRITE,PARTICLE_CFL,PARTICLE_CFL_MAX,PARTICLE_CFL_MIN,PERIODIC_TEST,POSITIVE_ERROR_TEST,& POROUS_FLOOR,PR,PROFILING,& @@ -5475,8 +5475,8 @@ SUBROUTINE PROC_REAC_2 ENDIF ELSE HOC_IF ! Heat of combustion not specified, use EPUMO2 or H_F is fuel is listed + SM => SPECIES_MIXTURE(1) LISTED_FUEL_IF: IF (.NOT. LISTED_FUEL ) THEN - SM => SPECIES_MIXTURE(1) IF (RN2%EPUMO2 < 0._EB) RN2%EPUMO2 = 13100000._EB ! J/kg RN2%HOC_COMPLETE = RN2%EPUMO2 * RN2%NU_O2 * SPECIES(O2_INDEX)%MW / SMF%MW IF (RN%N_SIMPLE_CHEMISTRY_REACTIONS==1) THEN @@ -5515,7 +5515,6 @@ SUBROUTINE PROC_REAC_2 ELSE EPUMO2_IF RN%HEAT_OF_COMBUSTION = SMF%H_F_HOC+RN%S*SPECIES_MIXTURE(1)%H_F_HOC - & (1._EB+RN%S)*SPECIES_MIXTURE(RN%PROD_SMIX_INDEX)%H_F_HOC - SM => SPECIES_MIXTURE(1) RN%EPUMO2 = RN%HEAT_OF_COMBUSTION*SMF%MW*RN%NU(RN%FUEL_SMIX_INDEX)/(RN%NU(1)*SM%MW*SM%MASS_FRACTION(O2_INDEX)) IF (SMF%H_F_HOC /= SMF%H_F) THEN REDEFINE_H_F(RN%FUEL_SMIX_INDEX) = .TRUE. @@ -5524,7 +5523,6 @@ SUBROUTINE PROC_REAC_2 IF (RN%N_SIMPLE_CHEMISTRY_REACTIONS==2) THEN RN2%HEAT_OF_COMBUSTION=SMF2%H_F+RN2%S*SPECIES_MIXTURE(1)%H_F_HOC - & (1._EB+RN2%S)*SPECIES_MIXTURE(RN2%PROD_SMIX_INDEX)%H_F_HOC - SM => SPECIES_MIXTURE(1) RN2%EPUMO2 = RN2%HEAT_OF_COMBUSTION*SMF2%MW*RN2%NU(RN2%FUEL_SMIX_INDEX)/& (RN2%NU(1)*SM%MW*SM%MASS_FRACTION(O2_INDEX)) RN%HOC_COMPLETE = RN%HEAT_OF_COMBUSTION + (1._EB+RN%S)*RN2%HEAT_OF_COMBUSTION @@ -5564,6 +5562,15 @@ SUBROUTINE PROC_REAC_2 RN%AIT_EXCLUSION_ZONE(IZ)%DEVC_INDEX,RN%AIT_EXCLUSION_ZONE(IZ)%CTRL_INDEX,IZ) ENDIF ENDDO + + IF (RN%SIMPLE_CHEMISTRY .AND. RN%HOC_COMPLETE > 2E8_EB) THEN + WRITE(MESSAGE,'(A,I0,A)') 'WARNING: The heat of combustion for REACtion ',NR,' exceeds 200,000 kJ/kg.' + IF (MY_RANK==0) WRITE(LU_ERR,'(A)') TRIM(MESSAGE) + ENDIF + IF (RN%SIMPLE_CHEMISTRY .AND. RN%EPUMO2 > 2E7_EB) THEN + WRITE(MESSAGE,'(A,I0,A)') 'WARNING: The EPUMO2 for REACtion ',NR,' exceeds 20,000 kJ/kg.' + IF (MY_RANK==0) WRITE(LU_ERR,'(A)') TRIM(MESSAGE) + ENDIF ENDDO REAC_LOOP @@ -5960,15 +5967,16 @@ SUBROUTINE READ_PART IF (ANY(ABS(ORIENTATION(1:3))>TWO_EPSILON_EB)) LPC%N_ORIENTATION = LPC%N_ORIENTATION + 1 IF (LPC%N_ORIENTATION>0) THEN + ORIENTED_PARTICLES = .TRUE. LPC%INCLUDE_BOUNDARY_RADIA_TYPE = .TRUE. N_ORIENTATION_VECTOR = N_ORIENTATION_VECTOR + 1 LPC%ORIENTATION_INDEX = N_ORIENTATION_VECTOR IF (N_ORIENTATION_VECTOR>UBOUND(ORIENTATION_VECTOR,DIM=2)) THEN ORIENTATION_VECTOR => REALLOCATE2D(ORIENTATION_VECTOR,1,3,0,N_ORIENTATION_VECTOR+10) - ORIENTATION_VIEW_ANGLE => REALLOCATE(ORIENTATION_VIEW_ANGLE,0,N_ORIENTATION_VECTOR+10) + COS_HALF_VIEW_ANGLE => REALLOCATE(COS_HALF_VIEW_ANGLE,0,N_ORIENTATION_VECTOR+10) ENDIF ORIENTATION_VECTOR(1:3,N_ORIENTATION_VECTOR) = ORIENTATION(1:3)/ NORM2(ORIENTATION) - ORIENTATION_VIEW_ANGLE(N_ORIENTATION_VECTOR) = 0._EB + COS_HALF_VIEW_ANGLE(N_ORIENTATION_VECTOR) = 0._EB ENDIF LPC%FREE_AREA_FRACTION = FREE_AREA_FRACTION LPC%POROUS_VOLUME_FRACTION = POROUS_VOLUME_FRACTION @@ -7560,7 +7568,10 @@ SUBROUTINE PROC_MATL DO NS=1,N_TRACKED_SPECIES IF (ML%NU_GAS(NS,1) > 0._EB) THEN ! drr/dT = H_V/(R T_BOIL) and rr = 1 - ML%RENODE_DELTA_T = ML%REAC_RATE_DELTA/(SPECIES(NS)%H_V(INT(ML%TMP_REF(1)))*SPECIES(NS)%MW/(R0*ML%TMP_REF(1)**2)) + ML%RENODE_DELTA_T = ML%H_R(1,INT(ML%TMP_REF(1)))*SPECIES(NS)%MW/R0 + ML%RENODE_DELTA_T = ML%RENODE_DELTA_T * EXP(ML%RENODE_DELTA_T/ML%TMP_REF(1)) * & + EXP(-ML%RENODE_DELTA_T/ML%TMP_REF(1))/ML%TMP_REF(1)**2 + ML%RENODE_DELTA_T = ML%REAC_RATE_DELTA/ML%RENODE_DELTA_T EXIT ENDIF ENDDO @@ -13319,7 +13330,7 @@ SUBROUTINE READ_DEVC IF (ABS(ORIENTATION(1)-ORIENTATION_VECTOR(1,I))UBOUND(ORIENTATION_VECTOR,DIM=2)) THEN ORIENTATION_VECTOR => REALLOCATE2D(ORIENTATION_VECTOR,1,3,0,N_ORIENTATION_VECTOR+10) - ORIENTATION_VIEW_ANGLE => REALLOCATE(ORIENTATION_VIEW_ANGLE,0,N_ORIENTATION_VECTOR+10) + COS_HALF_VIEW_ANGLE => REALLOCATE(COS_HALF_VIEW_ANGLE,0,N_ORIENTATION_VECTOR+10) ENDIF IF (ALL(ABS(ORIENTATION(1:3)) 0) THEN IF (PROPERTY(DV%PROP_INDEX)%VIEW_ANGLE < 180._EB) & - ORIENTATION_VIEW_ANGLE(DV%ORIENTATION_INDEX) = COS(PROPERTY(DV%PROP_INDEX)%VIEW_ANGLE/360._EB * PI) + COS_HALF_VIEW_ANGLE(DV%ORIENTATION_INDEX) = COS(PROPERTY(DV%PROP_INDEX)%VIEW_ANGLE/360._EB * PI) ENDIF ! Create an auto-ignition exclusion zone (AIT) in the cell containing a SPARK diff --git a/Source/type.f90 b/Source/type.f90 index 114f2699bd0..237ea7c01a6 100644 --- a/Source/type.f90 +++ b/Source/type.f90 @@ -388,7 +388,6 @@ MODULE TYPES INTEGER :: BR_INDEX=0 !< Variables devoted to radiation intensities INTEGER :: TAG !< Unique integer identifier for the particle INTEGER :: CLASS_INDEX=0 !< LAGRANGIAN_PARTICLE_CLASS of particle - INTEGER :: INITIALIZATION_INDEX=0 !< Index for INIT that placed the particle INTEGER :: ORIENTATION_INDEX=0 !< Index in the array of all ORIENTATIONs INTEGER :: WALL_INDEX=0 !< If liquid droplet has stuck to a wall, this is the WALL cell index INTEGER :: DUCT_INDEX=0 !< Index of duct diff --git a/Source/wall.f90 b/Source/wall.f90 index c4c3de57662..be71fdea4c5 100644 --- a/Source/wall.f90 +++ b/Source/wall.f90 @@ -1716,16 +1716,16 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX, VOLSUM,KAPSUM,DXF,DXB,HTCF,HTCB,Q_RAD_OUT,Q_RAD_OUT_OLD,Q_CON_F,Q_CON_B,& Q_WATER_F,Q_WATER_B,LAYER_DIVIDE,TMP_GAS_BACK,GEOM_FACTOR,DT_BC_SUB_OLD,& DEL_DOT_Q_SC,Q_DOT_G_PP,Q_DOT_G_PP_NET,Q_DOT_O2_PP,Q_DOT_O2_PP_NET,R_SURF,U_SURF,V_SURF,W_SURF,T_BC_SUB,DT_BC_SUB,& - Q_NET_F,Q_NET_B,TMP_RATIO,KODXF,KODXB,H_S,T_NODE,C_S,H_NODE,VOL,T_BOIL_EFF,& + Q_NET_F,Q_NET_B,TMP_RATIO,KODXF,KODXB,H_S,T_NODE,C_S,H_NODE,VOL,& RADIUS,HTC_LIMIT,CP1,CP2,DENOM,SF_HTC_F,SF_HTC_B,THICKNESS,DT_FO,DDSUM,NODE_RDT(NWP_MAX) REAL(EB), DIMENSION(N_TRACKED_SPECIES) :: M_DOT_G_PP_ADJUST,M_DOT_G_PP_ADJUST_NET,M_DOT_G_PP_ACTUAL,M_DOT_G_PP_ACTUAL_NET -REAL(EB), DIMENSION(MAX_MATERIALS) :: M_DOT_S_PP,M_DOT_S_PP_NET +REAL(EB), DIMENSION(MAX_MATERIALS) :: M_DOT_S_PP,M_DOT_S_PP_NET,T_BOIL_EFF REAL(EB), DIMENSION(MAX_LPC) :: Q_DOT_PART_S,M_DOT_PART_S REAL(EB), DIMENSION(NWP_MAX) :: TMP_S,RHO_H_S REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: RHO_DOT,INT_WGT REAL(EB), DIMENSION(MAX_LAYERS) :: DX_MIN REAL(EB), DIMENSION(MAX_LAYERS,MAX_MATERIALS) :: RHO_ADJUSTED -REAL(EB), DIMENSION(NWP_MAX) :: AAS,BBS,CCS,DDS,DDT,Q_S,TWO_DX_KAPPA_S,DX_S,MF_FRAC,REGRID_FACTOR +REAL(EB), DIMENSION(NWP_MAX) :: AAS,BBS,CCS,DDS,DDT,Q_S,Q_IR,Q_ADD,TWO_DX_KAPPA_S,DX_S,MF_FRAC,REGRID_FACTOR REAL(EB), DIMENSION(0:NWP_MAX+1) :: RHO_S,DELTA_TMP,RDX_S REAL(EB), DIMENSION(0:NWP_MAX) :: X_S_NEW,RDXN_S,R_S,R_S_NEW,DX_WGT_S INTEGER, DIMENSION(0:NWP_MAX+1) :: LAYER_INDEX @@ -2145,6 +2145,7 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX, PYROLYSIS_PREDICTED_IF: IF (ONE_D%PYROLYSIS_MODEL==PYROLYSIS_PREDICTED) THEN + T_BOIL_EFF = TMPA CALL PERFORM_PYROLYSIS ELSEIF (ONE_D%PYROLYSIS_MODEL==PYROLYSIS_SPECIFIED) THEN PYROLYSIS_PREDICTED_IF @@ -2156,10 +2157,15 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX, ENDIF PYROLYSIS_PREDICTED_IF + ! Determine additional heat sources + + Q_ADD = 0._EB + Q_IR = 0._EB + ! 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))) + Q_ADD(I) = ONE_D%HEAT_SOURCE(LAYER_INDEX(I))*EVALUATE_RAMP(T-T_BEGIN,ONE_D%RAMP_IHS_INDEX(LAYER_INDEX(I))) ENDDO ! Add special convection term for Boundary Fuel Model @@ -2179,7 +2185,7 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX, HTCF = HEAT_TRANSFER_COEFFICIENT(NM,DTMP,SF_HTC_F,SF,CFACE_INDEX_IN=CFACE_INDEX) ENDIF DEL_DOT_Q_SC = HTCF*DTMP - Q_S(I) = Q_S(I) + SF%SURFACE_VOLUME_RATIO(LAYER_INDEX(I))*SF%PACKING_RATIO(LAYER_INDEX(I))*DEL_DOT_Q_SC + Q_ADD(I) = Q_ADD(I) + SF%SURFACE_VOLUME_RATIO(LAYER_INDEX(I))*SF%PACKING_RATIO(LAYER_INDEX(I))*DEL_DOT_Q_SC ! Track average h_c for computing h_m in SURFACE_OXIDATION_MODEL B1%HEAT_TRANS_COEF = B1%HEAT_TRANS_COEF + HTCF N = N + 1 @@ -2210,21 +2216,22 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX, RFLUX_UP = B1%Q_RAD_IN + (1._EB-B1%EMISSIVITY)*Q_RAD_OUT_OLD/(B1%EMISSIVITY+1.0E-10_EB) DO I=1,NWP RFLUX_DOWN = ( RFLUX_UP + TWO_DX_KAPPA_S(I)*SIGMA*ONE_D%TMP(I)**4 ) / (1._EB + TWO_DX_KAPPA_S(I)) - Q_S(I) = Q_S(I) + (RFLUX_UP - RFLUX_DOWN)*RDX_S(I) + Q_IR(I) = Q_IR(I) + (RFLUX_UP - RFLUX_DOWN)*RDX_S(I) RFLUX_UP = RFLUX_DOWN ENDDO ! solution outwards RFLUX_UP = Q_RAD_IN_B + (1._EB-E_WALLB)*RFLUX_UP DO I=NWP,1,-1 RFLUX_DOWN = ( RFLUX_UP + TWO_DX_KAPPA_S(I)*SIGMA*ONE_D%TMP(I)**4 ) / (1._EB + TWO_DX_KAPPA_S(I)) - Q_S(I) = Q_S(I) + (RFLUX_UP - RFLUX_DOWN)*RDX_S(I) + Q_IR(I) = Q_IR(I) + (RFLUX_UP - RFLUX_DOWN)*RDX_S(I) RFLUX_UP = RFLUX_DOWN ENDDO Q_RAD_OUT = B1%EMISSIVITY*RFLUX_DOWN ENDIF + ! Add internal radiation and additional heat sources to pyrolysis + Q_S = Q_S + Q_IR + Q_ADD ! If the 3D solver is used, divide Q_S by 3 - Q_S = Q_S/REAL(SF%HT_DIM,EB) ! Explicitly update the temperature field and adjust time step if the change in temperature exceeds DELTA_TMP_MAX @@ -2244,7 +2251,15 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX, DT_BC_SUB_OLD = DT_BC_SUB DT_BC_SUB = DT_BC/REAL(MIN(NINT(SF%TIME_STEP_FACTOR*WALL_INCREMENT),MAX(1,NINT(TMP_RATIO))),EB) DT_BC_SUB = MIN( DT_BC-T_BC_SUB , DT_BC_SUB , DT_FO ) - IF (ONE_D%PYROLYSIS_MODEL==PYROLYSIS_PREDICTED .AND. DT_BC_SUB_OLD/=DT_BC_SUB) CALL PERFORM_PYROLYSIS + ! If DT change, rebuild Q_S + IF (ONE_D%PYROLYSIS_MODEL==PYROLYSIS_PREDICTED .AND. DT_BC_SUB_OLD/=DT_BC_SUB) THEN + Q_S = 0._EB + CALL PERFORM_PYROLYSIS + ! Add internal radiation and additional heat sources to pyrolysis + Q_S = Q_S + Q_IR + Q_ADD + ! If the 3D solver is used, divide Q_S by 3 + Q_S = Q_S/REAL(SF%HT_DIM,EB) + ENDIF ENDIF T_BC_SUB = T_BC_SUB + DT_BC_SUB @@ -2300,7 +2315,7 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX, IF (RHO_DOT(N,I) > TWO_EPSILON_EB) NODE_RDT(I) = MIN(NODE_RDT(I),MATERIAL(ONE_D%MATL_INDEX(N))%RENODE_DELTA_T) REGRID_FACTOR(I) = REGRID_FACTOR(I) + ONE_D%MATL_COMP(N)%RHO(I)/RHO_ADJUSTED(LAYER_INDEX(I),N) ENDDO MATERIAL_LOOP1a - + ! If there is any non-shrinking material, the material matrix will remain, and no shrinking is allowed MATERIAL_LOOP1b: DO N=1,ONE_D%N_MATL @@ -3016,7 +3031,7 @@ SUBROUTINE PYROLYSIS(N_MATS,MATL_INDEX,SURF_INDEX,IIG,JJG,KKG,TMP_S,TMP_F,Y_O2_F REAL(EB), DIMENSION(:), INTENT(OUT) :: M_DOT_G_PPP_ADJUST(N_TRACKED_SPECIES),M_DOT_G_PPP_ACTUAL(N_TRACKED_SPECIES) REAL(EB), DIMENSION(:), INTENT(OUT) :: M_DOT_S_PPP(MAX_MATERIALS),Q_DOT_PART(MAX_LPC),M_DOT_PART(MAX_LPC) REAL(EB), INTENT(OUT) :: Q_DOT_S_PPP,Q_DOT_G_PPP,Q_DOT_O2_PPP,B_NUMBER -REAL(EB), INTENT(INOUT) :: T_BOIL_EFF +REAL(EB), INTENT(INOUT) :: T_BOIL_EFF(MAX_MATERIALS) INTEGER, INTENT(IN), DIMENSION(:) :: MATL_INDEX(N_MATS) INTEGER :: N,NN,NNN,J,NS,SMIX_INDEX(N_MATS),NWP,NP,NP2,ITMP TYPE(MATERIAL_TYPE), POINTER :: ML @@ -3027,7 +3042,7 @@ SUBROUTINE PYROLYSIS(N_MATS,MATL_INDEX,SURF_INDEX,IIG,JJG,KKG,TMP_S,TMP_F,Y_O2_F D_FILM,H_MASS,RE_L,SHERWOOD,MFLUX,MU_FILM,SC_FILM,TMP_FILM,TMP_G,U2,V2,W2,VEL,& DR,R_S_0,R_S_1,H_R,H_R_B,H_S_B,H_S,LENGTH_SCALE,SUM_Y_GAS,SUM_Y_SV,NU_O2_CHAR,Y_O2_S,& SUM_Y_SV_SMIX(N_TRACKED_SPECIES),X_L_SUM,RHO_DOT_EXTRA,MFLUX_MAX,RHO_FILM,CP_FILM,PR_FILM,K_FILM,EVAP_FILM_FAC,& - RHO_DOT,RHO_DOT_REAC(MAX_REACTIONS),RHO_DOT_REAC_SUM + RHO_DOT,RHO_DOT_REAC(MAX_REACTIONS),RHO_DOT_REAC_SUM,H_MASS_DNS LOGICAL :: LIQUID(N_MATS),SPEC_ID_ALREADY_USED(N_MATS),DO_EVAPORATION B_NUMBER = 0._EB @@ -3102,8 +3117,8 @@ SUBROUTINE PYROLYSIS(N_MATS,MATL_INDEX,SURF_INDEX,IIG,JJG,KKG,TMP_S,TMP_F,Y_O2_F ! Determine volume fraction of MATL N in the liquid and then the surface vapor layer - T_BOIL_EFF = ML%TMP_BOIL - CALL GET_EQUIL_DATA(MW(N),TMP_F,PBAR(KKG,PRESSURE_ZONE(IIG,JJG,KKG)),H_R,H_R_B,T_BOIL_EFF,X_SV(N),ML%H_R(1,:)) + T_BOIL_EFF(N) = ML%TMP_BOIL + CALL GET_EQUIL_DATA(MW(N),TMP_F,PBAR(KKG,PRESSURE_ZONE(IIG,JJG,KKG)),H_R,H_R_B,T_BOIL_EFF(N),X_SV(N),ML%H_R(1,:)) X_L(N) = RHO_S(N)/(ML%RHO_S*X_L_SUM) ! Volume fraction of MATL component N in the liquid X_SV(N) = X_L(N)*X_SV(N) ! Volume fraction of MATL component N in the surface vapor based on Raoult's law @@ -3159,37 +3174,44 @@ SUBROUTINE PYROLYSIS(N_MATS,MATL_INDEX,SURF_INDEX,IIG,JJG,KKG,TMP_S,TMP_F,Y_O2_F H_MASS = SF%HM_FIXED - ELSEIF (SIM_MODE==DNS_MODE) THEN H_MASS_IF + ELSE H_MASS_IF SELECT CASE(ABS(IOR)) - CASE(1); H_MASS = 2._EB*D_FILM*RDX(IIG) - CASE(2); H_MASS = 2._EB*D_FILM*RDY(JJG) - CASE(3); H_MASS = 2._EB*D_FILM*RDZ(KKG) + CASE(0); H_MASS_DNS = 0._EB + CASE(1); H_MASS_DNS = 2._EB*D_FILM*RDX(IIG) + CASE(2); H_MASS_DNS = 2._EB*D_FILM*RDY(JJG) + CASE(3); H_MASS_DNS = 2._EB*D_FILM*RDZ(KKG) END SELECT + + IF (SIM_MODE==DNS_MODE) THEN - ELSE H_MASS_IF + H_MASS = H_MASS_DNS - IF (PRESENT(LPU) .AND. PRESENT(LPV) .AND. PRESENT(LPW)) THEN - U2 = 0.5_EB*(U(IIG,JJG,KKG)+U(IIG-1,JJG,KKG)) - V2 = 0.5_EB*(V(IIG,JJG,KKG)+V(IIG,JJG-1,KKG)) - W2 = 0.5_EB*(W(IIG,JJG,KKG)+W(IIG,JJG,KKG-1)) - VEL = SQRT((U2-LPU)**2+(V2-LPV)**2+(W2-LPW)**2) - ELSE - VEL = SQRT(2._EB*KRES(IIG,JJG,KKG)) - ENDIF - CALL GET_VISCOSITY(ZZ_GET,MU_FILM,TMP_FILM) - IF (PRESENT(R_DROP)) THEN - LENGTH_SCALE = 2._EB*R_DROP ELSE - LENGTH_SCALE = SF%CONV_LENGTH + IF (PRESENT(LPU) .AND. PRESENT(LPV) .AND. PRESENT(LPW)) THEN + U2 = 0.5_EB*(U(IIG,JJG,KKG)+U(IIG-1,JJG,KKG)) + V2 = 0.5_EB*(V(IIG,JJG,KKG)+V(IIG,JJG-1,KKG)) + W2 = 0.5_EB*(W(IIG,JJG,KKG)+W(IIG,JJG,KKG-1)) + VEL = SQRT((U2-LPU)**2+(V2-LPV)**2+(W2-LPW)**2) + ELSE + VEL = SQRT(2._EB*KRES(IIG,JJG,KKG)) + ENDIF + CALL GET_VISCOSITY(ZZ_GET,MU_FILM,TMP_FILM) + IF (PRESENT(R_DROP)) THEN + LENGTH_SCALE = 2._EB*R_DROP + ELSE + + LENGTH_SCALE = SF%CONV_LENGTH + ENDIF + RE_L = RHO_FILM*VEL*LENGTH_SCALE/MU_FILM + SELECT CASE(SF%GEOMETRY) + CASE DEFAULT ; SHERWOOD = 0.037_EB*SC_FILM**ONTH*RE_L**0.8_EB + CASE(SURF_SPHERICAL) ; SHERWOOD = 2._EB + 0.6_EB*SC_FILM**ONTH*SQRT(RE_L) + END SELECT + + H_MASS = MAX(H_MASS_DNS,SHERWOOD*D_FILM/LENGTH_SCALE) ENDIF - RE_L = RHO_FILM*VEL*LENGTH_SCALE/MU_FILM - SELECT CASE(SF%GEOMETRY) - CASE DEFAULT ; SHERWOOD = 0.037_EB*SC_FILM**ONTH*RE_L**0.8_EB - CASE(SURF_SPHERICAL) ; SHERWOOD = 2._EB + 0.6_EB*SC_FILM**ONTH*SQRT(RE_L) - END SELECT - H_MASS = SHERWOOD*D_FILM/LENGTH_SCALE ENDIF H_MASS_IF ENDIF IF_DO_EVAPORATION @@ -3238,17 +3260,17 @@ SUBROUTINE PYROLYSIS(N_MATS,MATL_INDEX,SURF_INDEX,IIG,JJG,KKG,TMP_S,TMP_F,Y_O2_F IF (DX_S(SOLID_CELL_INDEX)>TWO_EPSILON_EB) THEN ! If the liquid temperature (TMP_S) is greater than the boiling temperature of the current liquid component - ! (ML%TMP_BOIL), calculate the additional mass loss rate of this component (RHO_DOT_EXTRA) necessary to bring + ! ((T_BOIL_EFF(N)), calculate the additional mass loss rate of this component (RHO_DOT_EXTRA) necessary to bring ! the liquid temperature back to the boiling temperature. RHO_DOT_EXTRA = 0._EB - IF (TMP_S>ML%TMP_BOIL) THEN + IF (TMP_S>T_BOIL_EFF(N)) THEN ITMP = MIN(I_MAX_TEMP,INT(TMP_S)) H_S = ML%H(ITMP) + (TMP_S-REAL(ITMP,EB))*(ML%H(ITMP+1)-ML%H(ITMP)) - ITMP = INT(ML%TMP_BOIL) - H_S = H_S - (ML%H(ITMP) + (ML%TMP_BOIL-REAL(ITMP,EB))*(ML%H(ITMP+1)-ML%H(ITMP))) + ITMP = INT(T_BOIL_EFF(N)) + H_S = H_S - (ML%H(ITMP) + (T_BOIL_EFF(N)-REAL(ITMP,EB))*(ML%H(ITMP+1)-ML%H(ITMP))) H_S = H_S * RHO_S(N) - H_R = ML%H_R(1,NINT(ML%TMP_BOIL)) + H_R = ML%H_R(1,NINT(T_BOIL_EFF(N))) RHO_DOT_EXTRA = H_S/(H_R*DT_BC) ! kg/m3/s ENDIF diff --git a/Utilities/Matlab/FDS_verification_dataplot_inputs.csv b/Utilities/Matlab/FDS_verification_dataplot_inputs.csv index a3a6c31aafa..08cb656bd3e 100644 --- a/Utilities/Matlab/FDS_verification_dataplot_inputs.csv +++ b/Utilities/Matlab/FDS_verification_dataplot_inputs.csv @@ -405,8 +405,8 @@ d,multiple_reac_hrrpua,Species/multiple_reac_hrrpua_git.txt,Species/multiple_rea d,multiple_reac_n_simple,Species/multiple_reac_n_simple_git.txt,Species/multiple_reac_n_simple.csv,1,2,Time,CH4_CO|CH4_H2,Ideal CO|Ideal H2,ko|ro,0,100000,,0,100000,-1.00E+09,1.00E+09,0,Species/multiple_reac_n_simple_devc.csv,2,3,Time,CH4_CO|CH4_H2,FDS CO|FDS H2,k|r,0,100000,,0,100000,-1.00E+09,1.00E+09,0,CH4 Species Mass,Time (s),Mass (kg),0,0.0001,1,0,0.006,1,no,0.03 0.90,East,,1,linear,FDS_User_Guide/SCRIPT_FIGURES/multiple_reac_n_simple_CH4,Relative Error,end,0.01,Species,yd,y,TeX d,multiple_reac_n_simple,Species/multiple_reac_n_simple_git.txt,Species/multiple_reac_n_simple.csv,1,2,Time,C3H8_CO|C3H8_H2O,Ideal CO|Ideal H2O,ko|ro,0,100000,,0,100000,-1.00E+09,1.00E+09,0,Species/multiple_reac_n_simple_devc.csv,2,3,Time,C3H8_CO|C3H8_H2O,FDS CO|FDS H2O,k|r,0,100000,,0,100000,-1.00E+09,1.00E+09,0,C3H8 Species Mass,Time (s),Mass (kg),0,0.0001,1,0,0.025,1,no,0.03 0.90,East,,1,linear,FDS_User_Guide/SCRIPT_FIGURES/multiple_reac_n_simple_C3H8,Relative Error,end,0.01,Species,ys,y,TeX d,multiple_reac_n_simple,Species/multiple_reac_n_simple_git.txt,Species/multiple_reac_n_simple.csv,1,2,Time,C2H6_CO|C2H6_H2,Ideal CO|Ideal H2O,ko|ro,0,100000,,0,100000,-1.00E+09,1.00E+09,0,Species/multiple_reac_n_simple_devc.csv,2,3,Time,C2H6_CO|C2H6_H2,FDS CO|FDS H2,k|r,0,100000,,0,100000,-1.00E+09,1.00E+09,0,C2H6 Species Mass,Time (s),Mass (kg),0,0.0001,1,0,0.035,1,no,0.03 0.90,East,,1,linear,FDS_User_Guide/SCRIPT_FIGURES/multiple_reac_n_simple_C2H6,Relative Error,end,0.01,Species,yd,y,TeX -d,methanol_evaporation,Pyrolysis/methanol_evaporation_git.txt,Pyrolysis/methanol_evaporation_devc.csv,2,3,Time,mdot,Computed Evaporation Rate (mdot),k-,0,100000,,0,100000,-1.00E+09,1.00E+09,0,Pyrolysis/methanol_evaporation_devc.csv,2,3,Time,mdot2,Ideal Evaporation Rate (mdot2),k--,0,100000,,0,100000,-1.00E+09,1.00E+09,0,Liquid Evaporation (methanol\_evaporation),Time (min),Mass Loss Rate (kg/m²/s),0,15,60,0,0.02,1,no,0.05 0.90,SouthEast,,1,linear,FDS_User_Guide/SCRIPT_FIGURES/methanol_evaporation_mdot,Relative Error,end,0.015,Pyrolysis,mx,m,TeX -d,methanol_evaporation,Pyrolysis/methanol_evaporation_git.txt,Pyrolysis/methanol_evaporation.csv,1,2,Time,Tb,Measured Boiling Temperature (Tb),ko,0,100000,,0,100000,-1.00E+09,1.00E+09,0,Pyrolysis/methanol_evaporation_devc.csv,2,3,Time,Tsurf,Surface Temperature (Tsurf),k-,0,100000,,0,100000,-1.00E+09,1.00E+09,0,Liquid Evaporation (methanol\_evaporation),Time (min),Temperature (°C),0,15,60,0,100,1,no,0.05 0.90,SouthEast,,1,linear,FDS_User_Guide/SCRIPT_FIGURES/methanol_evaporation_temp,Relative Error,end,0.015,Pyrolysis,mx,m,TeX +d,methanol_evaporation,Pyrolysis/methanol_evaporation_git.txt,Pyrolysis/methanol_evaporation_devc.csv,2,3,Time,mdot,Computed Evaporation Rate (mdot),k-,0,100000,,0,100000,-1.00E+09,1.00E+09,0,Pyrolysis/methanol_evaporation_devc.csv,2,3,Time,mdot2,Ideal Evaporation Rate (mdot2),k--,0,100000,,0,100000,-1.00E+09,1.00E+09,0,Liquid Evaporation (methanol\_evaporation),Time (min),Mass Loss Rate (kg/m²/s),0,6,60,0,0.02,1,no,0.05 0.90,SouthEast,,1,linear,FDS_User_Guide/SCRIPT_FIGURES/methanol_evaporation_mdot,Relative Error,end,0.02,Pyrolysis,mx,m,TeX +d,methanol_evaporation,Pyrolysis/methanol_evaporation_git.txt,Pyrolysis/methanol_evaporation.csv,1,2,Time,Tb,Measured Boiling Temperature (Tb),ko,0,100000,,0,100000,-1.00E+09,1.00E+09,0,Pyrolysis/methanol_evaporation_devc.csv,2,3,Time,Tsurf,Surface Temperature (Tsurf),k-,0,100000,,0,100000,-1.00E+09,1.00E+09,0,Liquid Evaporation (methanol\_evaporation),Time (min),Temperature (°C),0,6,60,0,100,1,no,0.05 0.90,SouthEast,,1,linear,FDS_User_Guide/SCRIPT_FIGURES/methanol_evaporation_temp,Relative Error,end,0.015,Pyrolysis,mx,m,TeX d,MO_velocity_profile_stable,Atmospheric_Effects/MO_velocity_profile_stable_git.txt,Atmospheric_Effects/MO_velocity_profile_stable.csv,1,2,z (m),u (m/s),MO profile,k,0,100000,,0,100000,-1.00E+09,1.00E+09,0,Atmospheric_Effects/MO_velocity_profile_stable_line.csv,2,3,z,u,FDS profile,k--,0,100000,,0,100000,-1.00E+09,1.00E+09,0,Monin-Obukhov profile stable,z (m),u (m/s),0,32,1,0,10,1,no,0.03 0.90,SouthEast,,1,linear,FDS_Verification_Guide/SCRIPT_FIGURES/MO_velocity_profile_stable,Relative Error,area,0.05,Flowfields,r>,r,TeX d,MO_velocity_profile_unstable,Atmospheric_Effects/MO_velocity_profile_unstable_git.txt,Atmospheric_Effects/MO_velocity_profile_unstable.csv,1,2,z (m),u (m/s),MO profile,k,0,100000,,0,100000,-1.00E+09,1.00E+09,0,Atmospheric_Effects/MO_velocity_profile_unstable_line.csv,2,3,z,u,FDS profile,k--,0,100000,,0,100000,-1.00E+09,1.00E+09,0,Monin-Obukhov profile unstable,z (m),u (m/s),0,32,1,0,15,1,no,0.03 0.90,SouthEast,,1,linear,FDS_Verification_Guide/SCRIPT_FIGURES/MO_velocity_profile_unstable,Relative Error,area,0.05,Flowfields,r>,r,TeX d,Morvan_TGA,WUI/Morvan_TGA_git.txt,WUI/Morvan_Data_Mass.csv,1,2,T (C),Normalized Mass (M/M0),Exp (Morvan 2004),k^,0,100000,,0,100000,-1.00E+09,1.00E+09,0,WUI/Morvan_TGA_tga.csv,2,3,Temp,Total Mass,FDS TGA (Total Mass),k-,0,100000,,0,100000,-1.00E+09,1.00E+09,0,Morvan TGA; 1.6 °C/min,Temperature (°C),Normalized Mass,0,700,1,0,1.2,1,no,0.05 0.90,East,,1,linear,FDS_Verification_Guide/SCRIPT_FIGURES/Morvan_TGA_Total_Mass,N/A,mean,0,Needle TGA,kd,k,TeX @@ -451,7 +451,7 @@ d,obst_coarse_fine_interface,Pressure_Effects/obst_coarse_fine_interface_git.txt d,opening_ulmat,Pressure_Solver/opening_ulmat_git.txt,Pressure_Solver/opening_pressure_error.csv,1,2,Time,Pressure Tolerance,Ideal (Pressure Tolerance),ko--,0,100000,,0,100000,-1.00E+09,1.00E+09,0,Pressure_Solver/opening_ulmat_devc.csv,2,3,Time,perr-max,FDS (p err max),k-,0,100000,,0,100000,-1.00E+09,1.00E+09,0,Pressure Error (opening\_ulmat),Time (s),Pressure Error (Pa),0,10,1,0,1.00E-06,1,no,0.05 0.90,SouthEast,,1,semilogy,FDS_User_Guide/SCRIPT_FIGURES/opening_ulmat,Absolute Error,tolerance,1.00E-10,Pressure Solver,k+,k,TeX d,parabolic_profile,Flowfields/parabolic_profile_git.txt,Flowfields/parabolic_profile.csv,1,2,Time,Pressure,Exact (Pressure),ko,0,100000,,0,100000,-1.00E+09,1.00E+09,0,Flowfields/parabolic_profile_devc.csv,2,3,Time,pres,FDS (pres),k-,0,100000,,0,100000,-1.00E+09,1.00E+09,0,Pressure (parabolic\_profile),Time (s),Pressure (Pa),0,60,1,0,2500,1,no,0.05 0.90,SouthEast,,1,linear,FDS_User_Guide/SCRIPT_FIGURES/parabolic_profile,Relative Error,end,0.01,Pressure Effects,k+,k,TeX d,part_attenuation,Radiation/part_attenuation_git.txt,Radiation/part_attenuation_devc.csv,2,3,Time,Ref,Reference,ko,0,100000,,0,100000,-1.00E+09,1.00E+09,0,Radiation/part_attenuation_devc.csv,2,3,Time,Transparent|Water|Fuel|Opaque,Transparent|Water|Fuel|Opaque,k-|k--|r-.|b-,0,100000,,0,100000,-1.00E+09,1.00E+09,0,Radiation attenuation (part\_attenuation),Time (s),Radiative heat flux (kW/m²),0,2,1,3.9,6,1,no,0.05 0.90,SouthEast,,1,linear,FDS_Verification_Guide/SCRIPT_FIGURES/part_attenuation,N/A,end,0,Radiation,kd,k,TeX -d,part_drag_stretched,WUI/part_drag_stretched_git.txt,WUI/part_drag_stretched_devc.csv,2,3,Time,Fx_p,Fx part,ko,0,100000,,0,100000,-1.00E+09,1.00E+09,0,WUI/part_drag_stretched_devc.csv,2,3,Time,Fx_g,Fx gas,k-,0,100000,,0,100000,-1.00E+09,1.00E+09,0,Particle drag (part\_drag\_stretched),Time (s),Drag force (N),0,30,1,-2e-5,0e-5,1,no,0.05 0.90,SouthEast,,1,linear,FDS_Verification_Guide/SCRIPT_FIGURES/part_drag_stretched,Relative Error,end,0.03,WUI,kd,k,TeX +d,part_drag_stretched,WUI/part_drag_stretched_git.txt,WUI/part_drag_stretched_devc.csv,2,3,Time,Fx_p,Fx part,ko,0,100000,,0,100000,-1.00E+09,1.00E+09,0,WUI/part_drag_stretched_devc.csv,2,3,Time,Fx_g,Fx gas,k-,0,100000,,0,100000,-1.00E+09,1.00E+09,0,Particle drag (part\_drag\_stretched),Time (s),Drag force (N),0,30,1,-2.00E-02,0.00E+00,1,no,0.05 0.90,SouthEast,,1,linear,FDS_Verification_Guide/SCRIPT_FIGURES/part_drag_stretched,Relative Error,end,0.03,WUI,kd,k,TeX d,particle_drag_U10_N16,Sprinklers_and_Sprays/particle_drag_U10_N16_git.txt,Sprinklers_and_Sprays/particle_drag_U10_N16.csv,1,2,T,U,Analytical,ko,0,100000,,0,100000,-1.00E+09,1.00E+09,0,Sprinklers_and_Sprays/particle_drag_U10_N16_devc.csv,2,3,Time,U-VEL,FDS,k-,0,100000,,0,100000,-1.00E+09,1.00E+09,0,Gas phase velocity (particle\_drag\_A),Time (s),Velocity (m/s),0,100,1,0,12,1,no,0.05 0.90,East,,1,linear,FDS_Verification_Guide/SCRIPT_FIGURES/particle_drag_A,Absolute Error,end,0.01,Sprinklers and Sprays,kd,k,TeX d,particle_drag_U50_N16,Sprinklers_and_Sprays/particle_drag_U50_N16_git.txt,Sprinklers_and_Sprays/particle_drag_U50_N16.csv,1,2,T,U,Analytical,ko,0,100000,,0,100000,-1.00E+09,1.00E+09,0,Sprinklers_and_Sprays/particle_drag_U50_N16_devc.csv,2,3,Time,U-VEL,FDS,k-,0,100000,,0,100000,-1.00E+09,1.00E+09,0,Gas phase velocity (particle\_drag\_B),Time (s),Velocity (m/s),0,50,1,0,60,1,no,0.05 0.90,East,,1,linear,FDS_Verification_Guide/SCRIPT_FIGURES/particle_drag_B,Absolute Error,end,0.01,Sprinklers and Sprays,kd,k,TeX d,particle_drag_U100_N16,Sprinklers_and_Sprays/particle_drag_U100_N16_git.txt,Sprinklers_and_Sprays/particle_drag_U100_N16.csv,1,2,T,U,Analytical,ko,0,100000,,0,100000,-1.00E+09,1.00E+09,0,Sprinklers_and_Sprays/particle_drag_U100_N16_devc.csv,2,3,Time,U-VEL,FDS,k-,0,100000,,0,100000,-1.00E+09,1.00E+09,0,Gas phase velocity (particle\_drag\_C),Time (s),Velocity (m/s),0,25,1,0,120,1,no,0.05 0.90,East,,1,linear,FDS_Verification_Guide/SCRIPT_FIGURES/particle_drag_C,Absolute Error,end,0.01,Sprinklers and Sprays,kd,k,TeX diff --git a/Utilities/Matlab/FDS_verification_script.m b/Utilities/Matlab/FDS_verification_script.m index 0e1b682ef6f..6a272d843a7 100644 --- a/Utilities/Matlab/FDS_verification_script.m +++ b/Utilities/Matlab/FDS_verification_script.m @@ -107,5 +107,6 @@ disp('nat_conv_hot_plate...'); nat_conv_hot_plate disp('tree_shapes...'); tree_shapes disp('impinging_jet...'); impinging_jet +disp('part_drag_profile...'); part_drag_profile display('verification scripts completed successfully!') diff --git a/Utilities/Matlab/scripts/part_drag_profile.m b/Utilities/Matlab/scripts/part_drag_profile.m new file mode 100644 index 00000000000..59d4e48454b --- /dev/null +++ b/Utilities/Matlab/scripts/part_drag_profile.m @@ -0,0 +1,101 @@ +% McDermott +% 8-27-24 +% part_drag_profile.m + +close all +clear all + +plot_style + +figure +set(gcf,'Visible',Figure_Visibility); +set(gca,'Units',Plot_Units) +set(gca,'Position',[Plot_X Plot_Y Plot_Width Plot_Height]) + +% prescribed velocity profile + +% mpv = 4.; % kg/m^3, mass_per_volume (from FDS input file) +% v_xb = 10^3; % volume of XB region on init line in FDS input file +% nppc = 10; % number of particles per cell +% n = 5*5*20*nppc; % number of particles +% rho_p = 400; % density of grass, kg/m^3 +r_p = 0.001; % radius, m +l_p = 0.02; % length, m +v_p = pi*(r_p)^2*l_p; % volume of a single particle, m^3 +shape_factor = 0.25; % assumes random orientation of cylinders +a_p = shape_factor*l_p*(2*pi*r_p); % projected area, m^2 +% m_p = rho_p*v_p; % mass of single particle, kg +% pwt = mpv*v_xb/(n*m_p) % particle weight factor + +z = linspace(0,10,20); +u_z = z; +c_d = 2.8; % from FDS input file (specified) +rho_g = 1.195; % from FDS out file +f_x = c_d * a_p * 0.5*rho_g*(u_z.^2); % drag experienced by a single particle + +H(1)=plot(z,-f_x,'k-'); hold on + +set(gca,'FontName',Font_Name) +set(gca,'FontSize',Label_Font_Size) + +ddir='../../Verification/WUI/'; +chid={'part_drag_prof_ux','part_drag_prof_uy','part_drag_prof_uz',... + 'part_drag_prof_vx','part_drag_prof_vy','part_drag_prof_vz',... + 'part_drag_prof_wx','part_drag_prof_wy','part_drag_prof_wz'}; +j={1,2,3,1,2,3,1,2,3}; % coordinate direction (x=1, y=2, z=3) + +for i=1:length(chid) % chid_for + + skip_case = 0; + if ~exist([ddir,chid{i},'_1.prt5']) + display(['Error: File ' [ddir,chid{i},'_1.prt5'] ' does not exist. Skipping case.']) + skip_case = 1; + end + + if skip_case + return + end + + [STIME, XP, YP, ZP, QP] = read_prt5([ddir,chid{i},'_1.prt5'],'real*4'); + + switch j{i} + case 1 + H(2)=plot(XP(end,:),QP(end,:,1,1)./QP(end,:,1,2),'b.'); + v = abs( c_d * a_p * 0.5*rho_g*(XP(end,:).^2) - QP(end,:,1,1)./QP(end,:,1,2) ); + case 2 + H(2)=plot(YP(end,:),QP(end,:,1,1)./QP(end,:,1,2),'b.'); + v = abs( c_d * a_p * 0.5*rho_g*(YP(end,:).^2) - QP(end,:,1,1)./QP(end,:,1,2) ); + case 3 + H(2)=plot(ZP(end,:),QP(end,:,1,1)./QP(end,:,1,2),'b.'); + v = abs( c_d * a_p * 0.5*rho_g*(ZP(end,:).^2) - QP(end,:,1,1)./QP(end,:,1,2) ); + end + + err = norm(v)/length(v); + if err>1e-4 + display(['Error: Case ' [ddir,chid{i}] ' error = ' num2str(err)]) + end + +end % chid_for + +xlabel('Position (m)','FontSize',Label_Font_Size) +ylabel('Drag Force (N)','FontSize',Label_Font_Size) +lh=legend(H,'exact','FDS part'); +set(lh,'FontName',Font_Name,'FontSize',Key_Font_Size) + +Git_Filename = [ddir,chid{1},'_git.txt']; +addverstr(gca,Git_Filename,'linear') + +set(gcf,'Visible',Figure_Visibility); +set(gcf,'Units',Paper_Units); +set(gcf,'PaperSize',[Paper_Width Paper_Height]); +set(gcf,'Position',[0 0 Paper_Width Paper_Height]); +print(gcf,'-dpdf','../../Manuals/FDS_Verification_Guide/SCRIPT_FIGURES/part_drag_profile'); + + + + + + + + + diff --git a/Utilities/Scripts/qfds.sh b/Utilities/Scripts/qfds.sh index c77ab87b872..33bd3ad42d1 100755 --- a/Utilities/Scripts/qfds.sh +++ b/Utilities/Scripts/qfds.sh @@ -52,9 +52,10 @@ function usage { echo "Other options:" echo " -b email_address - send an email to email_address when jobs starts, aborts and finishes" echo " -d dir - specify directory where the case is found [default: .]" + echo " -G use GNU OpenMPI version of fds" echo " -I use Intel MPI version of fds" echo " -j prefix - specify a job prefix" - echo " -L use Open MPI version of fds" + echo " -L use Intel Fortran with Open MPI version of fds" echo " -n n - number of MPI processes per node [default: 1]" echo " -P use PBS/Torque" echo " -s - stop job" @@ -124,6 +125,7 @@ n_openmp_threads=1 use_debug= use_devel= use_intel_mpi=1 +use_gnu_openmpi= EMAIL= casedir= use_default_casedir= @@ -143,7 +145,7 @@ commandline=`echo $* | sed 's/-V//' | sed 's/-v//'` #*** read in parameters from command line -while getopts 'b:d:e:hHIj:Ln:o:Pp:q:stT:U:vw:y:Y' OPTION +while getopts 'b:d:e:GhHIj:Ln:o:Pp:q:stT:U:vw:y:Y' OPTION do case $OPTION in b) @@ -155,6 +157,9 @@ case $OPTION in e) exe="$OPTARG" ;; + G) + use_gnu_openmpi=1 + ;; h) usage exit @@ -265,6 +270,9 @@ if [ "$use_intel_mpi" == "1" ]; then exe=$FDSROOT/fds/Build/impi_intel_linux_openmp$DB/fds_impi_intel_linux_openmp$DB fi fi +if [ "$use_gnu_openmpi" == "1" ]; then + exe=$FDSROOT/fds/Build/ompi_gnu_linux$DB/fds_ompi_gnu_linux$DB +fi if [ "$exe" == "" ]; then exe=$FDSROOT/fds/Build/ompi_intel_linux$DB/fds_ompi_intel_linux$DB fi @@ -371,12 +379,12 @@ cat << EOF >> $scriptfile #SBATCH --partition=$queue #SBATCH --ntasks=$n_mpi_processes #SBATCH --cpus-per-task=$n_openmp_threads +#SBATCH --nodes=$nodes #SBATCH --time=$walltime EOF if [[ $n_openmp_threads -gt 1 ]] || [[ $max_mpi_processes_per_node -lt 1000 ]] ; then cat << EOF >> $scriptfile -#SBATCH --nodes=$nodes #SBATCH --ntasks-per-node=$n_mpi_processes_per_node EOF fi diff --git a/Verification/Pyrolysis/methanol_evaporation.csv b/Verification/Pyrolysis/methanol_evaporation.csv index 16a2e699af3..00a1576c901 100644 --- a/Verification/Pyrolysis/methanol_evaporation.csv +++ b/Verification/Pyrolysis/methanol_evaporation.csv @@ -1,11 +1,8 @@ Time,Tb 0,64.65 -100,64.65 -200,64.65 +60,64.65 +120,64.65 +180,64.65 +240,64.65 300,64.65 -400,64.65 -500,64.65 -600,64.65 -700,64.65 -800,64.65 -900,64.65 +360,64.65 diff --git a/Verification/Pyrolysis/methanol_evaporation.fds b/Verification/Pyrolysis/methanol_evaporation.fds index 510840aff06..85f2694dc7f 100644 --- a/Verification/Pyrolysis/methanol_evaporation.fds +++ b/Verification/Pyrolysis/methanol_evaporation.fds @@ -3,7 +3,7 @@ &MESH IJK=12,12,12, XB=-0.6,0,-0.6,0,0,0.6, MULT_ID='m1'/ &MULT ID='m1', DX=0.6, DY=0.6, DZ=0.6, I_UPPER=1, J_UPPER=1, K_UPPER=1 / 8 meshes -&TIME T_END=900. / +&TIME T_END=360. / &DUMP FLUSH_FILE_BUFFERS=T, DT_PROF=5., DT_DEVC=5., DT_HRR=5. / @@ -14,7 +14,7 @@ NU_SPEC = 1. SPEC_ID = 'METHANOL' HEAT_OF_REACTION = 1099 - CONDUCTIVITY = 0.2 + CONDUCTIVITY = 100 SPECIFIC_HEAT = 2.48 DENSITY = 796 ABSORPTION_COEFFICIENT = 1500 @@ -24,7 +24,8 @@ EMISSIVITY = 1. COLOR = 'YELLOW' MATL_ID = 'METHANOL LIQUID' - THICKNESS = 0.1 + THICKNESS = 0.05 + BACKING = 'INSULATED' EXTERNAL_FLUX=20 / &MATL ID = 'STEEL' @@ -63,7 +64,7 @@ &DEVC XB=-0.50,0.50,-0.50,0.50,0.05,0.05, QUANTITY='TOTAL HEAT FLUX', SPATIAL_STATISTIC='MEAN', ID='qdot' / &DEVC XB=-0.50,0.50,-0.50,0.50,0.05,0.05, QUANTITY='MASS FLUX', SPEC_ID='METHANOL', SPATIAL_STATISTIC='MEAN', ID='mdot' / -&DEVC XB=-0.50,0.50,-0.50,0.50,0.05,0.05, QUANTITY='TOTAL HEAT FLUX', SPATIAL_STATISTIC='MEAN', ID='mdot2', CONVERSION_FACTOR=0.000826 / +&DEVC XB=-0.50,0.50,-0.50,0.50,0.05,0.05, QUANTITY='TOTAL HEAT FLUX', SPATIAL_STATISTIC='MEAN', ID='mdot2', CONVERSION_FACTOR=0.00091 / &DEVC XB=-0.50,0.50,-0.50,0.50,0.05,0.05, QUANTITY='WALL TEMPERATURE', SPATIAL_STATISTIC='MEAN', SURF_ID='METHANOL POOL', ID='Tsurf' / &DEVC XB=-0.50,0.50,-0.50,0.50,0.05,0.05, QUANTITY='WALL TEMPERATURE', SPATIAL_STATISTIC='MAX', SURF_ID='METHANOL POOL', ID='Tsurf_max' / diff --git a/Verification/WUI/part_drag_prof_ux.fds b/Verification/WUI/part_drag_prof_ux.fds index 8f33ec45781..75ce580b55f 100644 --- a/Verification/WUI/part_drag_prof_ux.fds +++ b/Verification/WUI/part_drag_prof_ux.fds @@ -4,7 +4,7 @@ &TIME T_END=1.E-6 / -&MISC FREEZE_VELOCITY=T, STRATIFICATION=F, NOISE=F, RAMP_UX='linear ramp'/ +&MISC FREEZE_VELOCITY=T, STRATIFICATION=F, NOISE=F, RAMP_UX='linear ramp', NEAR_WALL_PARTICLE_INTERPOLATION=T/ &SURF ID='right', NO_SLIP=T, VEL=-10, COLOR='RED'/ &SURF ID='left', NO_SLIP=T, COLOR='BLUE'/ diff --git a/Verification/WUI/part_drag_prof_uy.fds b/Verification/WUI/part_drag_prof_uy.fds index 8ab41efaca9..01475106ec1 100644 --- a/Verification/WUI/part_drag_prof_uy.fds +++ b/Verification/WUI/part_drag_prof_uy.fds @@ -4,7 +4,7 @@ &TIME T_END=10 / -&MISC FREEZE_VELOCITY=T, STRATIFICATION=F, NOISE=F, RAMP_UY='linear ramp'/ +&MISC FREEZE_VELOCITY=T, STRATIFICATION=F, NOISE=F, RAMP_UY='linear ramp', NEAR_WALL_PARTICLE_INTERPOLATION=T/ &SURF ID='far', NO_SLIP=T, VEL_T(1)=10, COLOR='RED'/ &SURF ID='near', NO_SLIP=T, COLOR='BLUE'/ diff --git a/Verification/WUI/part_drag_prof_uz.fds b/Verification/WUI/part_drag_prof_uz.fds index c1de2b290c0..574200eacca 100644 --- a/Verification/WUI/part_drag_prof_uz.fds +++ b/Verification/WUI/part_drag_prof_uz.fds @@ -4,7 +4,7 @@ &TIME T_END=10. / -&MISC FREEZE_VELOCITY=T, STRATIFICATION=F, NOISE=F, RAMP_UZ='linear ramp'/ +&MISC FREEZE_VELOCITY=T, STRATIFICATION=F, NOISE=F, RAMP_UZ='linear ramp', NEAR_WALL_PARTICLE_INTERPOLATION=T/ ! &WIND SPEED=1., RAMP_SPEED_Z='linear ramp', DIRECTION=-90/ equivalent to RAMP_UZ, but less general diff --git a/Verification/WUI/part_drag_prof_vx.fds b/Verification/WUI/part_drag_prof_vx.fds index 2b8850e5596..68e400612a9 100644 --- a/Verification/WUI/part_drag_prof_vx.fds +++ b/Verification/WUI/part_drag_prof_vx.fds @@ -4,7 +4,7 @@ &TIME T_END=10. / -&MISC FREEZE_VELOCITY=T, STRATIFICATION=F, NOISE=F, RAMP_VX='linear ramp'/ +&MISC FREEZE_VELOCITY=T, STRATIFICATION=F, NOISE=F, RAMP_VX='linear ramp', NEAR_WALL_PARTICLE_INTERPOLATION=T/ &SURF ID='right', NO_SLIP=T, VEL_T(1)=10, COLOR='RED'/ &SURF ID='left', NO_SLIP=T, COLOR='BLUE'/ diff --git a/Verification/WUI/part_drag_prof_vy.fds b/Verification/WUI/part_drag_prof_vy.fds index ebd053317ae..0af799faa12 100644 --- a/Verification/WUI/part_drag_prof_vy.fds +++ b/Verification/WUI/part_drag_prof_vy.fds @@ -4,7 +4,7 @@ &TIME T_END=1.E-6 / -&MISC FREEZE_VELOCITY=T, STRATIFICATION=F, NOISE=F, RAMP_VY='linear ramp'/ +&MISC FREEZE_VELOCITY=T, STRATIFICATION=F, NOISE=F, RAMP_VY='linear ramp', NEAR_WALL_PARTICLE_INTERPOLATION=T/ &SURF ID='far', NO_SLIP=T, VEL=-10, COLOR='RED'/ &SURF ID='near', NO_SLIP=T, COLOR='BLUE'/ diff --git a/Verification/WUI/part_drag_prof_vz.fds b/Verification/WUI/part_drag_prof_vz.fds index 6b7b4a711ee..3cdcd5e29f3 100644 --- a/Verification/WUI/part_drag_prof_vz.fds +++ b/Verification/WUI/part_drag_prof_vz.fds @@ -4,10 +4,10 @@ &TIME T_END=10. / -&MISC FREEZE_VELOCITY=T, STRATIFICATION=F, NOISE=F, RAMP_VZ='linear ramp'/ +&MISC FREEZE_VELOCITY=T, STRATIFICATION=F, NOISE=F, RAMP_VZ='linear ramp', NEAR_WALL_PARTICLE_INTERPOLATION=T/ -&SURF ID='top', NO_SLIP=T, VEL_T(2)=10, COLOR='BLUE'/ -&SURF ID='bot', NO_SLIP=T, COLOR='RED'/ +&SURF ID='top', NO_SLIP=T, VEL_T(2)=10, COLOR='RED'/ +&SURF ID='bot', NO_SLIP=T, COLOR='BLUE'/ &VENT DB='XMIN', SURF_ID='PERIODIC' / &VENT DB='XMAX', SURF_ID='PERIODIC' / diff --git a/Verification/WUI/part_drag_prof_wx.fds b/Verification/WUI/part_drag_prof_wx.fds index 136483ccb73..7f7e28f84f6 100644 --- a/Verification/WUI/part_drag_prof_wx.fds +++ b/Verification/WUI/part_drag_prof_wx.fds @@ -4,7 +4,7 @@ &TIME T_END=10. / -&MISC FREEZE_VELOCITY=T, STRATIFICATION=F, NOISE=F, RAMP_WX='linear ramp'/ +&MISC FREEZE_VELOCITY=T, STRATIFICATION=F, NOISE=F, RAMP_WX='linear ramp', NEAR_WALL_PARTICLE_INTERPOLATION=T/ &SURF ID='left', NO_SLIP=T, COLOR='BLUE'/ &SURF ID='right', NO_SLIP=T, COLOR='RED', VEL_T(2)=10/ diff --git a/Verification/WUI/part_drag_prof_wy.fds b/Verification/WUI/part_drag_prof_wy.fds index 73fdc240ac9..25632cc06ec 100644 --- a/Verification/WUI/part_drag_prof_wy.fds +++ b/Verification/WUI/part_drag_prof_wy.fds @@ -4,10 +4,10 @@ &TIME T_END=10. / -&MISC FREEZE_VELOCITY=T, STRATIFICATION=F, NOISE=F, RAMP_WY='linear ramp'/ +&MISC FREEZE_VELOCITY=T, STRATIFICATION=F, NOISE=F, RAMP_WY='linear ramp', NEAR_WALL_PARTICLE_INTERPOLATION=T/ +&SURF ID='far', NO_SLIP=T, VEL_T(2)=10, COLOR='RED'/ &SURF ID='near', NO_SLIP=T, COLOR='BLUE'/ -&SURF ID='far', NO_SLIP=T, COLOR='RED', VEL_T(2)=10/ &VENT DB='XMIN', SURF_ID='PERIODIC' / &VENT DB='XMAX', SURF_ID='PERIODIC' / diff --git a/Verification/WUI/part_drag_prof_wz.fds b/Verification/WUI/part_drag_prof_wz.fds index a5bee040669..5a06e860265 100644 --- a/Verification/WUI/part_drag_prof_wz.fds +++ b/Verification/WUI/part_drag_prof_wz.fds @@ -4,7 +4,7 @@ &TIME T_END=1.E-6 / -&MISC FREEZE_VELOCITY=T, STRATIFICATION=F, NOISE=F, RAMP_WZ='linear ramp'/ +&MISC FREEZE_VELOCITY=T, STRATIFICATION=F, NOISE=F, RAMP_WZ='linear ramp', NEAR_WALL_PARTICLE_INTERPOLATION=T/ &SURF ID='bot', NO_SLIP=T, COLOR='BLUE'/ &SURF ID='top', NO_SLIP=T, COLOR='RED', VEL=-10/