[RichardL]

Hi everyone,
Please accept my apologies for this post, but I am not really sure what to search for in this forum to find a solution for my issue. I am currently facing an issue with Ansys 18 Mechanical’s transient thermal module. For my research I am trying to simulate the effects surface temperature variations on the subsurface geothermal regime. However, when QCing my Ansys results I observe significant deviations from a 1D solution for this diffusion problem, even in a very simple case.
Prior to working with my actual model I am testing the transient thermal module in Ansys mechanical on a simple block model. As a geologist I naturally think in larger dimensions, so that my simple model extends over 20 000 x 20 000 x 6 000m. Within these bounds I modelled three layers, which each a height of 2000m. The three layers consist of the same material, with the following parameters:
Density: 3000 kg/m3
Thermal Conductivity: 2.5 W/ (m K)
Specific heat capacity: 1190 J / (kg K)
I meshed the block with patch conforming tetrahedrons and defined a refinement are in the very centre of the model, from which I extract the results along a path. The ‘body sizing’ option is located at the top-centre of the model (10 000, 6 000, 10 000) with a sphere size of 2 500m and an element size of 150m. From the same point I defined a path vertically down to 3000m.
To test my model against a 1D solution, I insulated the side walls (12 faces), provide a constant heat flux on the base and defined a surface temperature on the top face. Firstly, I ran my model as a static model, with a heat flux of 0.05 W/m2 and a top surface temperature of 0 degC. The result looked convincing to me: a constant heat flux of 50mW/m2 through the entire model, a stable thermal gradient of 20degC/km and an evenly distributed temperature.
Next, I built the transient model using the input temperatures from my steady state model. For the transient model I kept the heat flux constant at 50mW/m2 for all time steps, but varied the surface temperature. I went for a simplified model mimicking the surface temperature changes during the last ice age. Please note, I used the duration of the events, not the exact times before present:
t(years) T(degC)
10^5 -4
10^4 3
6E+3 4
1E+3 3
2E+2 3.5
To bring this scheme into Ansys I converted years to seconds and added small time steps (1000 sec) to get a stepped temperature time profile:

0 -4

3153631536000 -4

3153631537000 3

34689915360003

3468991537000 4

3658207536000 4

3658207537000 3

3689743536000 3

3689743537000 3.5

3696050736000 3.5

A solution for the model is found after 5 Minutes and that’s where the trouble starts. For my 1D reference model I calculated HF according the equation below using the same input parameters as used above:
T(t, z) = T0 + dT*erfc(z/(2*sqrt(k*t)))
Where, T0 denotes the temperature of the previous time step; dT the temperature difference between the time steps, z depth below surface, k thermal diffusivity (m2/s) and t time in sec. To calculate T0 for the first step I used the thermal gradient from my static model, 20degC/km. Furthermore, I used a diffusivity of 0.7E-6 m2/s, which corresponds to the material properties defined in the Ansys model. Once the new temperatures were calculated I calculated the thermal gradient for each modelling step, i.e. every 10m, for all time steps. Finally, I calculated the heat flow by multiplying the thermal gradient with the thermal conductivity, i.e.2.5 W/(m K). In the steady state case this gives a ‘base heat flow’ of 50mW/m2.
I compared the data of my 1D model against the data calculated along the previously defined path in Ansys for each time step. For the first three time steps, i.e. steady state, 10^5 and 10^4 I get a perfect match between both models.
And now something really strange happens, for the time steps of 6000, 1000 and 200 years the Ansys solution moves to lower heat flow values with depth. The intersection between the 50mW/m2 is shifted 300 deeper, i.e. 1600m. Moreover, the heat flow values at the surface from are lower than the ones calculated with the 1D model.
To me it seems that for some reasons the material properties changed as even with a wrong timing of the temperature changes, the subsurface would more or less equilibrate at the same depth levels so that heat flow curves would meet much earlier. Is anybody there with an idea what happens here? I personally find only a few option in Ansys to play around with… Are there any suggestions on this side?
Many thanks in advance!
Cheers,
Richard