Modelling ice tributary in 2D flowline case

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

Modelling ice tributary in 2D flowline case

Post by joeatodd » 09 Apr 2012, 16:56

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
Site Admin
Posts: 3264
Joined: 22 Aug 2009, 11:57
Antispam: Yes
Location: Espoo, Finland
Contact:

Re: Modelling ice tributary in 2D flowline case

Post by raback » 09 Apr 2012, 21:55

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

Post by joeatodd » 10 Apr 2012, 19:19

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
  Temp Load = Real 1.23
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

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

Constants

End

Simulation
  Coordinate System = "Cartesian 2D" 
  Simulation Type = Steady-State
  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
  Gradient = Real -1.0E00
  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

  Gradient = Real 1.0E00

!  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
Site Admin
Posts: 3264
Joined: 22 Aug 2009, 11:57
Antispam: Yes
Location: Espoo, Finland
Contact:

Re: Modelling ice tributary in 2D flowline case

Post by raback » 10 Apr 2012, 21:29

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

Post by joeatodd » 11 Apr 2012, 00:23

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

Post by joeatodd » 11 Apr 2012, 03:39

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

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

Re: Modelling ice tributary in 2D flowline case

Post by raback » 11 Apr 2012, 14:43

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

Post by joeatodd » 11 Apr 2012, 15:36

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

Joe

Post Reply