Fully Coupled Thermal-Stress Analysis

Fully coupled thermal-stress analysis is needed when the stress analysis is dependent on the temperature distribution and the temperature distribution depends on the stress solution. For example, metalworking problems may include significant heating due to inelastic deformation of the material that, in turn, changes the material properties. In addition, contact conditions exist in some problems where the heat conducted between surfaces may depend strongly on the separation of the surfaces or the pressure transmitted across the surfaces (see Thermal Contact Properties). For such cases the thermal and mechanical solutions must be obtained simultaneously rather than sequentially. Coupled temperature-displacement elements are provided for this purpose in both Abaqus/Standard and Abaqus/Explicit; however, each program uses different algorithms to solve coupled thermal-stress problems.

Fully Coupled Thermal-Stress Analysis in Abaqus/Standard

In Abaqus/Standard the temperatures are integrated using a backward-difference scheme, and the nonlinear coupled system is solved using Newton's method. Abaqus/Standard offers an exact as well as an approximate implementation of Newton's method for fully coupled temperature-displacement analysis.

Exact Implementation

An exact implementation of Newton's method involves a nonsymmetric Jacobian matrix as is illustrated in the following matrix representation of the coupled equations:

[\begin{matrix} K_{u u} & K_{u θ} \\ K_{θ u} & K_{θ θ} \end{matrix}] {\begin{matrix} Δ u \\ Δ θ \end{matrix}} = {\begin{matrix} R_{u} \\ R_{θ} \end{matrix}},

where $Δ u$ and $Δ θ$ are the respective corrections to the incremental displacement and temperature, $K_{i j}$ are submatrices of the fully coupled Jacobian matrix, and $R_{u}$ and $R_{θ}$ are the mechanical and thermal residual vectors, respectively.

Solving this system of equations requires the use of the unsymmetric matrix storage and solution scheme. In addition, the mechanical and thermal equations must be solved simultaneously. The method provides quadratic convergence when the solution estimate is within the radius of convergence of the algorithm. The exact implementation is used by default.

Approximate Implementation

Some problems require a fully coupled analysis in the sense that the mechanical and thermal solutions evolve simultaneously, but with a weak coupling between the two solutions. In other words, the components in the off-diagonal submatrices $K_{u θ}, K_{θ u}$ are small compared to the components in the diagonal submatrices $K_{u u}, K_{θ θ}$ . An example of such a situation is the disc brake problem (Thermal-stress analysis of a disc brake). For these problems a less costly solution may be obtained by setting the off-diagonal submatrices to zero so that we obtain an approximate set of equations:

[\begin{matrix} K_{u u} & 0 \\ 0 & K_{θ θ} \end{matrix}] {\begin{matrix} Δ u \\ Δ θ \end{matrix}} = {\begin{matrix} R_{u} \\ R_{θ} \end{matrix}} .

As a result of this approximation the thermal and mechanical equations can be solved separately, with fewer equations to consider in each subproblem. The savings due to this approximation, measured as solver time per iteration, will be of the order of a factor of two, with similar significant savings in solver storage of the factored stiffness matrix. Further, in many situations the subproblems may be fully symmetric or approximated as symmetric, so that the less costly symmetric storage and solution scheme can be used. The solver time savings for a symmetric solution is an additional factor of two. Unless you explicitly choose the unsymmetric matrix storage and solution scheme, selection of the scheme depends on other details of the problem (see Defining an Analysis).

This modified form of Newton's method does not affect solution accuracy since the fully coupled effect is considered through the residual vector $R_{j}$ at each increment in time. However, the rate of convergence is no longer quadratic and depends strongly on the magnitude of the coupling effect, so more iterations are generally needed to achieve equilibrium than with the exact implementation of Newton's method. When the coupling is significant, the convergence rate becomes very slow and may prohibit obtaining a solution. In such cases the exact implementation of Newton's method is required. In cases where it is possible to use this approximation, the convergence in an increment depends strongly on the quality of the first guess to the incremental solution, which you can control by selecting the extrapolation method used for the step (see Defining an Analysis).

Steady-State Analysis

A steady-state coupled temperature-displacement analysis can be performed in Abaqus/Standard. In steady-state cases you should assign an arbitrary “time” scale to the step: you specify a “time” period and “time” incrementation parameters. This time scale is convenient for changing loads and boundary conditions through the step and for obtaining solutions to highly nonlinear (but steady-state) cases; however, for the latter purpose, transient analysis often provides a natural way of coping with the nonlinearity.

Frictional slip heat generation is normally neglected in for the steady-state case. However, it can still be accounted for if user subroutine FRIC provides the incremental frictional dissipation through the variable SFD. If frictional heat generation is present, the heat flux into the two contact surfaces depends on the slip rate of the surfaces. The “time” scale in this case cannot be described as arbitrary, and a transient analysis should be performed.

Transient Analysis

Alternatively, you can perform a transient coupled temperature-displacement analysis. You can control the time incrementation in a transient analysis directly, or Abaqus/Standard can control it automatically. Automatic time incrementation is generally preferred.

Automatic Incrementation Controlled by a Maximum Allowable Temperature Change

The time increments can be selected automatically based on a user-prescribed maximum allowable nodal temperature change in an increment, $Δ θ_{m a x}$ . Abaqus/Standard restricts the time increments to ensure that this value is not exceeded at any node (except nodes with boundary conditions) during any increment of the analysis (see Time Integration Accuracy in Transient Problems).

Fixed Incrementation

If you do not specify $Δ θ_{m a x}$ , fixed time increments equal to the user-specified initial time increment, $Δ t_{0}$ , is used throughout the analysis, except when the explicit creep integration scheme is used. In this case Abaqus/Standard might decrease the time increment if the stability limit is exceeded.

Spurious Oscillations due to Small Time Increments

In transient analysis with second-order elements there is a relationship between the minimum usable time increment and the element size. A simple guideline is

Δ t > \frac{ρ c}{6 k} Δ ℓ^{2},

where $Δ t$ is the time increment, $ρ$ is the density, c is the specific heat, k is the thermal conductivity, and $Δ ℓ$ is a typical element dimension (such as the length of a side of an element). If time increments smaller than this value are used in a mesh of second-order elements, spurious oscillations can appear in the solution, in particular in the vicinity of boundaries with rapid temperature changes. These oscillations are nonphysical and may cause problems if temperature-dependent material properties are present. In transient analyses using first-order elements the heat capacity terms are lumped, which eliminates such oscillations but can lead to locally inaccurate solutions for small time increments. If smaller time increments are required, a finer mesh should be used in regions where the temperature changes rapidly.

There is no upper limit on the time increment size (the integration procedure is unconditionally stable) unless nonlinearities cause convergence problems.

Automatic Incrementation Controlled by the Creep Response

The accuracy of the integration of time-dependent (creep) material behavior is governed by the user-specified accuracy tolerance parameter, $t o l e r a n c e \geq ({\dot{\bar{ε}}}^{c r}_{t + Δ t} - {\dot{\bar{ε}}}^{c r}_{t}) Δ t$ . This parameter is used to prescribe the maximum strain rate change allowed at any point during an increment, as described in Rate-Dependent Plasticity: Creep and Swelling. The accuracy tolerance parameter can be specified together with the maximum allowable nodal temperature change in an increment, $Δ θ_{m a x}$ (described above); however, specifying the accuracy tolerance parameter activates automatic incrementation even if $Δ θ_{m a x}$ is not specified.

Selecting Explicit Creep Integration

Nonlinear creep problems (Rate-Dependent Plasticity: Creep and Swelling) that exhibit no other nonlinearities can be solved efficiently by forward-difference integration of the inelastic strains if the inelastic strain increments are smaller than the elastic strains. This explicit method is efficient computationally because, unlike implicit methods, iteration is not required if no other nonlinearities are present. Although this method is only conditionally stable, the numerical stability limit of the explicit operator is in many cases sufficiently large to allow the solution to be developed in a reasonable number of time increments.

For most coupled thermal-stress analyses, however, the unconditional stability of the backward difference operator (implicit method) is desirable. In such cases the implicit integration scheme may be invoked automatically by Abaqus/Standard.

Explicit integration can be less expensive computationally and simplifies implementation of user-defined creep laws in user subroutine CREEP; you can restrict Abaqus/Standard to using this method for creep problems (with or without geometric nonlinearity included). See Rate-Dependent Plasticity: Creep and Swelling for further details.

Excluding Creep and Viscoelastic Response

You can specify that no creep or viscoelastic response occurs during a step even if creep or viscoelastic material properties have been defined.

Unstable Problems

Some types of analyses may develop local instabilities, such as surface wrinkling, material instability, or local buckling. In such cases it may not be possible to obtain a quasi-static solution, even with the aid of automatic incrementation. Abaqus/Standard offers a method of stabilizing this class of problems by applying damping throughout the model in such a way that the viscous forces introduced are sufficiently large to prevent instantaneous buckling or collapse but small enough not to affect the behavior significantly while the problem is stable. The available automatic stabilization schemes are described in detail in Automatic Stabilization of Unstable Problems.

Units

In coupled problems where two different fields are active, take care when choosing the units of the problem. If the choice of units is such that the terms generated by the equations for each field are different by many orders of magnitude, the precision on some computers may be insufficient to resolve the numerical ill-conditioning of the coupled equations. Therefore, choose units that avoid ill-conditioned matrices. For example, consider using units of megapascal (MPa) instead of pascal (Pa) for the stress equilibrium equations to reduce the disparity between the magnitudes of the stress equilibrium equations and the heat flux continuity equations.

Fully Coupled Thermal-Stress Analysis in Abaqus/Explicit

In Abaqus/Explicit the heat transfer equations are integrated using the explicit forward-difference time integration rule

θ_{(i + 1)}^{N} = θ_{(i)}^{N} + Δ t_{(i + 1)} {\dot{θ}}_{(i)}^{N},

where $θ^{N}$ is the temperature at node N and the subscript i refers to the increment number in an explicit dynamic step. The forward-difference integration is explicit in the sense that no equations need to be solved when a lumped capacitance matrix is used. The current temperatures are obtained using known values of ${\dot{θ}}_{(i)}^{N}$ from the previous increment. The values of ${\dot{θ}}_{(i)}^{N}$ are computed at the beginning of the increment by

{\dot{θ}}_{(i)}^{N} = {(C^{N J})}^{- 1} (P_{(i)}^{J} - F_{(i)}^{J}),

where $C^{N J}$ is the lumped capacitance matrix, $P^{J}$ is the applied nodal source vector, and $F^{J}$ is the internal flux vector.

The mechanical solution response is obtained using the explicit central-difference integration rule with a lumped mass matrix as described in Explicit Dynamic Analysis. Since both the forward-difference and central-difference integrations are explicit, the heat transfer and mechanical solutions are obtained simultaneously by an explicit coupling. Therefore, no iterations or tangent stiffness matrices are required.

Explicit integration can be less expensive computationally and simplifies the treatment of contact. For a comparison of explicit and implicit direct-integration procedures, see About Dynamic Analysis Procedures.

Stability

The explicit procedure integrates through time by using many small time increments. The central-difference and forward-difference operators are conditionally stable. The stability limit for both operators (with no damping in the mechanical solution response) is obtained by choosing

Δ t \leq m i n (\frac{2}{ω_{m a x}}, \frac{2}{λ_{m a x}}),

where $ω_{m a x}$ is the highest frequency in the system of equations of the mechanical solution response and $λ_{m a x}$ is the largest eigenvalue in the system of equations of the thermal solution response.

Estimating the Time Increment Size

An approximation to the stability limit for the forward-difference operator in the thermal solution response is given by

Δ t \approx \frac{L_{m i n}^{2}}{2 α},

where $L_{m i n}$ is the smallest element dimension in the mesh and $α = \frac{k}{ρ c}$ is the thermal diffusivity of the material. The parameters k, $ρ$ , and c represent the material's thermal conductivity, density, and specific heat, respectively.

In most applications of explicit analysis the mechanical response governs the stability limit. The thermal response may govern the stability limit when material parameter values are nonphysical or a very large amount of mass scaling is used. The calculation of the time increment size for the mechanical solution response is discussed in Explicit Dynamic Analysis.

Stable Time Increment Report

Abaqus/Explicit writes a report to the status (.sta) file during the data check phase of the analysis that contains an estimate of the minimum stable time increment and a listing of the elements with the smallest stable time increments and their values. The initial minimum stable time increment accounts for the stability requirements of both the thermal and mechanical solution responses. The initial stable time increments listed do not include damping (bulk viscosity), mass scaling, or penalty contact effects in the mechanical solution response.

This listing is provided because often a few elements have much smaller stability limits than the rest of the elements in the mesh. The stable time increment can be increased by modifying the mesh to increase the size of the controlling element or by using appropriate mass scaling.

Time Incrementation

The time increment used in an analysis must be smaller than the stability limits of the central- and forward-difference operators. Failure to use such a time increment results in an unstable solution. When the solution becomes unstable, the time history response of solution variables, such as displacements, usually oscillates with increasing amplitudes. The total energy balance also changes significantly.

Abaqus/Explicit has two strategies for time incrementation control: fully automatic time incrementation (where the code accounts for changes in the stability limit) and fixed time incrementation.

Scaling the Time Increment

To reduce the chance of a solution going unstable, the stable time increment computed by Abaqus/Explicit can be adjusted by a constant scaling factor. This factor can be used to scale the default global time estimate, the element-by-element estimate, or the fixed time increment based on the initial element-by-element estimate; it cannot be used to scale a fixed time increment that you specified directly.

Automatic Time Incrementation

The default time incrementation scheme in Abaqus/Explicit is fully automatic and requires no user intervention. Two types of estimates are used to determine the stability limit: element-by-element for both the thermal and mechanical solution responses and global for the mechanical solution response. An analysis always starts by using the element-by-element estimation method and may switch to the global estimation method under certain circumstances, as explained in Explicit Dynamic Analysis.

In an analysis Abaqus/Explicit initially uses a stability limit based on the thermal and mechanical solution responses in the whole model. This element-by-element estimate is determined using the smallest time increment size due to the thermal and mechanical solution responses in each element.

The element-by-element estimate is conservative; it gives a smaller stable time increment than the true stability limit, which is based on the maximum frequency of the entire model. In general, constraints such as boundary conditions and kinematic contact have the effect of compressing the eigenvalue spectrum, and the element-by-element estimates do not take this into account (see Explicit Dynamic Analysis)

The stable time increment size due to the mechanical solution response is determined by the global estimator as the step proceeds unless the element-by-element estimator is chosen, fixed time incrementation is specified, or one of the conditions explained in Explicit Dynamic Analysis prevents the use of global estimation. The stable time increment size due to the thermal solution response is always determined by using an element-by-element estimation method. The switch to the global estimation method in mechanical solution response occurs once the algorithm determines that the accuracy of the global estimation method is acceptable. For details, see Explicit Dynamic Analysis

For three-dimensional continuum elements and elements with plane stress formulations (shell, membrane, and two-dimensional plane stress elements) an improved estimate of the element characteristic length is used by default. This improved method usually results in a larger element stable time increment than a more traditional method. For analyses using variable mass scaling, the total mass added to achieve a given stable time increment is less with the improved estimate.

Fixed Time Incrementation

A fixed time incrementation scheme is also available in Abaqus/Explicit. The fixed time increment size is determined either by the initial element-by-element stability estimate for the step or by a user-specified time increment.

Fixed time incrementation may be useful when a more accurate representation of the higher mode response of a problem is required. In this case a time increment size smaller than the element-by-element estimates may be used. The element-by-element estimate can be obtained by running a data check analysis (see Abaqus/Standard and Abaqus/Explicit Execution).

When fixed time incrementation is used, Abaqus/Explicit does not check that the computed response is stable during the step. You should ensure that a valid response has been obtained by carefully checking the energy history and other response variables.

If you choose to use time increments the size of the initial element-by-element stability limit throughout a step, the dilatational wave speed and the thermal diffusivity in each element at the beginning of the step are used to compute the fixed time increment size. To reduce the chance of a solution going unstable, the initial stable time increment that Abaqus/Explicit computes can be adjusted by a constant scaling factor, as described above in Scaling the Time Increment. Alternatively, you can specify a time increment size directly.

Reducing the Computational Cost by Using Selective Subcycling

The selective subcycling method can be used in a coupled thermal-stress analysis exactly as in a pure mechanical analysis, as described in Explicit Dynamic Analysis and Selective Subcycling.

Monitoring Output Variables for Extreme Values

The extreme values defined as the element and nodal variables in a coupled thermal-stress analysis can be monitored exactly as described in Explicit Dynamic Analysis for a pure mechanical analysis.

Initial Conditions

By default, the initial temperature of all nodes is zero. You can specify nonzero initial temperatures. Initial stresses, field variables, etc. can also be defined; Initial Conditions describes all of the initial conditions that are available for a fully coupled thermal-stress analysis.

Boundary Conditions

Boundary conditions can be used to prescribe both temperatures (degree of freedom 11) and displacements/rotations (degrees of freedom 1–6) at nodes in fully coupled thermal-stress analysis (see Boundary Conditions). Shell elements in Abaqus/Standard have additional temperature degrees of freedom 12, 13, etc. through the thickness (see Conventions).

Boundary conditions can be specified as functions of time by referring to amplitude curves (Amplitude Curves).

Boundary conditions applied during a dynamic coupled temperature-displacement response step should use appropriate amplitude references (Amplitude Curves). If boundary conditions are specified for the step without amplitude references, they are applied instantaneously at the beginning of the step. Since Abaqus/Explicit does not admit jumps in displacement, the value of a nonzero displacement boundary condition that is specified without an amplitude reference is ignored, and a zero velocity boundary condition is enforced.

Loads

The following types of thermal loads can be prescribed in a fully coupled thermal-stress analysis, as described in Thermal Loads:

Concentrated heat fluxes.
Body fluxes and distributed surface fluxes.
Node-based film and radiation conditions.
Average-temperature radiation conditions.
Element and surface-based film and radiation conditions.

The following types of mechanical loads can be prescribed:

Concentrated nodal forces can be applied to the displacement degrees of freedom (1–6); see Concentrated Loads.
Distributed pressure forces or body forces can be applied; see Distributed Loads. The distributed load types available with particular elements are described in Abaqus Elements Guide.

Predefined Fields

Predefined temperature fields are not allowed in a fully coupled thermal-stress analysis. Boundary conditions should be used instead to prescribe temperature degree of freedom 11 (and 12, 13, etc. in Abaqus/Standard shell elements), as described earlier.

Other predefined field variables can be specified in a fully coupled thermal-stress analysis. These values affect only field-variable-dependent material properties, if any. See Predefined Fields.

Material Options

The materials in a fully coupled thermal-stress analysis must have both thermal properties, such as conductivity, and mechanical properties, such as elasticity, defined. See Abaqus Materials Guide for details on the material models available in Abaqus.

Thermal strain arises if thermal expansion (Thermal Expansion) is included in the material property definition.

In Abaqus/Standard a fully coupled temperature-displacement analysis can be used to analyze static creep and swelling problems, which generally occur over fairly long time periods (Rate-Dependent Plasticity: Creep and Swelling); viscoelastic materials (Time Domain Viscoelasticity); or viscoplastic materials (Rate-Dependent Yield).

Rate-Dependent Yield and Friction in Abaqus/Standard

In Abaqus/Standard you can control whether to consider or ignore the strain rate–dependence of the yield stress and the slip rate–dependence of the friction coefficient within the step.

Internal Heat Generation

In Abaqus/Standard analyses, you can define volumetric heat generation within a material in user subroutine HETVAL or user subroutine UMATHT. You can use these two user subroutines in the same analysis.

In Abaqus/Explicit analyses, you can define volumetric heat generation within a material in user subroutine VHETVAL or user subroutine VUMATHT. You can use these two user subroutines in the same analysis.

Defining Internal Heat Generation in User Subroutine HETVAL or VHETVAL

If you define internal heat generation in user subroutine HETVAL or VHETVAL, you must include heat generation in the material definition with the other thermal property definitions.

Heat generation might be associated with (relatively low) energy phase changes occurring during the solution. Such heat generation usually depends on state variables (such as the fraction transformed), which themselves evolve with the solution and are stored as solution-dependent state variables (see About User Subroutines and Utilities). The heat generation is computed in user subroutine HETVAL or VHETVAL, where any associated state variables can also be updated. The subroutine is called at all material calculation points for which the material definition includes heat generation.

Defining Internal Heat Generation in User Subroutine UMATHT or VUMATHT

If user subroutine UMATHT or VUMATHT is used to define internal heat generation, the constitutive thermal behavior must also be defined within the subroutine.

Inelastic Energy Dissipation as a Heat Source

You can specify an inelastic heat fraction in a fully coupled thermal-stress analysis to provide for inelastic energy dissipation as a heat source. The heat flux per unit volume, $r^{p l}$ , that is added into the thermal energy balance is computed using the equation

r^{p l} = η σ : {\dot{ε}}^{p l},

or, in the case when the nonlinear isotropic/kinematic hardening model is used, from the following equation:

r^{p l} = η (σ - α) : {\dot{ε}}^{p l},

where $η$ is a user-defined factor (assumed constant), $σ$ is the stress, $α$ is the backstress, and ${\dot{ε}}^{p l}$ is the rate of plastic straining.

Inelastic heat fractions are typically used in the simulation of high-speed manufacturing processes involving large amounts of inelastic strain, where the heating of the material caused by its deformation significantly influences temperature-dependent material properties. The generated heat is treated as a volumetric heat flux source term in the heat balance equation.

An inelastic heat fraction can be specified for materials with plastic behavior that use the Mises or Hill yield surface (Inelastic Behavior). It cannot be used with the combined isotropic/kinematic hardening model. The inelastic heat fraction can be specified for user-defined material behavior in Abaqus/Explicit and is multiplied by the inelastic energy dissipation coded in the user subroutine to obtain the heat flux. In Abaqus/Standard the inelastic heat fraction cannot be used with user-defined material behavior; in this case the heat flux that must be added to the thermal energy balance is computed directly in the user subroutine.

An inelastic heat fraction can also be specified for material definitions that include time-domain linear viscoelasticity (Time Domain Viscoelasticity) and time-domain nonlinear viscoelasticity defined within the parallel rheological framework (Parallel Rheological Framework), except in Abaqus/Explicit for large-strain linear viscoelasticity. For large-strain linear viscoelasticity in Abaqus/Standard (Time Domain Viscoelasticity), the energy dissipation is computed only approximately. Hence, the fraction of the dissipated energy converted into heat can be computed only approximately.

The default value of the inelastic heat fraction is 0.9. If you do not include the inelastic heat fraction behavior in the material definition, the heat generated by inelastic deformation is not included in the analysis.

Elements

Coupled temperature-displacement elements that have both displacements and temperatures as nodal variables are available in both Abaqus/Standard and Abaqus/Explicit (see Choosing the Appropriate Element for an Analysis Type). In Abaqus/Standard simultaneous temperature/displacement solution requires the use of such elements; pure displacement elements can be used in part of the model in the fully coupled thermal-stress procedure, but pure heat transfer elements cannot be used. In Abaqus/Explicit any of the available elements can be used in the fully coupled thermal-stress procedure; however, the thermal solution is obtained only at nodes where the temperature degree of freedom has been activated (that is, at nodes attached to coupled temperature-displacement elements).

The first-order coupled temperature-displacement elements in Abaqus use a constant temperature over the element to calculate thermal expansion. The second-order coupled temperature-displacement elements in Abaqus/Standard use a lower-order interpolation for temperature than for displacement (parabolic variation of displacements and linear variation of temperature) to obtain a compatible variation of thermal and mechanical strain.

Output

See Abaqus/Standard Output Variable Identifiers and Abaqus/Explicit Output Variable Identifiers for a complete list of output variables. The types of output available are described in About Output.

Input File Template

HEADING
…
** Specify the coupled temperature-displacement element type
ELEMENT, TYPE=CPS4T
…
**
STEP
COUPLED TEMPERATURE-DISPLACEMENT or
DYNAMIC TEMPERATURE-DISPLACEMENT, EXPLICIT
Data line to define incrementation
BOUNDARY
Data lines to define nonzero boundary conditions on displacement or
temperature degrees of freedom
CFLUX and/or CFILM and/or 
CRADIATE and/or DFLUX and/or 
DSFLUX and/or FILM and/or 
SFILM and/or RADIATE and/or 
SRADIATE
Data lines to define thermal loads
CLOAD and/or DLOAD and/or DSLOAD
Data lines to define mechanical loads
FIELD
Data lines to define field variable values
END STEP