## Modelling ice tributary in 2D flowline case

Numerical methods and mathematical models of Elmer
joeatodd
Posts: 34
Joined: 02 Feb 2012, 18:49
Antispam: Yes

### Modelling ice tributary in 2D flowline case

Hi,

I've been working with Elmer this year to develop a 2D flowline model of a glacier in Greenland. I've hit a bit of a stumbling block with regards to modelling the additional mass input at the point where a tributary glacier joins the main channel.

Is there a way to artificially create a 'mass source' by adding mass at individual nodes in my 2D mesh to account for the confluence of this tributary?

I thought about simply adding an accumulation flux at the surface BC, but I don't think this would produce realistic results.

Any help is greatly appreciated,

Thanks,

Joe

raback
Posts: 3321
Joined: 22 Aug 2009, 11:57
Antispam: Yes
Location: Espoo, Finland
Contact:

### Re: Modelling ice tributary in 2D flowline case

Hi Joe

Well, if you want a nodal load on node then basically you can do this by the following construct in 2D. It sets a nodal load at the node closest to origin for the continuity equation.

Code: Select all

``````Boundary Condition 5
Name = "Nodal Mass Source"
Target Coordinates(1,2) = 0.0 0.0
Flow Solution Loads 3 = Real 1.23
End
``````
-Peter

joeatodd
Posts: 34
Joined: 02 Feb 2012, 18:49
Antispam: Yes

### Re: Modelling ice tributary in 2D flowline case

Hi Peter,

Thanks for your quick reply. I'm having some trouble implementing your proposed solution, unfortunately. Would you mind having a quick look over my SIF to see if you can spot any obvious mistakes?

I've been trying to pinpoint the problem. I'm able to implement a nodal temperature load using

Code: Select all

``````Boundary Condition 5
Name = "Nodal Mass Source"
Target Coordinates(1,2) = 70000.0 0.0
End ``````
but, for some reason, it seems to ignore "Flow Solution Loads 3"

Code: Select all

``````\$Rho = 9.03760022565934E-019
!\$RhoW = 1.00417780285104E-018
\$RhoW = 1.03430313693657E-018
\$g = 9.85098424596868E015

!For Glens
\$GlenN = 3.0

!For Width
\$N = 3.0
\$Bforce = 6.0E9

!For basal slip
\$As  = 4.0E-000
\$q   = 2.0
\$C   = 0.42
\$ut0 = 50.0
\$m   = 3.0

!Piezometric stuff
\$L = 109550.0
\$Piez = 0.015

Mesh DB "." "mesh"
Include Path ""
Results Directory "." "Results"
End

Constants

End

Simulation
Coordinate System = "Cartesian 2D"
Coordinate Mapping(3) = 1 2 3
Steady State Min Iterations = 1
Steady State Max Iterations = 5
Set Dirichlet BCs By BC Numbering = Logical True
max output level = 3
End

Body 1
name = All
Equation = 1
Material = 1
Body force = 1
Initial Condition = 1
End

Material 1

Density = Real \$ Rho
Viscosity Model = Power Law
Viscosity Exponent = \$ 1.0/3.0
Critical Shear Rate = \$ 3.0e-09*31556926.0
Viscosity = Variable Temp
Real MATC "glen(tx) * 31556926.0^(-1.0/3.0)*1.0E-06"
! the heat capacity as a MATC function of temperature itself
!-----------------------------------------------------------
Temp Heat Capacity = Variable Temp
Real MATC "capacity(tx)*(31556926.0)^(2.0)"
! the heat conductivity as a MATC function of temperature itself
!--------------------------------------------------------------
Temp Heat Conductivity = Variable Temp
Real MATC "conductivity(tx) * (31556926) * 1.0E-06"
! Upper limit - pressure melting point
!  as a MATC function of the pressure (what else?)
!-------------------------------------------------
Temp Upper Limit = Variable Pressure
Real MATC "pressuremeltingpoint(tx)"
! lower limit (to be save) as 0 K
!--------------------------------
Temp Lower Limit = Real 0.0

Cauchy = Logical False

End

!!!!! SOLVERS

Solver 1
Equation = "Depth"
Exec Solver = "Before TimeStep"

Procedure = File "src/DepthSolver" "DepthSolver"
Variable = String "Depth"
Variable DOFs = 1
Linear System Direct Method = "Mumps"
Linear System Solver = Iterative
Linear System Max Iterations = 300
Linear System Iterative Method = "BiCGStab"
Linear System Preconditioning = ILU0
Linear System Convergence Tolerance = Real 1.0E-6
Linear System Abort Not Converged = False
Linear System Residual Output = 0
End

Solver 2
Equation = "Elevation"
Exec Solver = "Before TimeStep"

Procedure = File "src/DepthSolver" "DepthSolver"
Variable = String "Elevation"
Variable DOFs = 1

!  Linear System Solver = "Direct"
Linear System Direct Method = "Mumps"
!  Linear System Direct Method = umfpack

Linear System Solver = Iterative
Linear System Max Iterations = 300
Linear System Iterative Method = "BiCGStab"
Linear System Preconditioning = ILU0
Linear System Convergence Tolerance = Real 1.0E-6
Linear System Abort Not Converged = False
Linear System Residual Output = 0
End

Solver 3
Exec Solver = "Before Simulation"
Equation = "Normal vector"
Variable = "Normal Vector"
! in 3dimensional simulations we have 3 entries
Variable DOFs = 2
!NB: does not need to actually solve a matrix
!    hence no BW optimization needed
Optimize Bandwidth = Logical False
Procedure = "src/ComputeNormal" "ComputeNormalSolver"
! if set to True, all boundary normals would be computed by default
ComputeAll = Logical True
End

Solver 4
Equation = "Navier-Stokes"
Procedure = "FlowSolve" "FlowSolver"
Variable = Flow Solution[Velocity:2 Pressure:1]
Exec Solver = Always
Flow model = String "Stokes"
Stabilize = True
! Bubbles = True
Linear System Direct Method = "MUMPS"

Linear System Solver = "Iterative"
Linear System Iterative Method = "BiCGStabl"
BiCGStabl Polynomial Degree = 4
Linear System Max Iterations = 500
Linear System Convergence Tolerance = Real 1.0E-6
Linear System Abort Not Converged = False
Linear System Preconditioning = "ILU1"
Linear System Residual Output = 0

Nonlinear System Max Iterations = 50
Nonlinear System Convergence Tolerance = 1.0E-5
Nonlinear System Newton After Iterations = 3
Nonlinear System Newton After Tolerance = 1.0e-4
Nonlinear System Relaxation Factor = 0.999999 !0.999999
Steady State Relaxation Factor = Real 0.9
Steady State Convergence Tolerance = 1.0d-3

Exported Variable 1 = Width
Exported Variable 2 = KCoef
Exported Variable 3 = Flow BodyForce 1
Exported Variable 4 = Draggy
Exported Variable 5 = Slip Coefficient 2
Exported Variable 6 = Overburden
Update Exported Variables = Logical True
Nonlinear Update Exported Variables = Logical True

End

Solver 5
Equation = String "StressSolver"
Procedure =  File "src/ComputeDevStressNS" "ComputeDevStress"
! this is just a dummy, hence no output is needed
!-----------------------------------------------------------------------
Variable = -nooutput "Sij"
Variable DOFs = 1
! the name of the variable containing the flow solution (U,V,W,Pressure)
!-----------------------------------------------------------------------
Flow Solver Name = String "Flow Solution"
Exported Variable 1 = "Stress" ! [Sxx, Syy, Szz, Sxy] in 2D
! [Sxx, Syy, Szz, Sxy, Syz, Szx] in 3D
Exported Variable 1 DOFs = 4   ! 4 in 2D, 6 in 3D
Linear System Solver = "Iterative"
Linear System Iterative Method = "BiCGStab"
Linear System Max Iterations = 300
Linear System Convergence Tolerance = 1.0E-09
Linear System Abort Not Converged = True
Linear System Preconditioning = "ILU0"
Linear System Residual Output = 0
Cauchy = Logical False
End

Solver 6
!  Exec Solver = "After TimeStep"
Exec Solver = "After Saving"
!  Exec Solver = "After All"

Procedure = "SaveData" "SaveLine"
Filename =  "surface.dat"
Optimize Node Ordering = Logical False
File Append = Logical False
Variable 1 = Coordinate 1
Variable 2 = Velocity 1
Variable 3 = Slip Coefficient 2
End

Solver 7
Equation = "Result Output"
!  Exec Solver = "After Saving"
Exec Solver = "After All"

Procedure = "ResultOutputSolve" "ResultOutputSolver"
Output File Name = "results"

Output Format = String "vtu"
Vector Field 1 = String "Velocity"
Vector Field 2 = String "Stress"
Scalar Field 1 = String "Depth"
Scalar Field 2 = String "Elevation"
Scalar Field 3 = String "Pressure"
Scalar Field 4 = String "Temp"
Scalar Field 5 = String "Temp Homologous"
Scalar Field 6 = String "Temp Residual"
Scalar Field 7 = String "Velocity 1"
Scalar Field 8 = String "Width"
Scalar Field 9 = String "KCoef"
Scalar Field 10 = String "Flow BodyForce 1"
Scalar Field 11 = String "Draggy"
End

Solver 8
Equation = String "Homologous Temperature Equation"
Procedure =  File "src/TemperateIce" "TemperateIceSolver"
! Before Linsolve = "EliminateDirichlet" "EliminateDirichlet"
Variable = String "Temp"
Variable DOFs = 1
Linear System Solver = "Iterative"
Linear System Iterative Method = "BiCGStab"
Linear System Max Iterations = 10000
Linear System Convergence Tolerance = 1.0E-06
Linear System Abort Not Converged = False
Linear System Preconditioning = "ILU4"
Linear System Residual Output = 1
Steady State Convergence Tolerance = 1.0E-04
Nonlinear System Convergence Tolerance = 1.0E-05
Nonlinear System Max Iterations = 50
Nonlinear System Relaxation Factor = Real 0.999 !0.9999
Steady State Relaxation Factor = Real 0.9
Apply Dirichlet = Logical True
Stabilize = True
Exported Variable 1 = String "Temp Homologous"
Exported Variable 1 DOFs = 1
Exported Variable 2 = String "Temp Residual"
Exported Variable 2 DOFs = 1
End

!!!!! EQUATION

Equation 1
Active Solvers (8) = 1 2 3 8 4 5 6 7
Flow Solution Name = String "Flow Solution"
Convection = Computed
NS Convect = True
End

!!!!! BOUNDARY CONDITIONS

Boundary Condition 1
name = Base
Target Boundaries = 2

Flow Force BC = Logical True
Normal-Tangential Velocity = Logical True
Velocity 1 = 0.0

External Pressure = Variable Coordinate 1, Coordinate 2
Real MATC "getPw(tx(0),tx(1))"
Slip Coefficient 2 = Variable Coordinate 1
Real Procedure "./src/USF_Sliding" "Friction_jgr2007"

Overburden = Variable Depth
Real MATC "tx*Rho*g"

Friction Law Sliding Coefficient = Real \$As
Friction Law Post-Peak Exponent = Real \$q
Friction Law Maximum Value = Real \$C
Friction Law Linear Velocity = Real \$ut0
Friction Law PowerLaw Exponent = Real \$m

Temp Flux BC = Logical True
Temp Heat Flux = Real 55.0E-03

Elevation = Real 0.0
Save Line = Logical True
End

Boundary Condition 2
name = Surface
Target Boundaries = 3
Temp = Real 253.15
Depth = Real 0.0
Save Line = Logical True
End

Boundary Condition 3
name = Right
Target Boundaries = 1

Flow Force BC = Logical True
External Pressure = Variable Coordinate 2
Real MATC "if(tx>10.0){0.0}else{if(tx<0.0){RhoW*g*tx}else{0.0}}"
! Real MATC "if(tx>0.0){0.0}else{RhoW*g*tx}"
!  External Pressure = Variable Coordinate 2, Depth, Elevation
!    Real MATC "-Rho*g*(tx(1))+0.5*Rho*g*(1.0-Rho/RhoW)*(tx(1)+tx(2))"

Save Line = Logical True
End

Boundary Condition 4
name = Left
Target Boundaries = 4

Flow Force BC = Logical True
External Pressure = Variable Depth, Elevation
Real MATC "-0.9635 * Rho*g*tx(0)"
!Real MATC "-1.0E-5 * tx(1)/(tx(0)+tx(1))"

Save Line = Logical True
End

Boundary Condition 5
Name = "Nodal Mass Source"
Target Coordinates(1,2) = 70000.0 0.0
Flow Solution Loads 3 = Real 5.0E10
End

Body Force 1

Width = Variable Coordinate 1
Real
include Halfwidth.dat
End

KCoef = Variable Width, Temp
Real MATC "((N+1)^(1/N))/((tx(0)^((N+1)/N))*(2*getArrheniusFactor(tx(1))^(1/N)))"

Draggy = Variable Velocity 1
Real MATC "if(tx==0){0.0}else{tx * abs(tx)^((1.0/N)-1.0)}"

Flow BodyForce 1 = Variable KCoef, Draggy
Real MATC "Bforce*-tx(0)*tx(1)"
!!2.8E9 works.
Flow BodyForce 2 = Real \$ -g
Temp Volume Source = Real 1.0E00
End

Initial Condition 1

Temp = Real 263.15

End

!Functions

!! pressuremeltingpoint (in SI units)
\$ function pressuremeltingpoint(PIN) {\
P = PIN;\
if (P<0.0) P=0.0;\
beta=9.8E-08*1.0E06;\
_pressuremeltingpoint=273.15-(beta*P);\
}

!! in SI units, input in Kelvin
\$ function conductivity(T)  { _conductivity=9.828*exp(-5.7E-03*T)}

!! in SI units, input in Kelvin
\$ function capacity(T) { _capacity=146.3+(7.253*T)}

!! Glen's flow law
!-----------------
\$ function glen(Th) \
import GlenN {\
EF = 1.0;\
AF = getArrheniusFactor(Th);\
_glen = (2.0*EF*AF)^(-1.0/GlenN);\
}

!! Arrhenius factor needed by glen
!! (in SI units)
!---------------------------------
\$ function getArrheniusFactor(Th){ \
Th1 = Th - 273.15;\
if (Th1<-10) {_getArrheniusFactor=3.985E-13 * exp( -60.0E03/(8.314 * (273.15 + Th1)));}\
else {\
if (Th1>0) _getArrheniusFactor=1.916E03 * exp( -139.0E03/(8.314 *  (273.15)));\
else _getArrheniusFactor=1.916E03 * exp( -139.0E03/(8.314 *  (273.15 + Th1)));}\
}

!! Basal Water Pressure - calculates based on a given piezometric gradient.
!-----------------
!\$ function getPw(C1,Depth) \
!  import L,Base,RhoW,g,Divide { \
! x = (L - C1)/L;\
!y = Depth * (Base - ((Base - Divide)*x));\
!if (y>0) {_getPw = RhoW*g*y;}\
! else {_getPw = 0.0;}\
!}

\$ function getPw(C1,C2) \
import L,Piez,RhoW,g { \
x = L - C1;\
y = -C2 + (x*Piez);\
if (y>0) {_getPw = RhoW*g*y;}\
else {_getPw = 0.0;}\
}
``````
Thanks again for your help. Much appreciated!

Joe

raback
Posts: 3321
Joined: 22 Aug 2009, 11:57
Antispam: Yes
Location: Espoo, Finland
Contact:

### Re: Modelling ice tributary in 2D flowline case

Hi

Ok, my fault. I mixed two different conventions: the name of computed nodal loads and the name of given nodal loads. These should rather be same, but unfortunately they are not Both features are rather rarely used so it is not a major problem but still a slightly disturbing. So try

Code: Select all

``````Flow Solution 3 Load = Real 5.0e10
``````
-Peter

joeatodd
Posts: 34
Joined: 02 Feb 2012, 18:49
Antispam: Yes

### Re: Modelling ice tributary in 2D flowline case

Hi Peter,

Still no luck I'm afraid; am I maybe missing a keyword in the NS Solver section or something like that?

Cheers,

Joe

joeatodd
Posts: 34
Joined: 02 Feb 2012, 18:49
Antispam: Yes

### Re: Modelling ice tributary in 2D flowline case

Interestingly, Displacement 1 Load doesn't seem to do anything either. Is there something wrong with my NS Equation settings?

raback
Posts: 3321
Joined: 22 Aug 2009, 11:57
Antispam: Yes
Location: Espoo, Finland
Contact:

### Re: Modelling ice tributary in 2D flowline case

Hi

Ok, if the vector variable has components that have been named differently, as is the case with FlowSolve, then the name of the Load is derived from that name. I.e. use

Code: Select all

``````Pressure Load = Real 5.0e10
``````
There is a new consistency test for this as "pointload2". Sorry for the confusing instructions.

-Peter

joeatodd
Posts: 34
Joined: 02 Feb 2012, 18:49
Antispam: Yes

### Re: Modelling ice tributary in 2D flowline case

Yep that did it! Thanks very much for your help, Peter.

Joe