Culverts issue

[mikilterribile]

Hi everyone, I am trying to modelling some stuff with the Anuga's culvert operator. I use the box operator with a rectangular section 10m * 1m. The problem is that the flow rate is lower than what expected either if the culvert is gravity led or under pressure.

The amount of water expected, infact, is something like 100 m^3/2 (calculated with the uniform flow formula) while the value calculated by the culvert operator is about 50 m^3/s. I have tryed also to use e pipe culvert operator, instead the box operator, but the flow rate decrease to 40 m^/s.

NB: I use the operator with two points (start and end) and not with two lines.

[stoiver]

@mikilterribile our unittests do try to test against expected flows. But maybe you have found an example which we haven't tested! Actually just looking at the unittest code (test_boyd_box_operator) it looks like we have only tested using exchange lines! I will try to find some time to adapt the current tests to cover using end points as well.

Another possible problem is that geometry of your input regions (as generated by the end points) might be intersecting with the bank of your river.

And finally another problem is that as the operator extracts water from the inlet regions, this can reduce the depth of water in that region, and so affect the flow calculation. We provide the option of setting an enquiry_point (which should be a little upstream of the culvert inlet) to use when calculating the discharge through the culvert.

[stoiver]

@mikilterribile I had a closer look at our unittests and we already have tests using end_points, so hopefully that means at least the underlying calculation of the discharge should be ok (or at least following Michael Boyd's formulas).

So maybe the problem is with how your inlets are being embedded into the topography, or where the enquiry points are situated.

[mikilterribile]

Stephen, thank you for your answer. I don't know the Michael Boyd formula... I did a simple uniform flow calculation, with the assumption that the culvert is not going under pressure. So, I calculated the slope of the culvert as the difference from the initial and the final elevation and I have: if = 2.6/148 = 0.0177. Then I know the maximum flux rate : Q = 70 m^3/s.

Finally the size of the culvert is width b = 10 m and height h = 1 m and the manning coefficient I have choose is M = 0.01. The formula for a rectangular channel is Q = b Y 1/M ((bY)/(b + 2Y))^(2/3) * sqrt(iF) where the unknown variable is Y (the water depth).

Solving this implicit equation for Y finds: Y = 0.71. So the water depth inside the culvert in the moment of maximum flux rate should be 71 cm, while the culvert is 1 met ![Uploading image_for_stephen_0.png…]() er height. So I was expecting all the incoming fulx to be drained by the culvert operator without troubles. But this not happend, with the culvert draining less than 50 m^3/s and the water depth level increase of meters upstream of the culvert and conseguently going out of the stream bed (esondation).

[stoiver]

@mikilterribile I see you have a tricky problem. My guess is that there is a positive feedback which is conspiring to reduce the flow in the culvert, ie the flow hits your riverwall and that slows down the flow (and increases the depth) and consequently the flow around the riverwall (out of the stream).

By the way, I assume you have set losses to 0.0?

You should be able to improve the flow by increasing the enquiry_gap so the culvert operator samples the flow upstream where the velocity is higher. The operator should then take the higher velocity into account when calculating the expected flow through the culvert.

An unsatisfying remedy would be to artificially increase the width of the culvert to match the expected discharge.

I'll try to create a simple example which mimics your situation and see if I can track down the problem. I do think it is to do with the way velocity to used in the formulas.

[mikilterribile]

Ok thank you very much. Yes, I have set losses to 0.0. Best regards, Michele Zucchelli.

[mikilterribile]

A thing that maybe can help is that the culvert flow rate calculated by anuga seems to be completely indipendent from the roughness (the manning coefficient). I have tryed to set manning = 0.000001 and manning = 1 with the same results, despite in my code I pass this coefficient: "culvert_operator = anuga.Boyd_box_operator(domain, losses=losses, width=size[0], height=size[1], end_points=[ep0, ep1], use_momentum_jet=True, # False to switch it off use_velocity_head=True, # False to switch it off manning=manning, # culvert manning’s n barrels=n_barrels, blockage=blockage, logging=True, # False to switch it off label=name, verbose=False)"

[stoiver]

@mikilterribile Could you test your example using lines. The problem may be due to the boyd box operator spreading out water to the triangles making up the exchange region which are determined by the triangles that intersect the line that is constructed when when using endpoints. The length of the line will be 10m centred on the end point. That may be how the water is escaping the channel. When using exchange lines you can make the lines short enough to ensure the triangles constituting the exchange region are all within the channel. Keep the width of the culvert at 10m for the flow calculations.

[mikilterribile]

Hi Stephen, I have tried with the exchange lines, both with 10 meters length and 4 meters length, centred, but nothing has changed. The flow rate remain about to 50 m^/s when it should be much greater. I copy the log file below, maybe it can help...

===============================================
Parallel Structure Operator: culvert_duina_1_P0
===============================================
Structure Type: boyd_box
Description

Culvert   Height: 1.0
Culvert    Width: 10.0
Culvert Blockage: 0.0
No.  of  barrels: 1.0
-------------------------------------
Inlet 0
-------------------------------------
=====================================
Inlet
=====================================
======> inlet triangle indices and centres and elevation at P0
[ 381  960 1004]
[[529.74679372 355.30901605]
 [526.0886286  356.74457056]
 [532.85158162 357.16628503]]
[393.77996392 393.78365624 393.83232999]
Elevation range of 0.0523660779069246Warning: Non-constant inlet elevation can lead to well-balancing problems
Enquiry point:[532.85158162 357.16628503]
Enquiry Index:1004
line
True
-------------------------------------
Inlet 1
-------------------------------------
=====================================
Inlet
=====================================
======> inlet triangle indices and centres and elevation at P0
[879]
[[598.32357788 498.61281284]]
[390.88276001]

Enquiry point:[598.32357788 498.61281284]
Enquiry Index:879
line
True
=====================================

time,discharge_instantaneous,discharge_abs_timemean,velocity_instantaneous,driving_energy_instantaneous,delta_total_energy_instantaneous
72.00000, 0.00000, 0.00000, 0.00000, 0.00000, 2.94957
144.00000, 0.00000, 0.00000, 0.00000, 0.00000, 2.94957
216.00000, 4.97705, 1.43380, 1.71018, 0.44784, 3.16023
288.00000, 11.32609, 8.57376, 2.27405, 0.79185, 3.29345
360.00000, 13.35363, 12.37612, 2.36627, 0.85737, 3.28335
432.00000, 14.65278, 13.96322, 2.45798, 0.92512, 3.32559
504.00000, 15.93456, 15.30282, 2.56130, 1.00453, 3.39316
576.00000, 17.07407, 16.57026, 2.56074, 1.00409, 3.32427
648.00000, 18.38174, 17.80915, 2.66075, 1.08405, 3.39792
720.00000, 19.70719, 19.09726, 2.72153, 1.13415, 3.41623
792.00000, 20.82740, 20.36953, 2.76170, 1.16787, 3.41615
864.00000, 22.09996, 21.55829, 2.84057, 1.23553, 3.47272
936.00000, 23.14616, 22.77393, 2.85650, 1.24943, 3.44152
1008.00000, 24.56336, 23.95898, 2.95208, 1.33836, 3.53288
1080.00000, 25.40500, 24.79515, 2.97768, 1.39642, 3.57124
1152.00000, 26.25414, 25.66168, 3.01486, 1.48428, 3.64491
1224.00000, 27.72216, 26.61309, 3.02453, 1.50784, 3.61780
1296.00000, 28.33364, 27.59160, 3.07576, 1.63769, 3.75157
1368.00000, 29.56300, 28.61727, 3.10264, 1.70931, 3.79065
1440.00000, 30.48052, 29.60492, 3.15794, 1.81185, 3.86731
1512.00000, 31.21601, 30.56624, 3.28543, 1.93330, 3.96983
1584.00000, 32.26978, 31.53116, 3.37295, 2.01844, 4.01486
1656.00000, 33.47753, 32.48850, 3.44213, 2.08675, 4.03074
1728.00000, 34.60353, 33.46535, 3.51243, 2.15707, 4.05039
1800.00000, 35.38664, 34.42570, 3.63615, 2.28302, 4.15405
1872.00000, 35.99842, 35.33325, 3.77535, 2.42805, 4.28767
1944.00000, 37.25882, 36.40189, 3.87679, 2.53592, 4.34103
2016.00000, 38.27939, 37.45721, 3.95595, 2.62135, 4.38073
2088.00000, 39.62065, 38.39246, 3.99279, 2.66150, 4.34444
2160.00000, 40.54045, 39.36550, 4.08908, 2.76752, 4.41133
2232.00000, 41.46744, 40.26702, 4.16698, 2.85448, 4.45364
2304.00000, 42.16638, 41.16869, 4.28734, 2.99088, 4.57060
2376.00000, 43.15928, 42.01146, 4.33242, 3.04261, 4.56354
2448.00000, 43.72279, 42.75776, 4.44022, 3.16770, 4.67586
2520.00000, 44.16657, 43.44978, 4.52902, 3.27222, 4.76947
2592.00000, 45.01178, 44.03320, 4.53014, 3.27355, 4.70948
2664.00000, 44.93077, 44.62428, 4.68649, 3.46079, 4.93904
2736.00000, 45.80369, 45.21145, 4.70069, 3.47801, 4.89918
2808.00000, 46.42891, 45.68797, 4.71382, 3.49395, 4.87242
2880.00000, 46.84833, 46.07878, 4.73588, 3.52080, 4.87308
2952.00000, 46.56091, 46.38735, 4.83981, 3.64834, 5.04699
3024.00000, 46.96915, 46.65051, 4.84598, 3.65597, 5.02724
3096.00000, 47.49161, 46.88960, 4.83099, 3.63745, 4.96438
3168.00000, 47.81815, 47.11535, 4.83366, 3.64074, 4.93810
3240.00000, 48.05279, 47.32123, 4.84426, 3.65385, 4.93065
3312.00000, 48.14154, 47.52392, 4.88385, 3.70292, 4.98202
3384.00000, 48.40099, 47.70476, 4.88677, 3.70655, 4.96252
3456.00000, 48.66420, 47.88113, 4.88712, 3.70699, 4.93654
3528.00000, 48.81804, 48.05185, 4.90603, 3.73054, 4.95023
3600.00000, 48.33537, 48.21431, 5.01205, 3.86361, 5.14902
3672.00000, 48.56224, 48.36496, 5.01163, 3.86307, 5.12784
3744.00000, 49.11767, 48.40752, 4.92359, 3.75244, 4.94454
3816.00000, 48.83425, 48.35652, 4.95320, 3.78952, 5.01482
3888.00000, 48.84083, 48.25840, 4.92075, 3.74889, 4.96431
3960.00000, 48.45630, 48.13365, 4.94892, 3.78415, 5.04188
4032.00000, 48.04258, 47.99205, 4.96671, 3.80647, 5.10188
4104.00000, 47.76334, 47.83987, 4.96296, 3.80176, 5.12034
4176.00000, 47.86234, 47.63102, 4.91854, 3.74613, 5.05005
4248.00000, 47.99848, 47.49238, 4.85754, 3.67027, 4.94869
4320.00000, 47.57076, 47.33722, 4.87911, 3.69703, 5.01794
4392.00000, 47.27711, 47.16274, 4.88295, 3.70180, 5.04800
4464.00000, 47.23937, 46.97355, 4.83764, 3.64566, 4.98606
4536.00000, 47.03385, 46.78720, 4.82268, 3.62719, 4.98216
4608.00000, 47.10518, 46.58076, 4.76013, 3.55039, 4.88518
4680.00000, 47.03528, 46.37485, 4.71453, 3.49481, 4.82373
4752.00000, 46.63767, 46.16884, 4.73116, 3.51504, 4.88277
4824.00000, 46.44290, 45.96285, 4.70612, 3.48459, 4.86039
4896.00000, 45.99562, 45.74296, 4.71739, 3.49828, 4.91026
4968.00000, 45.90502, 45.50400, 4.67021, 3.44112, 4.84803
5040.00000, 45.79009, 45.24673, 4.61760, 3.37779, 4.77880
5112.00000, 44.93695, 44.96774, 4.66590, 3.43591, 4.91000
5184.00000, 44.99797, 44.66556, 4.59005, 3.34482, 4.79816
5256.00000, 44.66022, 44.33398, 4.55554, 3.30369, 4.77296
5328.00000, 44.52679, 43.96627, 4.47773, 3.21169, 4.67062
5400.00000, 43.67437, 43.55701, 4.49821, 3.23581, 4.76056
5472.00000, 43.36422, 43.09439, 4.42831, 3.15379, 4.68427
5544.00000, 42.74904, 42.60299, 4.39112, 3.11048, 4.67436
5616.00000, 42.39748, 42.07443, 4.32409, 3.03303, 4.60570
5688.00000, 41.47123, 41.51099, 4.27888, 2.98122, 4.60532
5760.00000, 41.03814, 40.78347, 4.16945, 2.85725, 4.48621
5832.00000, 40.41068, 40.03684, 4.04686, 2.72083, 4.36312
5904.00000, 38.99552, 39.06543, 4.01303, 2.68364, 4.41052
5976.00000, 37.52924, 38.02392, 3.94710, 2.61175, 4.41286
6048.00000, 36.76810, 36.85139, 3.76973, 2.42212, 4.23711
6120.00000, 35.02888, 35.56422, 3.68790, 2.33654, 4.23477
6192.00000, 34.43871, 34.27381, 3.47369, 2.11821, 4.01382
6264.00000, 32.84459, 32.85914, 3.34566, 1.99174, 3.95335
6336.00000, 31.20841, 31.48464, 3.25031, 1.89953, 3.93128
6408.00000, 29.44449, 30.04759, 3.12581, 1.77300, 3.87148
6480.00000, 28.30699, 28.45925, 3.04620, 1.56174, 3.66042
6552.00000, 26.52974, 26.78967, 2.98429, 1.41172, 3.54783
6624.00000, 23.53344, 24.95556, 2.85588, 1.24888, 3.41533
6696.00000, 21.45655, 22.59024, 2.81390, 1.21243, 3.46319
6768.00000, 19.52579, 20.57863, 2.68466, 1.10362, 3.36649
6840.00000, 18.31560, 19.04672, 2.64497, 1.07123, 3.37707
6912.00000, 17.23064, 17.82610, 2.62870, 1.05809, 3.41552
6984.00000, 16.05194, 16.65357, 2.54134, 0.98894, 3.35777
7056.00000, 14.92180, 15.51599, 2.48070, 0.94231, 3.33975
7128.00000, 13.79359, 14.35746, 2.41822, 0.89544, 3.32307
7200.00000, 12.70055, 13.21960, 2.33758, 0.83671, 3.28418

[stoiver]

@mikilterribile Michele, I notice that the enquiry points are associated with triangles that are also exchange triangles. This might imply a bug in the way we calculate the location of the enquiry points. But for now, could you try manually locating the enquiry points some distance upstream of the riverwall and down stream of the outfall.

I'm thinking that the water is removed from the exchange triangles and thus affects the calculation of the flow based on state at the enquiry point.

[stoiver]

@mikilterribile Michele, it is also interesting that your exchange lines just intersect with 3 and 1 triangle. That probably means that a slight refinement of the mesh near these exchange lines would be sensible.

Could you provide a plot showing the triangles in this region.

[mikilterribile]

This are the images of the mesh zooming in progressivly. I have tryed moving the in/out points/exchange lines but nothing has changed. Now I try to refine the mesh in this regions, and I will see... Thank you for your kind support, Michele.

[mikilterribile]

I have tryed to create a high resolution zone around the inlet-outlet of the culvert but it leds to the same results...