Solving Creep Problems with FLAC3D

The major difference between creep and other constitutive models in FLAC3D is the concept of problem time in the simulation. For creep runs, the problem time and timestep represent real time, while for static analyses in the other constitutive models, the timestep is an artificial quantity, used only as a means of stepping to a steady-state condition.

Creep Timestep

For time-dependent phenomena such as creep, FLAC3D allows the user to define a timestep. The default for this timestep is zero, in which case the program treats the material as linearly elastic (viscoelastic models) or elasto-plastic (viscoplastic models), as appropriate. (The command model creep active off can also be used to stop the creep calculations.) This can be used to attain equilibrium before starting a creep simulation. The constitutive laws for creep make use of the timestep in their equations, so timestep may affect the response.

Although the user may set the timestep, it is not arbitrary. For a system to always be in quasi-static mechanical equilibrium (as in a creep simulation), the time-dependent stress changes produced by the constitutive law must not be large compared to the strain-dependent stress changes. Otherwise, out-of-balance forces will be large, and inertial effects (which are theoretically absent) may affect the solution.

The creep processes are governed by the deviatoric stress state. An estimate for the maximum creep timestep for numerical accuracy can be expressed as the ratio of the material viscosity to the shear modulus,

(1)\[\Delta t^{cr}_{max} = {\eta \over G}\]

For a power-type creep model, the viscosity may be estimated as the ratio of the stress magnitude, \(\bar\sigma\), to the creep rate, \(\dot\epsilon_{cr}\). Using this equation as in the power model, the maximum creep timestep is

(2)\[\Delta t^{cr}_{max} = {\bar\sigma^{1-n} \over {A G}}\]

For a WIPP-type model, the viscosity may be estimated as the ratio of \(\bar\sigma\) to the secondary creep rate, \(\dot\epsilon_s\), and using this equation as in the WIPP model, the maximum creep timestep is

(3)\[\Delta t^{cr}_{max} = {e^{Q/RT} \over {G \ D \ \bar\sigma^{n-1}}}\]

For a Burgers-type model, Equation (1) must be interpreted as

(4)\[\Delta t_{\max }^{cr}=\min \left( {{\eta^K} \over {G^K}}, {{\eta^M} \over {G^M}} \right)\]

where the superscripts \(.^K\) and \(.^M\) refer to Kelvin and Maxwell properties, respectively.

The timestep limitation for creep compaction (e.g., WIPP-salt model) involves the volumetric response of the system, and is estimated as the ratio of viscosity to bulk modulus. This viscosity may be expressed as the ratio of \(\bar\sigma\) to the volumetric creep compaction rate, \(\dot\epsilon^c_v\). Using this equation as in the WIPP-salt model, the maximum creep timestep for creep compaction is

(5)\[\Delta t^{cr}_{max} = {\vert\sigma\vert \rho \over {K B_0 \left[e^{B_1 \vert\sigma\vert } - 1 \right] e^{B_2 \rho}}}\]

It is recommended that a creep analysis with FLAC3D begin with an initial creep timestep approximately two to three orders of magnitude smaller than \(\Delta t^{cr}_{max}\), as calculated from the appropriate formula above. By invoking model creep timestep automatic on, use can then be made of the automatic timestep adjustment. As a rule, the maximum value for the timestep (model creep timestep maximum) should not exceed the value derived for \(\Delta t^{cr}_{max}\). See Verifications and Example Problems for example applications.

The stress magnitude, \(\bar\sigma\), used in the calculation for \(\Delta t^{cr}_{max}\), can be determined from the initial stress state before the creep process begins. \(\bar\sigma\), also known as the von Mises stress invariant, can be calculated from the FISH function. The maximum \(\bar\sigma\) in the FLAC3D model should be used to calculate \(\Delta t^{cr}_{max}\).

;  calculate maximum von Mises stress in pipe
fish define mises
    local max = 0.0
    loop foreach local elem struct.list
        if struct.pos.y(elem) < 2.0 then
            local j2  = tensor.j2(struct.shell.stress(elem,0))
            local sr_vmstr  = math.sqrt(3.0*j2)
            max = math.max(max,sr_vmstr)
        endif
    endloop
    mises = max
end

Automatic Adjustment of the Creep Timestep

The timestep may be set by the user to a constant value, or controlled by FLAC3D to change automatically. If the timestep is changed automatically, it can be decreased whenever the maximum unbalanced force exceeds some threshold, and increased whenever it goes below some other level. The threshold is defined as the ratio of the maximum unbalanced force to the average gridpoint force.

Typical out-of-balance force criteria for the problem being solved can be determined by observing the out-of-balance force that occurs near equilibrium in the initial stage of the problem when only elastic effects are present. In many cases, a good performance can be obtained by using a gradual increase or decrease of timestep (e.g., with the default ratios lower-multiplier = 1.01 and upper-multiplier = 0.90).

In some cases, it may be preferable to avoid a continuous adjustment of the timestep, which may create “noise.” For this purpose, after a timestep change has occurred, there is a user-defined “latency period” (e.g., 100 steps) during which no further adjustments are made, allowing the system to settle. Normally, the timestep will start at a small value, to accommodate transients such as excavation, and then increase as the simulation proceeds. If a new transient is introduced, it may be desirable to reduce the timestep manually, and then let it increase again automatically.

The model creep timestep command is used to set the timestep and the parameters required to allow timestep to change automatically. The keywords are listed in the command.

Temperature Dependency

For the WIPP model, the creep rate is temperature-dependent. Temperatures may be supplied to the WIPP-type models in one of two ways: they may be specified as a property of the WIPP-type models; or they may be calculated using the thermal option of the code. With the first approach, temperatures are assigned with the zone property command and do not change during the calculation. In the second approach, the model configure thermal command must first be given and temperatures can be initialized with the zone gridpoint initialize temperature command. In the latter case, temperatures can change during cycling. A temperature gradient may be specified with either approach.