Fluidity

will bedload sediments go into suspension load?

Asked by liu chin chi

when use "reentrainment " (the boundary condition sediment ).
will the bedload part (decrease in bedloadVolumeFraction) go into suspension load(increase in sediment concentration) ?

Question information

Language:
English Edit question
Status:
Solved
For:
Fluidity Edit question
Assignee:
No assignee Edit question
Solved by:
Samuel Parkinson
Solved:
Last query:
Last reply:

This question was reopened

Revision history for this message
Samuel Parkinson (s-parkinson11) said : 		#1

Hello Liu,

Yes, eroded sediment will go back into suspension. The "reentrainment" boundary is implemented as a Neumann condition, or flux of sediment, which varies depending upon the bed shear stress and properties of the sediment in the bed that is being eroded.

Cheers,

Sam

Revision history for this message
liu chin chi (caite8001) said : 		#2

Thanks Samuel Parkinson, that solved my question.

Revision history for this message
liu chin chi (caite8001) said : 		#3

In the 2D domain of model, X was set from 0 to 4, and Y was set from 0 to 0.5.
Then we set some suspended sediments in the initial condition.
As time passed, the sediment settled to the bottom, and the thickness of bedload became thicker.
As suspended sediments settled totally to the bottom, the both side boundaries were given vertical average velocity = 0.1 m s-1. in X dirction.then the thickness of bedload drcrease close to zero, but suspended sediments don't increase.

the case seems eroded sediment doesn't go back to suspension again.

Would you please give me some comments?
Thank you very much.

Revision history for this message
Samuel Parkinson (s-parkinson11) said : 		#4

Hi Liu,

Can you provide a copy of your flml file? This is quite hard for me to answer without having a full understanding of your model set up.

Kind regards,

Sam

Revision history for this message
liu chin chi (caite8001) said : 		#5

Sir, I've posted the setup in flml file below.

I'm totally a newbie, your help and patience is very much appreciated.

---------------------------------------------------------------------------------------------------------------

<?xml version='1.0' encoding='utf-8'?>
<fluidity_options>
  <simulation_name>
    <string_value lines="1">iwsink2d</string_value>
    <comment>2D lock exchange</comment>
  </simulation_name>
  <problem_type>
    <string_value lines="1">oceans</string_value>
  </problem_type>
  <geometry>
    <dimension>
      <integer_value rank="0">2</integer_value>
    </dimension>
    <mesh name="CoordinateMesh">
      <from_file file_name="iwsink2d">
        <format name="triangle"/>
        <stat>
          <include_in_stat/>
        </stat>
      </from_file>
    </mesh>
    <mesh name="VelocityMesh">
      <from_mesh>
        <mesh name="CoordinateMesh"/>
        <stat>
          <exclude_from_stat/>
        </stat>
      </from_mesh>
    </mesh>
    <mesh name="PressureMesh">
      <from_mesh>
        <mesh name="CoordinateMesh"/>
        <stat>
          <exclude_from_stat/>
        </stat>
      </from_mesh>
    </mesh>
    <quadrature>
      <degree>
        <integer_value rank="0">3</integer_value>
      </degree>
      <surface_degree>
        <integer_value rank="0">3</integer_value>
      </surface_degree>
    </quadrature>
  </geometry>
  <io>
    <dump_format>
      <string_value>vtk</string_value>
    </dump_format>
    <dump_period>
      <constant>
        <real_value rank="0">1.0</real_value>
      </constant>
    </dump_period>
    <output_mesh name="VelocityMesh"/>
    <stat>
      <output_at_start/>
    </stat>
  </io>
  <timestepping>
    <current_time>
      <real_value rank="0">0.0</real_value>
    </current_time>
    <timestep>
      <real_value rank="0">0.05</real_value>
    </timestep>
    <finish_time>
      <real_value rank="0">20.0</real_value>
      <comment>10.0</comment>
    </finish_time>
    <nonlinear_iterations>
      <integer_value rank="0">2</integer_value>
    </nonlinear_iterations>
  </timestepping>
  <physical_parameters>
    <gravity>
      <magnitude>
        <real_value rank="0">9.81</real_value>
      </magnitude>
      <vector_field name="GravityDirection" rank="1">
        <prescribed>
          <mesh name="CoordinateMesh"/>
          <value name="WholeMesh">
            <constant>
              <real_value shape="2" dim1="dim" rank="1">0.0 -1.0</real_value>
            </constant>
          </value>
          <output/>
          <stat>
            <include_in_stat/>
          </stat>
          <detectors>
            <exclude_from_detectors/>
          </detectors>
        </prescribed>
      </vector_field>
    </gravity>
  </physical_parameters>
  <material_phase name="fluid">
    <equation_of_state>
      <fluids>
        <linear>
          <reference_density>
            <real_value rank="0">1.0</real_value>
          </reference_density>
          <subtract_out_hydrostatic_level/>
        </linear>
      </fluids>
    </equation_of_state>
    <scalar_field name="Pressure" rank="0">
      <prognostic>
        <mesh name="PressureMesh"/>
        <spatial_discretisation>
          <continuous_galerkin/>
        </spatial_discretisation>
        <reference_node>
          <integer_value rank="0">10</integer_value>
        </reference_node>
        <scheme>
          <poisson_pressure_solution>
            <string_value lines="1">only first timestep</string_value>
          </poisson_pressure_solution>
          <use_projection_method/>
        </scheme>
        <solver>
          <iterative_method name="cg"/>
          <preconditioner name="mg"/>
          <relative_error>
            <real_value rank="0">1.0e-6</real_value>
          </relative_error>
          <max_iterations>
            <integer_value rank="0">3000</integer_value>
          </max_iterations>
          <never_ignore_solver_failures/>
          <diagnostics>
            <monitors/>
          </diagnostics>
        </solver>
        <output/>
        <stat>
          <exclude_from_stat/>
        </stat>
        <convergence>
          <include_in_convergence/>
        </convergence>
        <detectors>
          <exclude_from_detectors/>
        </detectors>
        <steady_state>
          <include_in_steady_state/>
        </steady_state>
        <no_interpolation/>
      </prognostic>
    </scalar_field>
    <scalar_field name="Density" rank="0">
      <diagnostic>
        <algorithm name="Internal" material_phase_support="multiple"/>
        <mesh name="VelocityMesh"/>
        <output/>
        <stat/>
        <convergence>
          <include_in_convergence/>
        </convergence>
        <detectors>
          <include_in_detectors/>
        </detectors>
        <steady_state>
          <include_in_steady_state/>
        </steady_state>
      </diagnostic>
    </scalar_field>
    <vector_field name="Velocity" rank="1">
      <prognostic>
        <mesh name="VelocityMesh"/>
        <equation name="Boussinesq"/>
        <spatial_discretisation>
          <continuous_galerkin>
            <stabilisation>
              <no_stabilisation/>
            </stabilisation>
            <mass_terms>
              <lump_mass_matrix/>
            </mass_terms>
            <advection_terms/>
            <stress_terms>
              <tensor_form/>
            </stress_terms>
          </continuous_galerkin>
          <conservative_advection>
            <real_value rank="0">0.0</real_value>
          </conservative_advection>
        </spatial_discretisation>
        <temporal_discretisation>
          <theta>
            <real_value rank="0">0.5</real_value>
          </theta>
          <relaxation>
            <real_value rank="0">0.5</real_value>
          </relaxation>
        </temporal_discretisation>
        <solver>
          <iterative_method name="cg"/>
          <preconditioner name="sor"/>
          <relative_error>
            <real_value rank="0">1.0e-6</real_value>
            <comment>1.0e-6</comment>
          </relative_error>
          <max_iterations>
            <integer_value rank="0">3000</integer_value>
          </max_iterations>
          <never_ignore_solver_failures/>
          <diagnostics>
            <monitors/>
          </diagnostics>
        </solver>
        <initial_condition name="WholeMesh">
          <constant>
            <real_value shape="2" dim1="dim" rank="1">0.0 0.0</real_value>
          </constant>
        </initial_condition>
        <boundary_conditions name="top">
          <surface_ids>
            <integer_value shape="1" rank="1">9</integer_value>
          </surface_ids>
          <type name="dirichlet">
            <align_bc_with_cartesian>
              <y_component>
                <constant>
                  <real_value rank="0">0.0</real_value>
                </constant>
              </y_component>
            </align_bc_with_cartesian>
          </type>
        </boundary_conditions>
        <boundary_conditions name="sides">
          <surface_ids>
            <integer_value shape="2" rank="1">8 10</integer_value>
          </surface_ids>
          <type name="dirichlet">
            <align_bc_with_cartesian>
              <x_component>
                <python>
                  <string_value lines="20" type="code" language="python">def val(X,t):
 if(t&lt;10):
  return 0.025*(t-0)/10
 else:
  return 0.025</string_value>
                  <comment># elif(t&gt;60.0):
# return 0.2*(t-60)*2/60</comment>
                </python>
              </x_component>
              <y_component>
                <constant>
                  <real_value rank="0">0</real_value>
                </constant>
              </y_component>
            </align_bc_with_cartesian>
          </type>
        </boundary_conditions>
        <boundary_conditions name="bottom">
          <surface_ids>
            <integer_value shape="1" rank="1">7</integer_value>
          </surface_ids>
          <type name="dirichlet">
            <align_bc_with_cartesian>
              <y_component>
                <constant>
                  <real_value rank="0">0</real_value>
                </constant>
              </y_component>
            </align_bc_with_cartesian>
          </type>
        </boundary_conditions>
        <output/>
        <stat>
          <exclude_from_stat/>
          <previous_time_step>
            <exclude_from_stat/>
          </previous_time_step>
          <nonlinear_field>
            <exclude_from_stat/>
          </nonlinear_field>
        </stat>
        <convergence>
          <include_in_convergence/>
        </convergence>
        <detectors>
          <include_in_detectors/>
        </detectors>
        <steady_state>
          <include_in_steady_state/>
        </steady_state>
        <adaptivity_options>
          <absolute_measure>
            <vector_field name="InterpolationErrorBound" rank="1">
              <prescribed>
                <value name="WholeMesh">
                  <constant>
                    <real_value shape="2" dim1="dim" rank="1">0.001 0.001</real_value>
                    <comment>def val(X, t):

# exponential function of form p = Ae^(aq) + B,
# with q = q_1, p=p_1 and q = q_2, p=p_2
# giving p = p_1 + Aexp{aq_1}(exp{a(q-q_1)} -1)
# with A = (p_2-p_1)/(exp{aq_1}(exp{a(q_2-q_1)} -1))
# so p = p_1 + (p_2-p_1)*(exp{a(q-q_1)} -1) /(exp{a(q_2-q_1)} -1)
# let C = (p_2-p_1)/(exp{a(q_2-q_1)} -1)

  import math
  a = 100.0
  q_1 = 0.0
  p_1 = 0.00001
  q_2 = 0.05
  p_2 = 0.001

  if(X[1]&amp;amp;lt;=0.05):
    C = (p_2 - p_1)/((math.exp(a*(q_2 - q_1)))-1)
    q = X[1]
    p = p_1 + C*((math.exp(a*(q - q_1)))-1)
    u_val = p
  else:
    C = (p_2 - p_1)/((math.exp(a*(q_2 - q_1)))-1)
    q = 0.1-X[1]
    p = p_1 + C*((math.exp(a*(q - q_1)))-1)
    u_val = p

  return [u_val,0.001]</comment>
                  </constant>
                </value>
                <output/>
                <stat>
                  <include_in_stat/>
                </stat>
                <detectors>
                  <exclude_from_detectors/>
                </detectors>
              </prescribed>
            </vector_field>
          </absolute_measure>
        </adaptivity_options>
        <consistent_interpolation/>
      </prognostic>
    </vector_field>
    <vector_field name="BedShearStress" rank="1">
      <diagnostic>
        <algorithm name="Internal" material_phase_support="multiple"/>
        <mesh name="VelocityMesh"/>
        <density>
          <real_value rank="0">1000</real_value>
        </density>
        <calculation_method>
          <drag_coefficient>
            <real_value rank="0">0.02</real_value>
          </drag_coefficient>
        </calculation_method>
        <output/>
        <stat>
          <include_in_stat/>
        </stat>
        <detectors>
          <include_in_detectors/>
        </detectors>
      </diagnostic>
    </vector_field>
    <sediment>
      <scalar_field name="qqq" rank="0">
        <prognostic>
          <mesh name="VelocityMesh"/>
          <equation name="AdvectionDiffusion"/>
          <spatial_discretisation>
            <control_volumes>
              <face_value name="FirstOrderUpwind"/>
              <diffusion_scheme name="ElementGradient"/>
            </control_volumes>
            <conservative_advection>
              <real_value rank="0">0.0</real_value>
            </conservative_advection>
          </spatial_discretisation>
          <temporal_discretisation>
            <theta>
              <real_value rank="0">0.5</real_value>
            </theta>
            <control_volumes/>
          </temporal_discretisation>
          <solver>
            <iterative_method name="gmres">
              <restart>
                <integer_value rank="0">30</integer_value>
              </restart>
            </iterative_method>
            <preconditioner name="sor"/>
            <relative_error>
              <real_value rank="0">1.0e-7</real_value>
            </relative_error>
            <max_iterations>
              <integer_value rank="0">3000</integer_value>
            </max_iterations>
            <never_ignore_solver_failures/>
            <diagnostics>
              <monitors/>
            </diagnostics>
          </solver>
          <initial_condition name="WholeMesh">
            <constant>
              <real_value rank="0">0.0</real_value>
            </constant>
          </initial_condition>
          <boundary_conditions name="bottom">
            <surface_ids>
              <integer_value shape="1" rank="1">7</integer_value>
            </surface_ids>
            <type name="sediment_reentrainment">
              <algorithm>Generic</algorithm>
              <viscosity>
                <real_value rank="0">0.000001</real_value>
                <comment>0.000001</comment>
              </viscosity>
            </type>
          </boundary_conditions>
          <boundary_conditions name="top">
            <surface_ids>
              <integer_value shape="1" rank="1">9</integer_value>
            </surface_ids>
            <type name="zero_flux"/>
          </boundary_conditions>
          <boundary_conditions name="sides">
            <surface_ids>
              <integer_value shape="2" rank="1">8 10</integer_value>
            </surface_ids>
            <type name="zero_flux"/>
          </boundary_conditions>
          <scalar_field name="SinkingVelocity" rank="0">
            <prescribed>
              <value name="WholeMesh">
                <constant>
                  <real_value rank="0">0.0001</real_value>
                </constant>
              </value>
              <output/>
              <stat/>
              <detectors>
                <exclude_from_detectors/>
              </detectors>
            </prescribed>
          </scalar_field>
          <output/>
          <stat/>
          <convergence>
            <include_in_convergence/>
          </convergence>
          <detectors>
            <include_in_detectors/>
          </detectors>
          <steady_state>
            <include_in_steady_state/>
          </steady_state>
          <consistent_interpolation/>
          <scalar_field name="Bedload" rank="0">
            <prognostic>
              <surface_ids>
                <integer_value shape="1" rank="1">7</integer_value>
              </surface_ids>
              <equation name="SedimentBedload"/>
              <spatial_discretisation>
                <continuous_galerkin>
                  <stabilisation>
                    <no_stabilisation/>
                  </stabilisation>
                  <advection_terms/>
                  <mass_terms/>
                </continuous_galerkin>
                <conservative_advection>
                  <real_value rank="0">0.0</real_value>
                </conservative_advection>
              </spatial_discretisation>
              <temporal_discretisation>
                <theta>
                  <real_value rank="0">0.5</real_value>
                </theta>
              </temporal_discretisation>
              <solver>
                <iterative_method name="gmres">
                  <restart>
                    <integer_value rank="0">30</integer_value>
                  </restart>
                </iterative_method>
                <preconditioner name="sor"/>
                <relative_error>
                  <real_value rank="0">1e-7</real_value>
                </relative_error>
                <max_iterations>
                  <integer_value rank="0">3000</integer_value>
                </max_iterations>
                <never_ignore_solver_failures/>
                <diagnostics>
                  <monitors/>
                </diagnostics>
              </solver>
              <initial_condition name="WholeMesh">
                <python>
                  <string_value lines="20" type="code" language="python">def val(X,t):
 if(X[1]&lt;0.025):
  return 0.001
 else:
  return 0.0</string_value>
                </python>
              </initial_condition>
              <output/>
              <stat/>
              <convergence>
                <include_in_convergence/>
              </convergence>
              <detectors>
                <include_in_detectors/>
              </detectors>
              <steady_state>
                <include_in_steady_state/>
              </steady_state>
              <consistent_interpolation/>
            </prognostic>
          </scalar_field>
          <scalar_field name="BedloadVolumeFraction" rank="0">
            <diagnostic>
              <output/>
              <stat/>
              <convergence>
                <include_in_convergence/>
              </convergence>
              <detectors>
                <include_in_detectors/>
              </detectors>
              <steady_state>
                <include_in_steady_state/>
              </steady_state>
            </diagnostic>
          </scalar_field>
          <scalar_field name="BedloadDepositRate" rank="0">
            <diagnostic>
              <output/>
              <stat/>
              <convergence>
                <include_in_convergence/>
              </convergence>
              <detectors>
                <include_in_detectors/>
              </detectors>
              <steady_state>
                <include_in_steady_state/>
              </steady_state>
            </diagnostic>
          </scalar_field>
          <scalar_field name="BedloadErosionRate" rank="0">
            <diagnostic>
              <output/>
              <stat/>
              <convergence>
                <include_in_convergence/>
              </convergence>
              <detectors>
                <include_in_detectors/>
              </detectors>
              <steady_state>
                <include_in_steady_state/>
              </steady_state>
            </diagnostic>
          </scalar_field>
          <submerged_specific_gravity>
            <real_value rank="0">1.65</real_value>
          </submerged_specific_gravity>
          <diameter>
            <real_value rank="0">0.0005</real_value>
            <comment>0.0005</comment>
          </diameter>
          <bed_porosity>
            <real_value rank="0">0.3</real_value>
            <comment>0.3</comment>
          </bed_porosity>
          <erodability>
            <real_value rank="0">1</real_value>
          </erodability>
          <critical_shear_stress>
            <real_value rank="0">0.00000001</real_value>
          </critical_shear_stress>
        </prognostic>
      </scalar_field>
      <scalar_field name="SedimentBedActiveLayerD50" rank="0">
        <diagnostic>
          <surface_ids>
            <integer_value shape="1" rank="1">7</integer_value>
          </surface_ids>
          <algorithm name="Internal" material_phase_support="multiple"/>
          <mesh name="VelocityMesh"/>
          <output/>
          <stat/>
          <convergence>
            <include_in_convergence/>
          </convergence>
          <detectors>
            <include_in_detectors/>
          </detectors>
          <steady_state>
            <include_in_steady_state/>
          </steady_state>
        </diagnostic>
      </scalar_field>
      <scalar_field name="SedimentBedActiveLayerSigma" rank="0">
        <diagnostic>
          <surface_ids>
            <integer_value shape="1" rank="1">7</integer_value>
          </surface_ids>
          <algorithm name="Internal" material_phase_support="multiple"/>
          <mesh name="VelocityMesh"/>
          <output/>
          <stat/>
          <convergence>
            <include_in_convergence/>
          </convergence>
          <detectors>
            <include_in_detectors/>
          </detectors>
          <steady_state>
            <include_in_steady_state/>
          </steady_state>
        </diagnostic>
      </scalar_field>
    </sediment>
  </material_phase>
</fluidity_options>

Revision history for this message
Best Samuel Parkinson (s-parkinson11) said : 		#6

Hi Liu,

You have no diffusivity in your sediment field. Reentrainment is set by using a Neumann bc. Without any diffusivity Neumann boundary conditions do not work, however the bedload field does not know this.

In the code the following procedure occurs:
1) a reentrainment rate is determined
2) the reentrainment boundary condition is set
3) the bed depth is modified
with no diffusivity the bed depth is still modified (3) but reentrainment (2) doesn't work and hence the sediment is lost. This is a bug. There should be a check to ensure that there is a non-zero diffusivity before reducing the sediment bed depth.

This is not a straight forward problem and is something I am looking at in detail at the moment. Literature states that diffusion of sediment is negligible for most grain sizes, which leads to the idea that there should be zero diffusion - hence no reentrainment. Additionally, if we think about what the boundary condition means physically we realise that reentrainment happens due to turbulence in the boundary layer generated by fluid flow over a rough bed surface (which we are not modelling).

My current view is that we need a turbulence model that includes some parameterisation of a boundary layer over a rough surface. I believe that such a turbulence model would provide some diffusion of sediment near the boundary. This type of turbulence model does not currently exist in Fluidity and so is not an option. I plan to start looking at developing this capability over the coming months.

Another option is to include a small diffusion of the sediment field throughout the domain that is lower than the numerical diffusion in your problem (i.e. 1e-10). The rate of sediment entrainment will be correct but I have found that near bed concentrations of sediment can become unbounded leading to instabilities. It is also dubious whether this approach will provide physically realistic results.

For my current work I have actually removed sediment reentrainment due to these issues and would recommend you think about whether you really need this capability and consider not using it. I am very interested to know the problem you are trying to tackle and any views you have on this.

A very long and complicated answer. Sorry about this, I am still learning how to apply this properly myself. Let me know anything that you don't understand.

Cheers,

Sam

Revision history for this message
liu chin chi (caite8001) said : 		#7

Sir,

It works as the diffusion option is opened.
I'd try to simulate an internal wave passing through a slope in a tank, and see if small sand waves will be formed.(Maybe resuspension of sediments will settle down in other area and won't take a large affection in original erosion area.)

Appreciate your comments very much.

Regards,

ChinChi

Revision history for this message
liu chin chi (caite8001) said : 		#8

Thanks Samuel Parkinson, that solved my question.