Correlation between the wall roughness and friction factor

[florinmtl]

Hello all,

I am presently in the situation of validating the design of the ventilation system for a road tunnel. The use of FDS has been challenged by one of our international consultant concerning the computing of the linear pressure losses, i.e. the correlation between the wall roughness and Darcy’s friction factor “f”?
So, in order to clarify the correlation between the roughness of the walls and Darcy’s friction factor value, I performed several simulations without fire, in the empty tunnel, and could not reach any conclusion (see 180720PressLossNoSlip5_devcLambda.fds).

Model’s input and simulation parameters:

Tunnels geometry and characteristics:

1. 1390 m (L) x12.8 m (W) x5 (H) m tunnel;
2. hydraulic diameter Dh=7191 mm
3. forced ventilation: 2 fans x 250 m3/s exhaust, situated at 950 m from the entry portal;
4. no vehicles;
5. no fire;
6. 24 meshes;
7. cell dimensions: 1 m * 0.512 m * 0.5 m;
8. Walls’ material: Concrete (characteristics from FDS library);
9. Wall roughness: 5 mm.
   Here it is the description of the surface:

&SURF ID='WallConc',
BACKING='INSULATED',
MATL_ID(1,1)='CONCRETE',
MATL_MASS_FRACTION(1,1)=1.0,
THICKNESS(1)=0.25,
NO_SLIP=.TRUE.,
ROUGHNESS=5.0E-3/

(The use of NO_SLIP=.TRUE. seems to be an error; if so, my understanding of the user guide was erroneous)

10. Standard ambient conditions:
    a. temperature 20 deg C;
    b. density rho=1.206 kg/m3;
    c. atmospheric pressure p= 1325E5 Pa,
    d. kinematic viscosity of the air nu = 1.506*10^-5 m2/s

11. Pressure sensors all over the tunnel, at 0.1 m from the wall and 2.5 m high spaced at 100 m, with the following generic description:

&DEVC ID='PresNi', QUANTITY='PRESSURE', XYZ=50.0,0.1,2.5, ORIENTATION=0.0,1.0,0.0

FDS results:

1. upstream the exhaust fans positive velocity of 2.6 m/s;
2. downstream negative velocity negative 4.92 m/s (air is drawn from the exterior).

The mean values of the pressure calculated in a stabilized regime, have been used to compute the “f” factor:

f=2deltaPressure_i,i+1Dh /(2rhodeltaLength_i, i+1).

The values are (see 180720PressLossNoSlip5_devcLambda.xlsx):

Upstream (velocity 2.6 m/s):
“f” calculated with FDS results : 0.09400 ; 0.10904 ; 0.11788 ; 0.118980 ; 0.118578 ; 0.118098 ; 0.118413 ; 0.117669 for Re = 1251025

Downstream (velocity of 4.92 m/s):
“f” calculated with FDS results : 0.13456 ; 0.11691for Re : 2349254

These values have been compared with those calculated with the Colebrook equation:

f = {-2log10[e/(D3.7)+2.51/(Re*f^1/2)]}^-2

Computed values:

1. upstream: Re=1 251 025; f=0.01833564

2. downstream: Re=2 349 254; f=0.01818287

The differences are huge and the values of « f » computed from FDS results are way too big.

On the other hand, if there is no fire, the values of the static pressure measured by sensors must always be in my opinion negative because I assume that the devices are measuring gauge pressure, so inferior to the atmospheric pressure (In the simulation whose results are presented above the values were negative). This happens in some simulations but not in all.

I studied the examples given for the Moody validation, and I ran some of them.
For example in “z0=p0001_dpdx=-1_N8-press.fds” I modified the length of the mesh to 20m and I added some pressure devices.
In this case, I also find the pressure increases along the tunnel (so “f” results negative).

Thus being said, will you please help me to get in the good direction?
The examples given in the validation guide are all using WIND parameters, but pressure devices do not give the same pressure gradient (Pa/m).
Should I use “Statistic” cross-sections instead of devices?
On the other hand, if the fire is present, the friction factor varies downstream because of the Re and viscosity, so the constant pressure gradient hypothesis does not function anymore.
Also, the dimensions of a real tunnel are not comparable with those used in the simulations, the air/gases are flowing in both directions, etc.

So, for a road tunnel with mechanical ventilation (exhaust fans, jet fans, Saccardos, etc) I would like you to indicate the best practice to correlate the roughness and the linear pressure losses.
Also, if I have to compare the values of “f” calculated from Colebrooke and Darcy, which is the best way to measure the mean values of the pressure in different cross sections?

I thank in advance for your precious help.
Best regards,

Florin

180720PressLossNoSlip5.zip

[mcgratta]

Your case is too complicated to sort out these fundamental issues. Do this: create a very, very simple input file with a horizontal tunnel with approximately the same dimensions and shape as your real tunnel. Blow air through the tunnel and let's see if the model achieves a pressure drop that is comparable to empirical correlations.

[rmcdermo]

Also, you should not have NO_SLIP and ROUGHNESS specified on the same SURF. These are mutually exclusive parameters. We should throw an error for this input. I will fix this.

…

On Thu, Aug 2, 2018 at 11:12 AM, Kevin McGrattan ***@***.***> wrote: Your case is too complicated to sort out these fundamental issues. Do this: create a very, very simple input file with a horizontal tunnel with approximately the same dimensions and shape as your real tunnel. Blow air through the tunnel and let's see if the model achieves a pressure drop that is comparable to empirical correlations. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#6738 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AENrwUe59Asl47v7sFIb8tkMPMO9UH-kks5uMxbpgaJpZM4VsdNy> .

[rmcdermo]

Also, when you use WIND there is a mean force being applied to the right hand side of the FDS momentum calculation. This force acts like a mean pressure gradient. The DEVCs for pressure will not pick this up. What the DEVC reports is the local hydrodynamic pressure, which sits on top of the mean pressure gradient from the nudging scheme used by WIND. We should probably create a DEVC that integrates the mean nudging force in each component direction so that this could be included in the force balance. On Thu, Aug 2, 2018 at 11:18 AM, Randy McDermott <randy.mcdermott@gmail.com> wrote:

…

Also, you should not have NO_SLIP and ROUGHNESS specified on the same SURF. These are mutually exclusive parameters. We should throw an error for this input. I will fix this. On Thu, Aug 2, 2018 at 11:12 AM, Kevin McGrattan ***@***.*** > wrote: > Your case is too complicated to sort out these fundamental issues. Do > this: create a very, very simple input file with a horizontal tunnel with > approximately the same dimensions and shape as your real tunnel. Blow air > through the tunnel and let's see if the model achieves a pressure drop that > is comparable to empirical correlations. > > — > You are receiving this because you are subscribed to this thread. > Reply to this email directly, view it on GitHub > <#6738 (comment)>, > or mute the thread > <https://github.com/notifications/unsubscribe-auth/AENrwUe59Asl47v7sFIb8tkMPMO9UH-kks5uMxbpgaJpZM4VsdNy> > . >

[florinmtl]

Hello Kevin,

I shall do this.
Please look to the "z0=p0001_dpdx=-1_N8-press.fds" issue I described. The pressure is increasing in the direction of the flow. Could you please explain?

[rmcdermo]

Where is this input file? I do not see it in the zip.

[rmcdermo]

Regardless, I've already explained why the DEVC will not work for giving pressure loss when you use a mean force.

[florinmtl]

180720PressLossNoSlip5.zip
Sorry I forget to update it.

[rmcdermo]

Like I just said, DEVC for pressure will not work for this. Study the Matlab scripts that are used to post-process the Moody chart cases. You have to do an integral momentum balance based on the specified FORCE_VECTOR.

https://github.com/firemodels/fds/blob/master/Utilities/Matlab/scripts/fds_moody_chart.m

https://github.com/firemodels/fds/blob/master/Utilities/Matlab/scripts/friction_factor_calc.m

[florinmtl]

Hello Randy,

I do not want to use WIND. The previously studied wind influence on the tunnel is negligible.
We are using Saccardos and powerful extraction shafts and we're simulating real traffic/fire conditions: stalled vehicles and a very high power fire.
In these conditions the friction losses are significant. That is why the values of the flow rate and the aerodynamic performances of the equipment must be accurate.

On the other hand, the measure of the dynamic pressure is direction dependent and strictly positive.
If the DEVC category measures the dynamic pressure (rho*v2/2), it should be oriented opposite to u-vel. In these conditions, the measured value should be positive. I put pressure devices back to back (x=-1, x=1) at the same point or very close, on the same abscise and the results are not very different.

[mcgratta]

Then do not use WIND. Model your tunnel using specified values of VEL at boundaries.

[florinmtl]

Hello again Randy,

The “.m” files are clear and I suppose that these are the results used in the Technical Note 1640. However, my problem persists because the dpdx [Pa/m] was imposed as input, but in my case, I have to measure it in order to compute “f”.
I also remarked that the dynamic viscosity order of value is 10e-5 Pa/s or kg/m/s (“mu = 1.84e-5” in “.m” files) and the one calculated by fds in the validation examples (z=0….fds, moody*) is 10e-2 kg/m/s. While density unit is kg/m3, the dynamic and the kinematic viscosities should both be 10e-5.
In conclusion, please give me some advice on what I do should:

1. Use SURF and ROUGHNESS, FREE_SLIP=.TRUE. and measure the mean pressure values using STATISTIC planes?
2. Is it any precaution to take regarding the meshing near the wall, ratios, etc?
   Thanks

[rmcdermo]

I created a simple test case with a long tunnel, specifying the inlet velocity and using the sand-grain roughness you specified in your example. I added a line of DEVC for pressure. I compare the friction factors from FDS and the Colebrook equation and they are 0.0090 and 0.0093, respectively. The input file, output files, and Matlab scripts to process the case are attached in the zip.

Archive.zip

Things to be careful about:

1. If you are using multiple meshes, you need to make sure your mesh to mesh velocity errors are near zero, else you will get strong oscillations in the tunnel volume flow. This transient in mean velocity is indistinguishable from a pressure gradient as far as the momentum equation is concerned.
2. Probably not as big an issue, but watch out for entrance length effects as well. For high Reynolds number you may see L/D entrance of 50 or more.

[rmcdermo]

Do not set FREE_SLIP or NO_SLIP. The effect of these parameters is explained in the guides.

Your near wall grid resolution should put you at about y+ of 100 or less. You can monitor this with BNDF of 'VISCOUS WALL UNITS'.

[rmcdermo]

Actually, the upper range on y+ can be much higher. In my example, I get about 5000. But that is because Re is 6.6e6, which is very high.

[rmcdermo]

Also, the VISCOSITY that FDS outputs is the turbulent + molecular. This is not the value you should use the compute the Re for the friction factor. To get Re you only use the molecular viscosity. Either look in the .out file or use an output of 'MOLECULAR VISCOSITY'.

[florinmtl]

Hello Randy,

1. I could not find any documentation about the pressure DEVC line and don't know how it works. As far as I can see there are 160 (or 162) press devices on the mid-axis of the tunnel but for each of it, the pressure (s) should be constant in time. Will you be so kind and give me some references about this item?
2. Could the pressure sensors be replaced by statistic cross-sections?
3. I don't understand why the attached diagram is time-dependent when the simulation time is only 160s. I suppose that the plot should be against x.
4. Please clarify what pressure is measured. In the .m file there is a "gauge pressure" label, which by default must be static pressure. In these conditions, I cannot understand the values of the pressure. If it is static pressure, we have u-vel constant=10 m/s, then pdynamic = 59.8 Pa and the pressure values should be P(x) = [101325-59.8-frhoVel^2/(2HydraulicDia)]x [Pa]. So, what do these values represent?
5. Finally, I found 0.009388 from Colebrook. Calculating with the values within "tunnel_pressure_drop_line.csv" (values that I do not understand) the f is 0.00898 so 0.009.
6. Just for curiosity: why the density of the air is 1.196 kg/m3 when at 20 degC rho = 1.204 and at 40 % humidity rho = 1.2024?

Thanks for your patience and help and I wish you a great weekend.
Florin

[rmcdermo]

1. FDS User Guide, 20.2.2 Linear Array of Point Devices
2. Yes, and should be. I've continued to work on this case, and interestingly the f does not work out as nicely for higher resolution. I suspect that spatial and time averaging are problems here. I'm working on it.
3. Yes, should be x.
4. GAS_PHASE_OUTPUT_RES = PBAR(KK,PRESSURE_ZONE(II,JJ,KK)) + RHO(II,JJ,KK)*(H(II,JJ,KK)-KRES(II,JJ,KK)) - P_0(KK) PBAR is thermodynamic pressure in equation of state. It can rise in a sealed compartment. P_0 is the atmospheric pressure at discrete height Z(KK). In this case they cancel. RHO*(H - KRES) is the fluctuating hydrodynamics pressure. See section 2.4 Low Mach Number Approximation of the FDS Tech Guide.
5. Why do you not understand them? See 4.
6. See .out file for composition of AIR.

[drjfloyd]

6. Air at 20 C and 40 % relative humidity has an MW of 28.76 g/mol (as Randy indicated see the .out file). The default ambient pressure is 101325 Pa.

0.02876 * 101325 / (8.314472 * 293.15) = 1.196

[florinmtl]

Hello, gentlemen,

And thanks for the answers.
Will you please send me the .out file?
I'm looking forward to hearing from you about the new results concerning the correlation that has been the starting point of these exchanges.
I think that the release of an abstract gathering the best practices (meshing, sensors, etc.), as well as special instructions to add, “by hand”, in the fds files for those who, like me, are using software interfaces to generate fds files, will be extremely useful for the design engineers.
Once again many thanks for your help.
Best regards,

Florin

[rmcdermo]

The .out file is created by FDS for every run. You have it in your working directory.

[florinmtl]

Sorry,

I forgot I have run the simulation myself.
Thanks

[rmcdermo]

I am continuing to work on this and today noticed something interesting that was not obvious to me at first. By pushing flow in from the inlet we get different behavior for the PRESSURE DEVC output than if we suck the flow from the outlet side.

The reason for this is that the two boundary conditions for H are different. By specifying the velocity on the outlet side, we pin down the velocity, but allow the pressure to float to satisfy momentum, and the value of H gets pinned to a constant at the inlet. If you do not do this, the magnitude of H fluctuates in time and this leads to erroneous output for the PRESSURE DEVC.

We need to work this on our side, since the natural way to specify the flow is via the inlet.

I've now verified that by sucking the flow from the outlet I can get the same correct friction factor for both coarse and finer resolution, but I need to get the plots put together. Hopefully have something tomorrow.

[florinmtl]

Hello,

I followed Kevin's recommendations and I built one simple model, with one mesh and with 24 meshes (for 24 cores), blowing air through one supply vent. Alongside with the pressure_line _devc, I "installed" a number of statistics pressure planes.
The results are bizarre :

1. the pressure highly fluctuates on the planar surfaces in steady-state conditions. (see figure).
   
2. The values of the friction factor given by fds are 0.05629242 and 0.075082 and that computed by Colebrook eq. 0.02144.
   Friction.zip
   Just a small remark concerning the ventilation strategies: ideally is to exhaust the smoke at the source. However, for old or immersed tunnels this is pretty difficult and the exhaust shafts are a common solution. For the longitudinal ventilation is more popular to use jet-fans (only momentum) or Saccardos, whose nozzles are frequently located on the ceiling or on the lateral walls and not at the portals. So, the portals are in general "OPEN" vents.
   Best regards,

[rmcdermo]

@florinmtl I am still working on this, preparing a systematic verification case. I have learned a few things:

1. For multimesh cases, you must have very tight velocity tolerances. For now, I would recommend just using GLMAT solver.

2. At least use CFL_VELOCITY_NORM=2. The default value of 0 leads to a few instants where the pressures really can go unstable.

3. With regard to what strategy you use for pushing flow in your application, for now you must think of this exercise as a verification that FDS does the wall stress correctly. Once you start changing the pressure boundary conditions away from sucking from the outlet, all bets are off with the pressure devices. The magnitude of the pressure in the tunnel will adjust to satisfy the bcs. There is no pinned Dirichlet value for the pressure. Once you verify the pressure drop in a simple case, you have to use that as evidence that FDS is computing the stress correctly. The only term that shows up in the momentum equation at an instant in time is the pressure gradient---the pressure can be of arbitrary magnitude as long as the gradient is right.

[rmcdermo]

OK, I have to revise my conclusions somewhat. The good news is that I think we are fine to push the flow through the tunnel with DEVC for pressure. The problem I was seeing before was mainly due to erratic behavior caused by the CFL velocity norm resulting in a time step that was not truly stable. So, the keys are:

1. Use GLMAT solver for multi-mesh.
2. Use CFL_VELOCITY_NORM=2 (or better, e.g., 1 would also work)

With these settings I get f_10 = 0.0088 and f_20=0.0085 and f_Colebrook = 0.0093, where 10 represents 10 cells across the height the tunnel, etc.

[florinmtl]

Hello,

1. Does "10 cells across the height " mean that you computed "f" at different heights?
2. I ran another simulation (started yesterday, still running) with the same tunnel but with a smaller meshing (refine factor: 2), meaning that the cell is 0,50,2560,25, and the results are totally incomprehensible, for me at least. Tunnel straight, rectangular cross-section, one portal in supply at 4 m/s, no obstacles, all wall roughness 10 mm. Al other conditions standard. A pressure line devices on the centerline, pressure devices on the same line oriented in various directions (+x, -x, -z, +y).
   Please look at the results. I don't have any explanation. Please help me to understand.
   
   
   
   
   
   
   

180808PressLoss10 .fds.txt
180808PressLoss10.zip

[rmcdermo]

I said you need to use the GLMAT solver and set CFL_VELOCITY_NORM=1. Here is a sample input file.

tunnel_pressure_drop_10.fds.txt

And, no, there is only one "f" for the tunnel.

[florinmtl]

1. Do I have to understand though that you performed 2 simulations with &MESH IJK=10,10,10 and &MESH IJK=10,10,20?
2. L/D entrance = 60 means that the sensors are giving reliable results only if installed starting with L=60D (600 m in your case) to avoid effects due to high turbulence?
3. Why measure RHO, MU, U_TAU, Y_PLUS?
   Will this give ac future information for a better mesh pattern?

[rmcdermo]

1. Yes, I performed simulations at two grid resolutions to see what the effect was. It was not much.
2. At first I was worried about the entrance effects. But that was due to other issues regarding the CFL.
3. I am just double checking that I get what I expect.

I've finished writing up the section for the FDS User Guide. It will be available after tonights Firebot run. Here is section of the guide.

FDS_User_Guide_tunnel_pressure_drop.pdf

[florinmtl]

Thank you very much.
I hope that this solves and closes the problem. I shall try the simulation with a 10 mm roughness and let you know if the results are still correlated.
All the best!

[florinmtl]

Hello gentlemen,

I performed 2 simulations with 10 mm roughness (instead of 0.1 mm used in your examples), 1390(L)x12.8(W)x5(H) m tunnel, only one mesh of 1x0.512x0.5m and 0,5x0.512x0.5m. I used all the setup given in FDS_User_Guide_tunnel_pressure_drop.pdf and in the previous exchanges.
Conclusions:

1. The friction factor calculated from "pressure DEVC line" using the regression line given by Excel and pressure statistic planes are similar. However, the pressure plane values are different in the first and last 200 m of the tunnel (i.e. 0.18996, 0.17923 0 to 200m from entry, and 0.17215 at 200 m to exit ). Mean value: 0.01833.
2. The value calculated from dev line : 0.018338.
3. The value calculated with Colebrook: 0,021229.
4. Error: -13%. Still too big.
5. The mesh dimensions did not seem to influence the values. Is it because only the x dimension was reduced?
   I'm looking forward to hearing from you.

[rmcdermo]

Post the input file and I'll look at it when I get time. But if you think that CFD is going to give you a real-world tunnel pressure inside of 13% error, I think you are being unrealistic.

[florinmtl]

180810PressLoss10-NIST1MESH05.fds.txt
180809PressLoss10-NIST1MESH.fds.txt
Thanks. Even if the tolerance is 20 % if I knew it from the beginning is a step forward.
In the examples given in the “FDS_User_Guide_tunnel_pressure_drop” the error is under 10 % and reduces with the mesh dimension for roughness of 0.1 mm but increases for 100mm. I just wonder if it is not because of these extreme values of the roughness (for example, for concrete tunnel structures, both of them are out of normal range). Thanks once again. Best regards.

[rmcdermo]

I don't know what you mean by "normal range" for roughness in a tunnel. The correlations we use for friction factor were developed for pipes and pipelines. I've never been in a tunnel that I thought could be described with exactly the same correlation as a pipe. You have to use some engineering judgement.

[rmcdermo]

The results are extremely sensitive to how you choose to define the length scale in the Reynolds number. The hydraulic diameter for the case you sent me is about 7.2 m. However, the D of an equivalent area circle is about 9.0 m. An almost perfect match with the friction factor is achieved using 8.2 m. I leave the rest to you. I cannot keep spending time on this problem, as I do not see an error in FDS.

tunnel_pressure_drop_e_DH7p2.pdf
tunnel_pressure_drop_e_DH9p0.pdf
tunnel_pressure_drop_e_DH8p2.pdf

[rmcdermo]

I tried one more thing, I set the near wall eddy viscosity model to WALE.

&MISC CFL_VELOCITY_NORM=1, NEAR_WALL_TURBULENCE_MODEL='WALE'/

Here is the result using DH=7.2, an error 2.3% relative to Colebrook.

I also tracked down this paper, which is worth reading, in that it shows some of the scatter one can expect in the friction factor data.

WALE model will be under consideration for the default for next version. But we first have to run through all validation tests, which will take some time.