Dear all,
I am working on a coupled thermal-creep model. The coupling would be one way from thermal to creep through the thermal dependency of the creep parametrisation (I am using the power constitutive model). This basically through the temperature term in a Norton-Hoff creep law:
\dot \epsilon = A^{-Q\RT} \sigma^{n}
In my experience with other software, this is a fairly straightforward coupling, but I do not have much experience with code like FLAC that include the effect of thermal expansion in the coupling.
Therefore, before doing what I described above I am first testing the the simple thermo-mechanical coupling with just thermal expansion. as follow:
- Take a simulation that run with creep only.
- activate thermal and set the thermal parameters and boundary conditions
- I do not implement the coupling through the creep law
- set the thermal expansion term to a very low value (1e-12 K-1)
- solve to a set time with : "model solve time [time] "
To my understanding with this configuration I should have the same results as if I would run creep and thermal separately. Instead I get surprising results:
- The solve timestep is much lower - by one or two order of magnitude - than the timestep when solving with creep or thermal only.
- I observed oscillations in the velocity field when I track velocity history at key zones in the model that then rapidly dampen out.
Fig 1: velocity history for thermal expansion = 1e-12 K-1
Fig 2: velocity history for creep only solve
Fig 3: velocity history with fixed timestep at a larger value ( ~ 10 time the timestep in fig 1)
Fig 4: velocity history for thermal expansion = 1e-6 K-1
Fig 4: velocity history for thermal expansion = 1e-6 K-1, with relaxation phase before.
This is fairly surprising. As I don’t understand were these oscillations are coming from. Am I using the coupling incorrectly?
-
If I drop further the thermal expansion coefficient to 1e-18, I get similar timestep and oscillations as for 10-e12.
-
If in increase the thermal expansion to 1e-6 - which is around the expected value for rock salt - I get very large oscillations and the timestep get 1 order of magnitude smaller. This oscillation dampen fairly quickly. Yet, it feel a bit much, as the thermal expansion effect on stress should be relatively small. If I do a back of the envelop calculation, equation 36 of the this page of the documentation (Numerical Formulation — Itasca Software 9.0 documentation), I should expect a maximum possible stress from thermal expansion of ~0.6 MPa, this is about a tenth of the deviatoric stress I have from creep.
Then, I try I activated the coupling between creep and thermal in the Norton-Hoff creep law, by having a fish function update the parameters of the constitutive model at each timestep. This seem to work fine, I does not create any additional oscillations or slower convergence, which is what I expected, as the thermal effect cool down the rock is relatively slow and reduces the creep rate.
Then, I try to relax the system first, by running the thermo-elastic simulation (I deactivate creep) for a short time before running the thermal-creep model. It does not seem to work that well (Figure 5) and Frankly I am not sure I am configuring the solve correctly, I have the feeling there is a better way to do this
model thermal timestep maximum [3600*24]
zone thermal implicit offmodel creep active off
model solve-thermal-coupled time [3600*24] synchronizemodel thermal active off
model solve ratio 1e-6model creep on
model thermal onmodel solve time [3600 24 365]
Finally, I tried to just fix the timestep for both creep and thermal to a value that I know if stable for creep and thermal only for this problem using the following commands:
model creep timestep fix 15000
model thermal timestep fix 15000
This seem to work and it even dampen the oscillation. So I am I have the feeling that these oscillations that I observe at smaller timestep are a transient creep relaxation mode that do not impact the quality of my solution. If that is the case I am still very surprised it happen given the very small value I uses for thermal expansion.
Would someone have experience with thermal-creep coupling in FLAC3D and could give me a pointer on what to could be the effect or if I made a mistake could point be to were to look into?
Are my commands for the thermo-elastic relaxation step correct and doing what I thing it is or am I doing something strange?
Thank you I advance!
As an additional info, I the numerical stability of the creep only simulation is well test so it is unlikely to some numerical problem that is already there in the creep only solve.
Also the oscillation I show are worst offender, if I trace history of stress at the same point I barely see the oscillations.




