# Difference between revisions of "COMSOL activities"

(→Field equations) |
(→Field equations) |
||

(25 intermediate revisions by 6 users not shown) | |||

Line 20: | Line 20: | ||

====Field equations==== | ====Field equations==== | ||

This equation mode solves equations of the form | This equation mode solves equations of the form | ||

− | :<math>e_a\frac{\partial^2 H}{\partial t^2} + d_a \frac{\partial H}{\partial t} + \nabla \cdot \ | + | :<math>e_a\frac{\partial^2 H}{\partial t^2} + d_a \frac{\partial H}{\partial t} + \nabla \cdot \Gamma = F</math> |

Which is just what we want if we recognize that in our system <math>e_a</math>=0, <math>d_a</math>=1, <math>F</math>=<math>M</math>, and | Which is just what we want if we recognize that in our system <math>e_a</math>=0, <math>d_a</math>=1, <math>F</math>=<math>M</math>, and | ||

Line 40: | Line 40: | ||

==Shallow shelf approximation== | ==Shallow shelf approximation== | ||

+ | ===Field equations=== | ||

Now, consider the equations describing a flow that is vertically integrated. The equations are | Now, consider the equations describing a flow that is vertically integrated. The equations are | ||

:<math>\frac{\partial}{\partial x}\left ( 2 \eta H | :<math>\frac{\partial}{\partial x}\left ( 2 \eta H | ||

Line 45: | Line 46: | ||

+\frac{\partial}{\partial y}\left(\eta H\left( | +\frac{\partial}{\partial y}\left(\eta H\left( | ||

\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\right) | \frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\right) | ||

− | =\ | + | =\rho_w gH \frac{\partial s}{\partial x} |

</math> | </math> | ||

Line 53: | Line 54: | ||

+\frac{\partial}{\partial x}\left(\eta H\left( | +\frac{\partial}{\partial x}\left(\eta H\left( | ||

\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\right) | \frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\right) | ||

− | =\ | + | =\rho_w gH \frac{\partial s}{\partial y} |

</math> | </math> | ||

− | <math> \eta</math> is the non-linear, vertically averaged viscosity. | + | <math> \eta</math> is the non-linear, vertically averaged viscosity, <math>H</math> is the ice shelf thickness and <math>s</math> is the surface elevation. <math> \eta</math> will need to be entered as a '''scalar expression''', and is written |

− | + | ||

− | :<math>\eta = \frac{B}{2}\left[ \left(\frac{\partial u}{\partial x}\right)^2 + \left(\frac{\partial v}{\partial y}\right)^2 + \frac{1}{4} \left(\frac{\partial u}{\partial | + | :<math>\eta = \frac{B}{2}\left[ \left(\frac{\partial u}{\partial x}\right)^2 + \left(\frac{\partial v}{\partial y}\right)^2 + \frac{1}{4} \left(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}\right)^2 + \frac{\partial u}{\partial x}\frac{\partial v}{\partial y}\right]^{-(n-1)/2n}</math> |

+ | |||

+ | ===Boundary conditions=== | ||

+ | There are two flavors of boundary conditions applied in the typical ice sheet model. | ||

+ | ====Kinematic==== | ||

+ | First, the Dirichlet, or ''kinematic'' boundary condition specify the velocity | ||

+ | :<math> \mathbf{u} = \mathbf{u_b},~\forall \partial \Omega_k \in \partial \Omega.</math> | ||

+ | This is the boundary condition where the ice moving across the grounding line. In order to model an ice shelf, one determines (or estimates), and specifies that velocity. | ||

+ | ====Dynamic==== | ||

+ | This is the Neumann, or ''dynamic'' boundary condition that is applied along the ice front, | ||

+ | |||

+ | :<math>- 2 \eta H \left(2\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)\mathbf{n_x} | ||

+ | |||

+ | -\eta H\left( | ||

+ | \frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\mathbf{n_y} | ||

+ | = | ||

+ | -\rho_w g H^2\left(1-\frac{\rho_i}{\rho_w}\right) \mathbf{n_x} </math> | ||

+ | |||

+ | and | ||

+ | |||

+ | <math>- \eta H \left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\mathbf{n_x} | ||

+ | |||

+ | -2\eta H\left( | ||

+ | \frac{\partial u}{\partial x}+2\frac{\partial v}{\partial y}\right)\mathbf{n_y} | ||

+ | = | ||

+ | -\rho_w g H^2\left(1-\frac{\rho_i}{\rho_w}\right) \mathbf{n_y} </math> | ||

+ | |||

+ | |||

+ | :<math> ~\forall \partial \Omega_d \in \partial \Omega.</math> | ||

+ | |||

+ | Note that in case you will need to know how to tell Comsol the normal vectors. These are '''nx''', and '''ny'''. | ||

===Exercises=== | ===Exercises=== | ||

#As a first exercise in solving these equations, try the experiments described in the EISMINT ice shelf models, but never published [http://homepages.vub.ac.be/~phuybrec/eismint/iceshelf.html]. Get the '''self-descr.pdf''', or the first hyperlink on the page. Let's do experiments 3-4 on page 6 of the document (Note that we will see and work with the solution to these experiments again when we [[Adding a module to Glimmer I|do some exercises with the higher-order dynamics routines in Glimmer/CISM]]). | #As a first exercise in solving these equations, try the experiments described in the EISMINT ice shelf models, but never published [http://homepages.vub.ac.be/~phuybrec/eismint/iceshelf.html]. Get the '''self-descr.pdf''', or the first hyperlink on the page. Let's do experiments 3-4 on page 6 of the document (Note that we will see and work with the solution to these experiments again when we [[Adding a module to Glimmer I|do some exercises with the higher-order dynamics routines in Glimmer/CISM]]). | ||

− | # The above is neat, but ultimately not that useful because | + | # The above is neat, but ultimately not that useful because it just isn't a realistic geometry. Try out [[Media:Ross_shelf.zip|this]] model, that has the geometry and boundary conditions for the Ross Ice shelf, but otherwise the same equations solved in the previous exercise. This model is based on MacAyeal et. al. (1996) <ref name="EISMINT_ROSS"> MacAyeal, D.R., V. Rommelaere, Ph. Huybrechts, C.L. Hulbe, J. Determann, and C. Ritz (1996) [http://homepages.vub.ac.be/~phuybrec/pdf/MacAyeal.Ann.Glac.23.pdf pdf]</ref>. An ice-shelf model test based on the Ross ice shelf. Annals of Glaciology 23, 46-51 See the utility of this? Is the solution dependent on the mesh? Can you do anything with the solver to improve the time required for a solution. |

+ | |||

+ | ==Full Stoke's== | ||

+ | It is possible (even easy!) to represent the full Stoke's equations in Comsol, by starting from the fluid mechanics application mode. [[Media:arolla_tripped.zip|Here]] is a geometry for the Arolla glacier. See if you can work out the eta_nonlinear field, the appropriate boundary conditions, and the appropriate subdomain settings for full Stoke's flow. A useful reference on higher order modeling is | ||

+ | Pattyn, F (2003) <ref Name="pattyn_03"> ''J. Geophys. Res.'', '''108''', B8, 2382, {{doi|10.1029/2002JB002329}}</ref>. | ||

==References== | ==References== | ||

<references/> | <references/> |

## Latest revision as of 18:01, 7 August 2009

## Contents |

## Overview

Now, let's see if COMSOL can be used to solve problems of glaciological relevance. We'll look at shallow ice approximation flow, and shallow shelf approximation flows.

## Isothermal Shallow Ice Approximation

Begin with the often used shallow ice form for ice thickness evolution, which casts evolution as a non-linear diffusion problem

where

with boundary condition on the edge of the computational domain.

### Comsol Modeling

We will use the **PDE, General Form** transient mode to solve this equation. For convenience, make the dependent variable .

#### Geometry

You should not find it difficult to create a unit square. Once it's made, you can double click it to change it's size and do other transformations. Read below to find the appropriate domain.

#### Field equations

This equation mode solves equations of the form

Which is just what we want if we recognize that in our system =0, =1, =, and

Now the problem has been reduced to one of typing. It will make the COMSOL model easier to read if you create a scalar expression for . Then your and are very clear.

#### Boundary conditions

This type of problem requires a Dirchlet boundary condition. Set = 0 on all four sides.

#### Other

You'll also be needing to know how to tell COMSOL to use a derivative. That is **Hx**, **Hy**, and **Ht** for , , and respectively.

### Exercises

- Complete the model, and do the isothermal
*fixed margin*experiment Huybrechts (1996)^{[1]}. You'll find all values of constants there as well. Verify that your model is providing results consistent with those reported in the paper. - Now alter your model (the accumulation field) to do the isothermal
*moving margin*, again verify that it's at least just as wrong as the other models. You're going to have to come up with something to deal with the negative values of thickness that you'll get...

## Shallow shelf approximation

### Field equations

Now, consider the equations describing a flow that is vertically integrated. The equations are

is the non-linear, vertically averaged viscosity, is the ice shelf thickness and is the surface elevation. will need to be entered as a **scalar expression**, and is written

### Boundary conditions

There are two flavors of boundary conditions applied in the typical ice sheet model.

#### Kinematic

First, the Dirichlet, or *kinematic* boundary condition specify the velocity

This is the boundary condition where the ice moving across the grounding line. In order to model an ice shelf, one determines (or estimates), and specifies that velocity.

#### Dynamic

This is the Neumann, or *dynamic* boundary condition that is applied along the ice front,

and

Note that in case you will need to know how to tell Comsol the normal vectors. These are **nx**, and **ny**.

### Exercises

- As a first exercise in solving these equations, try the experiments described in the EISMINT ice shelf models, but never published [1]. Get the
**self-descr.pdf**, or the first hyperlink on the page. Let's do experiments 3-4 on page 6 of the document (Note that we will see and work with the solution to these experiments again when we do some exercises with the higher-order dynamics routines in Glimmer/CISM). - The above is neat, but ultimately not that useful because it just isn't a realistic geometry. Try out this model, that has the geometry and boundary conditions for the Ross Ice shelf, but otherwise the same equations solved in the previous exercise. This model is based on MacAyeal et. al. (1996)
^{[2]}. An ice-shelf model test based on the Ross ice shelf. Annals of Glaciology 23, 46-51 See the utility of this? Is the solution dependent on the mesh? Can you do anything with the solver to improve the time required for a solution.

## Full Stoke's

It is possible (even easy!) to represent the full Stoke's equations in Comsol, by starting from the fluid mechanics application mode. Here is a geometry for the Arolla glacier. See if you can work out the eta_nonlinear field, the appropriate boundary conditions, and the appropriate subdomain settings for full Stoke's flow. A useful reference on higher order modeling is
Pattyn, F (2003) ^{[3]}.