Repulsive force of thermal radiation on persons in EVAC

Alex Wang

unread,

Mar 15, 2022, 3:34:23 AM3/15/22

to FDS and Smokeview Discussions

Hello sirs,

I want to add the repulsive force caused by thermal radiation in the EVAC code, to simulate human avoidance of the fire source. However I don't know how to use the vector pointing from fire source to the person. Need I add SUBROUTINE "Find_fire" just like "Find_walls" in the code? Could you please give me some advice?

Best wishes~

TimoK

unread,

Mar 16, 2022, 9:39:27 AM3/16/22

to FDS and Smokeview Discussions

The "FED-file" contains:

WRITE (LU_EVACFED) REAL(FED_CO_CO2_O2,FB),REAL(SOOT_DENS,FB),REAL(TMP_G,FB),REAL(RADFLUX,FB)

SUBROUTINE EVAC_MESH_EXCHANGE(T,T_SAVE,I_MODE, ICYC, EXCHANGE_EVACUATION, MODE)

CALL GET_FIRE_CONDITIONS(NOM,I1,J1,K1,FED_CO_CO2_O2,FED_CO_CO2_O2_REST,FED_CO_CO2_O2_HARD, &
SOOT_DENS,TMP_G,RADFLUX,ZZ_GET)

SUBROUTINE GET_FIRE_CONDITIONS(NOM,I,J,K,fed_indx,fed_indx_rest,fed_indx_hard,soot_dens,gas_temp,rad_flux,ZZ_GET)

! Rad flux, ind=18, kW/m2 (no -sigma*Tamb^4 term)
rad_flux = MAX(MESHES(nom)%UII(I,J,K)/4.0_EB,SIGMA*TMPA4)

So, there is some information of the radiation that you can use in the evac.f90 routines. You could

use this (defined on the (x,y) of the current evacuation mesh, of course). You can calculate the

(local) gradient of rad_flux at the position of the agent, for example. Or just use "line-of_sight" and

see, if rad_flux is large in this path.

TimoK

Alex Wang

unread,

Mar 16, 2022, 8:48:51 PM3/16/22

to FDS and Smokeview Discussions

Tkank you for your reply, Timok. I will try it.

Alex Wang

unread,

Mar 25, 2022, 3:24:09 AM3/25/22

to FDS and Smokeview Discussions

Dear Timok,

I got a question in adding fire force in evac.f90. I tried to use the local gradient of rad_flux as you advised. I used the average values in each directions(+x, -x, +y, -y) for modelling the direction of the force from the fire. The codes I added as follows:

HR%F_fire = HEAT_FAC * HUMAN_GRID(II,JJ)%RADFLUX * (((HUMAN_GRID(II,JJ)%RADFLUX - HUMAN_GRID(II-1,JJ)%RADFLUX + HUMAN_GRID(II,JJ)%RADFLUX - HUMAN_GRID(II,JJ-1)%RADFLUX + HUMAN_GRID(II,JJ)%RADFLUX - HUMAN_GRID(II+1,JJ)%RADFLUX + HUMAN_GRID(II,JJ)%RADFLUX - HUMAN_GRID(II,JJ+1)%RADFLUX) / 4.0_EB) / SQRT (((( HUMAN_GRID(II,JJ)%RADFLUX - HUMAN_GRID(II-1,JJ)%RADFLUX + HUMAN_GRID(II,JJ)%RADFLUX - HUMAN_GRID(II,JJ-1)%RADFLUX + HUMAN_GRID(II,JJ)%RADFLUX - HUMAN_GRID(II+1,JJ)%RADFLUX + HUMAN_GRID(II,JJ)%RADFLUX - HUMAN_GRID(II,JJ+1)%RADFLUX) / 4.0_EB)**2)))

It can be compiled well, however the simulation result seems not so fine? Where am I wrong? Could you give me some advice?

Best wishes for you!

在2022年3月16日星期三 UTC+8 21:39:27<TimoK> 写道：

Alex Wang

unread,

Mar 26, 2022, 10:22:08 PM3/26/22

to FDS and Smokeview Discussions

Also, in the simulation, some evacuees stop moving when they encounter the fire area, it seems that the forces on the evacuees has reached an equilibrium state. Am I supposed to modify the random force?

在2022年3月16日星期三 UTC+8 21:39:27<TimoK> 写道：

TimoK

unread,

Apr 4, 2022, 4:27:29 AM4/4/22

to FDS and Smokeview Discussions

Well, your HR%F_fire should be a vector (in a xy-plane), so it should have two components: HR%F_fire_x and HR%F_fire_y. And then you add the x and y components of the fire force to:

! Save the forces for the loop evac_move_loop (next time step)
HR%F_X = P2P_U
HR%F_Y = P2P_V

And you have lines like:

FC_X = FC_X - FC_DAMPING*(U_TMP(III)-U_TMP(JJJ))*(X_TMP(III)-X_TMP(JJJ))/TIM_DIST
FC_Y = FC_Y - FC_DAMPING*(V_TMP(III)-V_TMP(JJJ))*(Y_TMP(III)-Y_TMP(JJJ))/TIM_DIST
CONTACT_F = CONTACT_F + SQRT(FC_X**2 + FC_Y**2)
CONTACT_FX = CONTACT_FX + FC_X
CONTACT_FY = CONTACT_FY + FC_Y
P2P_U = P2P_U + FC_X
P2P_V = P2P_V + FC_Y

P2P stands for "Person to person force" and it has both x and y components. Well, "person to person" includes also walls:

CALL WALL_CONTACTFORCES(X_TMP(1), Y_TMP(1), R_TMP(1), U_TMP(1), V_TMP(1), D_XY, &
P2P_U, P2P_V, P2P_TORQUE, CONTACT_F, D_WALLS, FOUNDWALL_XY, CONTACT_FX, CONTACT_FY)
CALL WALL_CONTACTFORCES(X_TMP(2), Y_TMP(2), R_TMP(2), U_TMP(2), V_TMP(2), D_XY, &
P2P_U, P2P_V, P2P_TORQUE, CONTACT_F, D_WALLS, FOUNDWALL_XY, CONTACT_FX, CONTACT_FY)
CALL WALL_CONTACTFORCES(X_TMP(3), Y_TMP(3), R_TMP(3), U_TMP(3), V_TMP(3), D_XY, &
P2P_U, P2P_V, P2P_TORQUE, CONTACT_F, D_WALLS, FOUNDWALL_XY, CONTACT_FX, CONTACT_FY)

Above one calculates three times wall forces for an agent, why? The reason is that the agent body is represented by three circles. For each circle one has to see if it touches a wall. In your case, it is enough to use the centre point of the agent, i.e., the cell II,JJ where it is and see, what is the gradient of the red.intensity there. And you need not to calculate any torque due to the force. You do not want to turn the agent body due to the fire force, you just want to change the movement as due to a repulsive force. The agent nose is pointing where the desired velocity V0 is pointing (the angular equation has "motive torgue" such that it makes the agent to turn its nose towards the desired walking velocity direction).

The other option would be to change the desired velocity vector V0 of the agent due to the "fire force". In this case, you would not add anything to P2P_U not HR%F_X, but you would turn the V0 vector to some other direction. See;

! ========================================================
! Collision avoidance, counterflow, etc.
! ========================================================

DO III = 1, N_SECTORS+1
IF (SUM_SUUNTA(III) > SUM_SUUNTA_MAX) THEN
UBAR = U_THETA(III)
VBAR = V_THETA(III)
SUM_SUUNTA_MAX = SUM_SUUNTA(III)
I_SUUNTA_MAX = III
END IF
END DO

So, the agents check three direction in head of it and sees, what is the best. And then it changes the V0 (x and y components UBAR and VBAR, respectively). One of the V0 directions in the "counter flow" model is the original V0 direction. So, you could modify the UBAR,VBAR vector so that the agent does not go towards the negative gradient of the rad.intensity (gradient points away from the fire, gradient "goes downhill"). This might be easier than to modify the force.

So, you should use the gradient of the rad.intensity. A gradient is a vector.

TimoK

Alex Wang

unread,

Apr 4, 2022, 10:19:14 PM4/4/22

to FDS and Smokeview Discussions

Dear Timok,

Thank you for your kindly response. I will try your advice.
By the way, I want to calculate the FED value caused by heat. Refered to SFPE user guide, I added these codes in evac.f90.

FED_HEAT = 0_EB

FED_HEAT = (HUMAN_GRID(II,JJ)%RADFLUX**1.33_EB/10.0_EB)+(HUMAN_GRID(II,JJ)%TMP_G**3.4_EB/5.0E7_EB) (! equations from SFPE user guide to caluculate FED value)

Then the total FED value was modified as

HR%INTDOSE =DTSP*HUMAN_GRID(II,JJ)%FED_CO_CO2_O2 + DTSP*FED_HEAT/60._EB + HR%INTDOSE

However, the simulated FED values were much greater than 1, I can't figure out where am I wrong. Could you give me some guidance? Your advice will be much appreciated.

TimoK

unread,

Apr 6, 2022, 3:53:07 AM4/6/22

to FDS and Smokeview Discussions

Well, you should add a separate FED counter for the agents, so HR%INTDOSE_TEMP. So, keep these separate.

Well, I should check the SFPE handbook, it might be that the toxic gas FED and temperature FED are combined there.

If it is so, then just one INTDOSE is enough.

The time step DTSP is in the units of seconds. Is the 60.0 here, because the SFPE equations use minutes?

DTSP*FED_HEAT/60._EB

The other units: RADFLUX in fds source is, probably, W/m2. The SFPE might use kW/m2. Or both might be kW/m2.

RADFLUX

In the dump.f90 you find lines like:

CASE(10) ! TOTAL HEAT FLUX
SOLID_PHASE_OUTPUT = (ONE_D%Q_RAD_IN-ONE_D%Q_RAD_OUT+ONE_D%Q_CON_F)*0.001_EB

CASE(18) ! INTEGRATED INTENSITY
GAS_PHASE_OUTPUT_RES = UII(II,JJ,KK)*0.001_EB

CASE(19) ! RADIATION LOSS
GAS_PHASE_OUTPUT_RES = QR(II,JJ,KK)*0.001_EB
So, it might be that fds source uses SI units (W,s,K,m,J,...). And in the output routines the values are changed

to Celsius and kW (and kJ). By checking the source code you should easily find the units of the radiative intensity.

And the TMP_G is Kelvin in the source code, the SFPE might use Celsius (or even Farenheit).

HUMAN_GRID(II,JJ)%TMP_G**3.4_EB/5.0E7_EB

See the human_type "fortran type" in type.f90:

TYPE HUMAN_TYPE

clip-clip-clip...

REAL(EB) :: SumForces=0._EB, IntDose=0._EB, DoseCrit1=0._EB, DoseCrit2=0._EB, SumForces2=0._EB

So, the varibales like HR%X, HR%INTDOSE, etc. are defined in the HUMAN_TYPE type. In evac.f90 the HUMAN_TYPE is

used to define new variables etc. So, you just copy the INTDOSE stuff to INTDOSE_TEMP, so the programming should

be quite easy. You just change the arithmetic formulas.

Alex Wang

unread,

Apr 7, 2022, 9:56:40 AM4/7/22

to FDS and Smokeview Discussions

Thank you Timok, your reply is really precious for me, and now it finally works.

Yes, the SFPE equations use minutes, so I apply the formula DTSP*FED_HEAT/60._EB.

But something goes wrong with the output data, while in the _evmc.csv file, "NaN" appears in the columns of FED_max and FED_max_alive. It must be because I changed the calculation formula of FED. So where am I supposed to modify, in order to get the FED_max (calculated by both gas and heat) I want in the .csv file. In addition, could I output an extra column named "FED_heat_max" just like FED_MAX does? Then what else should I do?

Best wishes~

Alex Wang

unread,

Apr 22, 2022, 10:12:27 PM4/22/22

to FDS and Smokeview Discussions

Dear Timok,

I got a question about the TEM_G unit you mentioned earlier. As you said the TMP_G is Kelvin in the source code, I changed the FED_HEAT formula as: (HUMAN_GRID(II,JJ)%TMP_G - 273.15_EB)**3.4_EB/5.0E7_EB.

However I got "NaN" in the csv file. It seems that the value (HUMAN_GRID(II,JJ)%TMP_G - 273.15_EB) less than zero. When I removed - 273.15_EB, the result in the csv file seems to be OK.

Besides, some codes like "TMP_G=TMPA", isn't TMPA the temperature value (Celsius) input by the user? If so, the unit of TEM_G should be Celsius, right? I don't know if my understanding is wrong, and looking forward to your reply.

Best wishes~

Alex

在2022年4月6日星期三 UTC+8 15:53:07<TimoK> 写道：

Kevin McGrattan

unread,

Apr 23, 2022, 8:33:54 AM4/23/22

to fds...@googlegroups.com

All temperatures are input to, and output from, FDS in degrees C. However, within the code these temperatures are all converted to K. So if you are modifying source code, assume that the temperature, like TMP_G, is in K.

TimoK

unread,

May 13, 2022, 6:48:18 AM5/13/22

to FDS and Smokeview Discussions

I have been a little bit busy lately, but now have some time.

Yes, you could add/modify evac.f90 and there:

SUBROUTINE DUMP_EVAC_CSV(Tin)
!
! Dump agent data to CHID_evac.csv
! This subroutine is called from the main program.

And earlier in evac.f90 the column headers are written in SUBROUTINE INITIALIZE_EVAC_DUMPS:

WRITE (tcform,'(a,i4.4,a)') "(",n_cols+3+(j_density-j_ntargets),"(a,','),a)"
WRITE (LU_EVACCSV,tcform) 's','AgentsInside', &
('AgentsInsideMesh', i=1,n_egrids), &
('AgentsInsideCorr', i=1,n_corrs), &
('ExitCounter', i=1,N_EXITS), &
('DoorCounter', i=1,N_DOORS), &
('TargetExitCounter', i=1,N_EXITS-n_co_exits), &
('TargetDoorCounter', i=1,N_DOORS), &
('DensityCounter', i=1,j_density-j_ntargets), &
'Agents','FED_Index','FED_Index'
WRITE (LU_EVACCSV,tcform) 'EVAC_Time','AllAgents', &
(TRIM(EVAC_Node_List(i)%GRID_NAME), i=1,n_egrids), &
(TRIM(EVAC_CORRS(i)%ID), i=1,n_corrs), &
(TRIM(EVAC_EXITS(i)%ID), i=1,N_EXITS), &
(TRIM(EVAC_DOORS(i)%ID), i=1,N_DOORS), &
(TRIM(CTEMP(i)), i=1,N_EXITS-n_co_exits), &
(TRIM(EVAC_DOORS(i)%ID), i=1,N_DOORS), &
(TRIM(CTEMP(i)), i=j_ntargets+1,j_density), &
'Number_of_Deads','FED_max','FED_max_alive'

Check that your are defined your fed_heat_max similarly as fed_max, fed_max_alive. They should be defined in the

same place as fed_max_alive and fed_max. These are in the common declarations of the whole module EVAC, i.e,

they are common to all subprograms in the evac.f90 file:

REAL(EB) :: fed_max_alive, fed_max

And search for fed_max_alive and fed_max to see, where they are initialized to zero etc., so put your

fed_heat_max similarly.

TimoK

Alex Wang

unread,

May 13, 2022, 9:02:08 AM5/13/22

to FDS and Smokeview Discussions

Hi Timo,

Much glad to receive your comments. However, I got an another question puzzled me recently. And I need your guidance.

I am not sure what RADFLUX means in evac.f90...I just want to calculate the FED value caused by fire radiation using RADFLUX. I want to calculate the FED_HEAT value when the heat radiation value is greater than 1, so I added the following command in evac.f90:

IF (HUMAN_GRID(II,JJ)%RADFLUX * 0.001_EB > 1.0_EB) THEN ! RADFLUX in fds source is W/m2

FED_HEAT =((HUMAN_GRID(II,JJ)%RADFLUX * 0.001_EB)**1.33_EB)/(600.0_EB)

ENDIF

HR%INTDOSE_HEAT=DTSP*FED_HEAT + HR%INTDOSE_HEAT

However, the calculation results can't consider the change of thermal radiation in the fire field. I even got zero as the agent pass through the fire source.(HRR = 2000KW). I'm wondering if I have missuse "RADFLUX" (is it represents the "INTEGRATED INTENSITY" in the current evacuation cell?), or am I missing some setting issues?