Enthalpy Thermomechanical Coupling of "Test Glacier"

Extension of Elmer in computational glaciology
Post Reply
andrewdnolan
Posts: 16
Joined: 26 Oct 2020, 18:33
Antispam: Yes

Enthalpy Thermomechanical Coupling of "Test Glacier"

Post by andrewdnolan »

Hello,

I've just wrapped up the beginners course over the last two weeks. Thanks again for putting the course on, it was great!

As a jump start on some research objectives I’ve been trying to rewrite the thermomechanically coupled diagnostic run from the testglacier_flowline module with an enthalpy formulation. I’ve been working from the Stokes_diagnostic_thermomech.sif file from the course folder and have changed Solver 4 to solve the “Enthalpy Equation” instead of “Homologous Temperature Equation.” Following the enthalpy solver documentation on the Wiki, the enthalpy test on the github repo, and some slides on the Wiki I’ve:
  • added the necessary enthalpy related constants to the Constants section
  • changed the material section to a “power law” Viscosity Model from “Glen”. Not totally sure about this, but all the examples I’ve come across have done this. Is this necessary because the built in “Glen” Viscosity Model is formulated as a function of temperature and not enthalpy?
  • added enthalpy related parameters to the Material section (i.e. Enthalpy Density, Enthalpy Heat Diffusivity , and Enthalpy Water Diffusivity ).
    Changed the basal boundary condition to be a Enthalpy Heat Flux, all of which are in consitent units (MPa - m - a, at least I think).
  • Changed the surface boundary condition to be a prescribed Enthalpy_h which I’ve assumed is the same prescribed temperature as before now converted to enthalpy. I got the formula from integrating the top condition of Equation 14 from Gilber et al 2014.
I’ve gone through my new .sif file and made sure all the base units are in MPa - m - a as was suggested the course. Unfortunately, all of this is to no avail and my results still disagree with the temperature field from "Stokes_diagnostic_thermomech.sif."

I was hoping someone could provide some guidance as to where I am might be going wrong. I've attached my most up to date .sif file for those that are interested.

Thanks in advance,
Andrew
Attachments
diagnostic_enthalpy.sif
(9.63 KiB) Downloaded 409 times
AdrienG
Posts: 7
Joined: 21 Nov 2018, 13:59
Antispam: Yes

Re: Enthalpy Thermomechanical Coupling of "Test Glacier"

Post by AdrienG »

Hi Andrew,

Everything seems correct to me. However if you solve the stoke solver with the "power law" viscosity like you do, you are not computing the viscosity as a function of temperature anymore, it could lead to different result then.

There is no issue with using the Glen viscosity when solving temperature through the enthalpy solver. The enthalpy solver export a variable "temperature" that can be directly used to calculate the viscosity.

If still some problem, send me the whole test case, I can try to run it and help you further.

Best regards,
Adrien
andrewdnolan
Posts: 16
Joined: 26 Oct 2020, 18:33
Antispam: Yes

Re: Enthalpy Thermomechanical Coupling of "Test Glacier"

Post by andrewdnolan »

Thanks for the timely reply!

Good to know about viscosity model, I've switched back to Glen since that's the file I'm trying to reproduce uses.

Taking a look at the outputted temperature field from the enthalpy run it looks plausible, but considerably different than the output of the run with the temperature solver.

I've attached a minimum working example with all the files need to reproduce my output. Any guidance would be much appreciated.

Greatly appreciate you taking the time to look at this.

Best,
Andrew
Attachments
testglacier_enthalpy.tar.gz
(205.26 KiB) Downloaded 384 times
AdrienG
Posts: 7
Joined: 21 Nov 2018, 13:59
Antispam: Yes

Re: Enthalpy Thermomechanical Coupling of "Test Glacier"

Post by AdrienG »

Hi Andrew,

Sorry for late reply. I think that the issue come from the fact that you do unit conversion in material :

Enthalpy Density = real #910.0 / (1.0e6*yearinsec^2) ! [kg m-3] --> [MPa a2 m-2]
Enthalpy Heat Diffusivity = real #(2.1/2050.0) * (1.0e6*yearinsec^-1) ! [kg m-1 s-1] --> [MPa a]
Enthalpy Water Diffusivity = real #1.045e-4 / (1.0e6*yearinsec) ! [kg m-1 s-1] --> [MPa a]

But not in the constant section :

T_ref_enthalpy = real 200.0 ![K]
L_heat = real 334000.0 ![J kg-1]
! Cp(T) = A*T + B
Enthalpy Heat Capacity A = real 7.253 ![J kg-1 K-2]
Enthalpy Heat Capacity B = real 146.3 ![J kg-1 K-1]
P_triple = real 0.061173 ![MPa] Triple point pressure for water (Gilbert et al. 2014)
P_surf = real 0.1013 ![MPa] Surface atmospheric pressure
beta_clapeyron = real 0.0974 ![K MPa-1] clausus clapeyron relationship (Gilbert et al. 2014)

It results in homogeneity problem. The constant are defined for enthalpy in J kg-1. With your definition done in material, you need to modify L_heat, Enthalpy Heat Capacity A , Enthalpy Heat Capacity B accordingly.

Cheers,
Adrien
Post Reply