## Abstract

We present a model for approximating residual stresses arising from thermally induced plasticity and apply it to a simple two-dimensional (plane strain) test problem of transient heating. The model uses stress fields derived from simple eigenstrains to perturb an initially elastic-only solution in a way that mimics plastic flow. The plane strain condition of the test problem gives rise to plastic strains in the *out-of-plane* direction, owing to constraints on thermal expansion. However, the model is able to cope with the challenges this poses and also encapsulates the time-dependent nature of the problem. We show that the model agrees well with finite-element simulations under moderate strength heat sources provided that appropriate plastic flow rules are used. There exists sufficient generality in the model to allow its extension to more realistic scenarios.

## 1. Introduction

Thermal processes continue to play an important role in the modern day mechanics of metallic materials, be it through joining (e.g. welding) or material processing (e.g. quenching, heat treatment). An accurate knowledge of the residual stress state resulting from such processes is necessary to critically assess the integrity, or fitness-for-purpose, of an engineering component or structure. Strong transient heating in particular can have many interrelated effects upon a material, e.g. partial melting, phase transformation, plastic strain caused by constraint to localized thermal expansion and partial erasure of previous plastic strain history through annealing. As these effects may occur concurrently, it is difficult to gain a clear conceptual understanding of which are dominant with respect to the final residual stress state in any given thermal treatment, which in turn can make it hard to understand variations between apparently similarly treated materials or to know which aspects to model most carefully when simulating a process with finite-element (FE) software.

A recent European FE modelling round robin highlighted these problems. The NeT Task Group 1 (Truman & Smith 2009) was formed to examine the benchmark problem of a single weld bead laid down on the top surface of an austenitic steel plate. This weld geometry produces a strongly three-dimensional residual stress distribution, with similar characteristics to a weld repair. The single weld pass is relatively straightforward to model using FE methods, and allows full moving heat source residual stress simulations to be performed in practical time scales. However, because the weld bead is laid onto a relatively thin plate, the predicted stresses are very sensitive to key modelling assumptions such as total heat input and the thermal transient time history. The results obtained from over 10 international universities, research institutes and companies displayed significant differences. Of particular concern was the inability to form definitive conclusions concerning whether residual stresses were tensile or compressive at several key locations within the specimen.

Accurate knowledge of residual stress profiles generated by welding is a key requirement in the application of structural integrity assessment codes, such as the R6 defect assessment code (R6 Revision 4 2005). As well as detailed FE modelling, such assessment typically involves experimental measurement of the component using at least two independent measurement techniques. Measurement results often show differences between the techniques and are further complicated by the fact that measurements from different techniques are usually made at different spatial locations. When comparing and combining these measurements to determine an overall residual stress profile, the current focus is on using advanced statistical tools such as Bayesian analysis to determine a best estimate (Nadri *et al.* 2009). However, to apply these statistical tools, a ‘baseline’ model is required that predicts the typical residual stress at a given location as a function of those parameters that are poorly measured or unknown to the experimenter (e.g. exact heat input, weld duration, etc.). Constructing such a baseline for use by these tools requires tabulating to very high detail the response of a given weld model to all of these parameters.

FE models of realistic welds can take a substantial amount of time to develop (a week to define an appropriate mesh, geometry, material hardening model and heat input, for example) and a punitive amount of time to actually run (a month’s CPU time is not uncommon), making the construction of such baseline data impractical for a completely realistic model and still difficult for less ambitious simplified models. The highly non-proportional loading conditions, on the other hand, mean that developing a closed-form analytical solution to use as the baseline would be very difficult to do without introducing so many simplifications that the model loses its predictive merits. A *semianalytical* model, which starts with an analytical or else pre-tabulated *elastic* solution and then perturbs it in a way to approximate plastic flow, however, offers several potential increases in computational speed over FE modelling without necessarily requiring as many drastic simplifications as a fully analytical model would.

The potential speed gains from a semianalytical approach over that of the FE method come from two separate areas. Firstly, an FE model must calculate at each time increment the heat flux, temperature change and resulting stresses and displacements for every element in the model. A semianalytical model on the other hand already has immediate access to these quantities for all those elements in the elastic regime: either from an analytical formula or from a pre-calculated table. In addition, since the entire history of the thermoelastic solution is known at the outset, the semianalytical model can be started at the time increment at which plastic yielding first occurs rather than having to start running calculations from the moment the heat source is applied. In the case of a pre-tabulated elastic solution rather than an entirely closed-form solution, the amount of pre-tabulation may preclude being able to measure the variation to all possible parameters; however, the heat source strength, which is one of the most uncertain quantities in real specimens, can always be varied since it corresponds simply to multiplying the formula/table by a fixed factor.

The second potential speed gain of the semianalytical approach comes from the handling of plastic flow. Because plastic flow is not being accurately calculated for each element but approximated by a series of perturbations there is the possibility for this process to be computationally faster than the full plastic flow calculations. This is helped by the fact that the effect on all the other model elements of applying a perturbation to a particular element is again already known—either from a closed-form solution or from a pre-calculated table of the perturbing functions. While the results of such a model clearly cannot be as accurate as a full FE computation, because of the approximations involved, the result may still be accurate enough to create useful baseline models for experimental analysis or for other parameter/sensitivity studies that might shed light on the complications of the welding process but where the resolution or range of parameter space that needs to be covered makes FE simulation too computationally expensive.

In this paper, we present a proof-of-concept model using this semianalytical approach. Our goal is simply to determine whether such a perturbative framework is in fact capable of operating under the difficult conditions of strong and transient thermal gradients that processes such as welding involve. Future work can then investigate adapting the approach to more physically relevant situations and determining what modifications and tunings might be required to realize the potential speed increases the approach may offer. For this initial work, we thus restrict our attention to the effects of constrained thermal expansion and consequent plasticity arising from a transient heat source at the surface of a material in a two-dimensional plane strain problem. Our perturbative approach is inspired by that used by Blomerus & Hills (1998) to model plasticity arising from mechanical loading. However, we take a rather different approach in our choice for the nature of the perturbing stress fields and, owing to the thermal nature of our problem, must additionally deal with time dependence and plastic strains that are directed out of the plane of the otherwise two-dimensional problem.

The plan of the paper is as follows. In §2, we briefly present the problem to be solved, together with an overall summary of the perturbative framework in which the solution is to be implemented. The specific procedures of the model itself are presented in §3, and then detailed descriptions of the derivation of the elastic and perturbative stress fields used in the model are given in §§4 and 5. The results of the model are then presented and discussed in §6, followed by concluding remarks in §7.

## 2. Description of the problem and solution framework

### (a) The model problem

The problem we aim to model in this paper is outlined in figure 1. We consider the infinite half-space *y*<0, and assume it to consist of a single uniform material. We then apply a line source of heat to its surface, so that the line is parallel to the *z*-axis and infinitely long but of some finite thickness *w* over the *x*-axis. This line source of heat is of uniform strength *Q* and is applied at time *t*=0 for some duration before being removed. Boundary conditions are to be kept as simple as possible—traction-free and non-conducting.

This model can be thought of as an extremely simplistic and idealized representation of a weld or other transient heating problem. The symmetry and infinite extent of the problem suggest it should be treated in plane strain: that is, the total normal strain in the out-of-plane direction *ε*_{zz}, given by the sum of the elastic () and plastic () strains, should remain zero at all times.

### (b) Solution framework

The problem is to be solved by a perturbative approach. The starting point for this approach is the *purely elastic* response to the prescribed heating (derived in §4), totally ignoring the effects of any plasticity. This elastic solution is then divided into time steps. Starting from *t*=0, the solution is advanced through these time steps, comparing the stress state at each step with some chosen plastic yield criterion. If the stress state in a given step violates this yield criterion, then the stress state is perturbed by superimposing upon it additional stress fields. The aim of this perturbing stress is to produce a new stress state that approximates the elastic state of the material that would arise following a period of plastic flow sufficient to prevent violation of the yield criterion (owing to the way such methods initially allow the stress state to violate the yield criterion, this process of perturbation is sometimes referred to as ‘returning the stress state to the yield surface’). These perturbative stresses have a cumulative effect upon the stress state of later time steps and, once the thermal load has been completely removed, will remain in the material, forming the residual stress.

The application of this method to any given problem thus requires three choices: (i) the choice of a suitable yield criterion, (ii) a set of plastic flow rules, governing the expected form of the post-plastic-flow stress state, and (iii) a suitable set of stress ‘basis functions’ with which to perform the perturbations. Below, we address each of these issues in turn.

#### (i) Yield criterion

Plastic flow results in permanent plastic strains being introduced into a material in regions where the stresses would otherwise have violated the material’s yield criterion. For isotropic polycrystalline materials, the criterion usually chosen to best match the material’s true behaviour is either that of *von Mises* or that of *Tresca*.

The von Mises criterion agrees well with experiment but is a nonlinear function of the stress-tensor components; the Tresca criterion is somewhat less accurate but, provided the direction of principal stresses does not significantly change, can be written as a linear function of the stress-tensor components. Since—for efficiency reasons—we wish to restrict our solution to linear (if incremental) steps, we choose to use the latter but with a normalization chosen for better accuracy.

The normalization between the Tresca and von Mises yield criteria is somewhat arbitrary but is traditionally taken such that both criteria predict the same stress at yield under uniaxial tension. However, Blomerus *et al.* (1999) point out that if Poisson’s ratio is given by *ν*=1/2, then normalizing the Tresca criterion to match that of von Mises under pure torsion instead will render the two criteria identical under plane strain conditions; they argue that typical materials with a Poisson’s ratio of *ν*=0.3 deviate only slightly from this behaviour and so this normalization should still be preferred, and they also demonstrate its effectiveness in their own Tresca-based models of limited plasticity when compared with von Mises-based FE simulations.

This equivalence of the two yield criteria under *ν*=0.5 occurs only if there are no thermal stresses or out-of-plane plastic deformations (namely deformations with a *z*-component as defined by figure 1), since these will alter the relationship between the principal out-of-plane stress and the in-plane (the *xy*-plane as defined by figure 1) stresses. In our model, not only are thermal stresses present, but they are also sufficiently strong to cause out-of-plane plastic flow (as will be shown later in this section)—indeed, the maximum shear used by the Tresca yield criterion is rarely, if ever, found to involve just the in-plane principal stresses. Interestingly we find that, despite this, the above normalization still seems to provide an advantage over the traditional normalization when we compare our results with FE simulations, perhaps indicating that most of the plastic flow in our models is driven by stress states in which torsion dominates, and so we will use this normalization for all of the results presented in this paper.

The yield criterion effectively defines a volume of allowed stress states, bounded by a *yield surface* on which lie those stress states that have just undergone or are on the threshold of undergoing plastic flow. In the rest of this paper, when we talk of stress states that ‘exceed’ or go ‘beyond’ the yield criterion, we are referring to those stress states that lie outside this prescribed volume and so would be physically prevented by plastic flow, but which in our model must be returned to the yield surface via additional perturbative stresses.

#### (ii) Plastic flow rules

In classical plasticity theory, the components of the plastic strain increment that results from attempting to stress a material beyond its yield criterion are given by Melan (1938)
where *f* is the yield function (which the current stress state *σ*_{ij} is assumed to satisfy), *g* is the ‘plastic potential’ and *h* is an arbitrary scalar function. Since *h* and d*f* are independent of the suffices *ij*, then, for the purposes of determining the *relative* strengths of the different components of the strain increment, we can just write
(2.1)
The magnitude can then be determined by demanding that after plastic flow the (initially over yield) stress state has been returned to one lying on the yield surface.

A particularly useful choice for the potential (e.g. Hill 1950) is that of *g*=*f*. In our case, this would signify equating *g* to the yield function for the Tresca yield criterion. If we choose to orient our coordinate system so that it is always aligned with the principal axes of the stress state in the time step under consideration, then we can write
where is the maximum shear stress before yield and the three different cases correspond to different ‘planes of maximum shear’.1 Inserting this expression into equation (2.1), we obtain
(2.2)

However, there is no particular reason why we cannot instead equate *g* to the von Mises yield function, especially as we are only treating the Tresca criterion as a more convenient approximation for this criterion. In this case, substitution into equation (2.1) recovers the familiar Reuss equations (Reuss 1930)
(2.3)
(using the global coordinate system defined in figure 1). We shall consider both sets of plastic flow rules when implementing our model and compare their effectiveness when presenting our results (§6).

#### (iii) Choice of basis functions

Since plastic flow results in permanent deformations of the material that are not necessarily compatible with its original reference configuration, we must choose as our basis functions for perturbing the elastic solution stress fields which themselves *do not* satisfy the compatibility conditions in order to reconstruct the final (post-plastic flow) stress state of the material.

We have several choices in this matter. The set of stress fields resulting from idealized material dislocations represent one choice since, by their very nature, dislocations represent a physical discontinuity in the displacement field; and these have indeed been used successfully to model limited plasticity in the past (Blomerus & Hills 1998). Another choice is the Kelvin point force stress solution2 (or equally a doublet of such forces), which violates compatibility at the point the force is applied, and has also been successfully applied to plasticity problems (Chen & Nisitani 1985). Perhaps the most direct route, however, is to use stress fields specifically derived from the insertion of a permanent strain into an otherwise uniform material—so-called *eigenstrain* solutions (e.g. Eshelby 1961; Mura 1991).

Calculations with FE simulations3 show that the model of figure 1 results in a residual stress state that has a significant out-of-plane elastic strain (as shown in figure 2). Since the assumed plane strain condition means that the *total* out-of-plane strain must be zero, this in turn implies that there must have been significant out-of-plane plastic strain. This is normally unusual in purely mechanical two-dimensional problems, since the value of Poisson’s ratio for most materials means that the out-of-plane stress during loading typically lies between the two other principal stresses (i.e. it is neither the largest nor smallest of the principal stresses), with the result that plastic yielding occurs in-plane. Under thermal loading, however, the out-of-plane stress is governed not just by Poisson’s ratio but also by the resistance to thermal expansion in the out-of-plane direction implied by the plane strain condition, allowing the out-of-plane stress to assume much larger values. This affects how easily the three different stress basis functions discussed above can be implemented.

Using eigenstrain stress solutions from the literature, it is relatively straightforward to devise basis functions that are both independent of *z* position *and* correspond to non-zero out-of-plane plastic strain, hence giving rise to out-of-plane elastic strains. For instance, cylinders of constant permanent strain orientated parallel to the *z*-axis, or periodic permanent strain fields that vary only in *x* and *y*, can be constructed with out-of-plane permanent strains. On the other hand, the task is much more difficult using either point forces or dislocations. In this case, neither the two-dimensional point force/force doublet nor the two-dimensional dislocation solutions (either edge or screw) produce normal out-of-plane elastic strains. Thus, one must start from a fully three-dimensional solution such as a general force doublet or an infinitesimal loop dislocation, and then attempt to somehow remove the dependence on the *z*-coordinate by, for example, constructing an infinite line of such sources parallel to the *z*-axis. However, since these solutions are incompatible only at their point of insertion, the only contribution to the out-of-plane strain from such sources occurs *on* the singular insertion point. To get a net contribution of out-of-plane strain in some finite region of the *xy*-plane would thus seem to require considering a *distribution* of such sources and so introduces the additional complexities of handling the resulting singular integrals.

For these reasons, we shall use the more easily adapted eigenstrain-based basis functions (which we derive in detail in §5) in the remainder of this work, deferring consideration of other basis functions to future work.

## 3. Method of solution

The procedure followed by our model is thus as follows. The half-space is divided into square elements (in fact because of the symmetry of the problem, we need consider only the quarter-space). The initial underlying purely elastic solution is then broken into small time steps. Starting from *t*=0, each time step is examined in detail to see whether the stress state in any of the square elements exceeds the yield criterion at that moment. Once a time step is reached, where yielding does occur in some of the elements, an eigenstrain (namely permanent plastic strain) with components determined by the plastic flow rule, but unknown magnitude, is assigned to each yielded element. (For simplicity, the stress is examined only at the centre of each element, and the inserted eigenstrain is uniform across the element.) A set of simultaneous equations is then constructed based on the demand that the new stress state at each yielded element be reduced to lie exactly on the yield surface as a result of the influence of the eigenstrain imposed at that and every other yielded element (which for the Tresca yield criterion means reducing the maximum shear stress to a critical value ). Since there are as many unknown eigenstrain magnitudes as there are yielded elements, and since the stress at any point owing to the eigenstrain in a given element is linearly proportional to the eigenstrain magnitude (see §5 for details), these equations are linear and can be solved by a single matrix equation. Specifically,
(3.1)
where is the (over-yield) shear stress felt by the *i*th element in the current time step, before any insertion of eigenstrain; the eigenstrain, of unknown magnitude *m*_{j}, which the plastic flow rule has determined for the *j*th grid element is *m*_{j}*ε*^{pl,j}; and is the perturbing shear stress felt at the centre of the *i*th element (*x*^{i}) owing to the insertion of a square of constant eigenstrain of form *ε*^{pl,j} in the *j*th element. The indices *i* and *j* range over all over-yield elements, so that the final shear stress in these elements after perturbation exactly equals the maximum allowable under the Tresca yield criterion, . The calculation must be carried out in the coordinate system *x*′,*y*′ in which the shear component of the stress tensor *σ*^{0,i} is maximized. This coordinate system is simply related to that of the principal stress axes (although in general *x*′,*y*′ will not lie in the *xy*-plane).

The solution is then advanced to the next time step, only this time the elastic solution is augmented by the previously determined eigenstrain-induced stress field. If any elements are now found to have exceeded the yield criterion (possibly for the second or more time), then the above process is repeated again to find a *new* eigenstrain-induced perturbative stress field that is added to the previously determined perturbations and also to the elastic solution which is then advanced to the next time step again and so on until the solution has been advanced sufficiently far in time that no further plastic flow will occur (i.e. when most or all of the applied heat has dissipated).

In practice, the nonlinear nature of the plastic flow process cannot be completely ignored. In particular the perturbative stress field may cause a significant change in the principal stress directions, which then invalidates the solution (since it will alter the orientation of the *x*′,*y*′ coordinate system required for the left-hand side of equation (3.1) to equal the maximum shear stress of the system). If it is detected during a time step that this has occurred, then the matrix solution is retried iteratively, using instead for the principal stress directions those directions that were found to be the direction of principal stress at the *end* (i.e. after insertion of the eigenstrain) of the previous pass. This process is continued until convergence is reached, as measured by there being little or no further change in these principal stress directions.

Another possibility that can require re-evaluation of the solution is if the calculated perturbing stresses cause neighbouring elements that were previously unyielded to exceed the yield criterion. In this case, we must increase the size of the matrix to additionally include these new elements and then re-solve the perturbative stress for the current time step using the new, larger set of elements.

The results achieved with this model are presented and discussed in §6, but first in the following two sections we present the explicit derivation of the thermoelastic solution and the perturbing basis functions used by the above algorithm.

## 4. Derivation: the thermoelastic solution

To calculate the thermoelastic solution for the problem posed in §2*a*, consisting of the temperature (*T*) and stress (*σ*_{ij}) fields as a function of time *t* and position (*x*,*y*), we start from the plane strain solution for an instantaneous and infinitely thin line source of heat in the *constrained* half-space (e.g. Barber & Martin-Moran 1982), which we use as an effective Green’s function,
written here with respect to a heat source placed along the line *x*=*y*=0 and activated at time *t*=0 with strength *Q*. *E* represents Young’s modulus of the material, *c* the specific heat capacity, *k* the thermal diffusivity, *α* the thermal expansivity, *ρ* the density, *ν* Poisson’s ratio, and we have written *r*^{2} for *x*^{2}+*y*^{2}. Since we are assuming all material properties to be independent of temperature in this simple model, the temperature field *T* can just be thought of as the temperature *difference* from ambient rather than being linked to a specific zero point. Lastly, we note that plane strain conditions mean that the out-of-plane stress *σ*_{zz} will be given by
(4.1)

To obtain the final required solution, that of a heat source of finite duration and finite width, we must integrate this solution over space and time as follows:
and
where *σ*^{cor} represents the additional correction factor (discussed below) needed to change from a constrained to an unconstrained half-space. It turns out that the integral over time can be evaluated in terms of the exponential integral function
for which quickly converging numerical series are available. This results in the form
(4.2)
where, for ,
and, for ,
with still being given by the plane strain conditions (equation (4.1)).

The correction terms *σ*^{cor} needed to move from the constrained solution, where the surface *y*=0 is constrained to remain exactly plane (i.e. no displacement), to the unconstrained solution, where the surface *y*=0 is traction free, can be derived by distributing point forces (given by the well-known Flamant solution; e.g. Barber 2002) along the *y*=0 surface with magnitudes that exactly cancel the tractions existing in the constrained solution. The constrained solution already has no shear tractions along *y*=0, so we need only to cancel the normal tractions in this way, giving
where is the traction present in the constrained solution, found by substituting *y*=0 into equation (4.2) *without* any correction term. We note that, although the integrals in equation (4.2) do not in general seem to have closed form solutions, for the specific case of *y*=0 we are able to find a closed-form for the tractions, namely for the case ,
where *x*_{↓}=*x*−*w*/2, and *x*_{↑}=*x*+*w*/2. (The case is similar but with like terms in being subtracted from the solution.)

In this way, all the terms of *σ*^{final} (and *T*^{final}) are reduced to terms involving only single integrals. While these remaining integrals do not seem to have closed-form solutions, they can nevertheless be efficiently calculated numerically—in our work, we have chosen to use an algorithm based on the double exponential integral formulation (e.g. Takahasi & Mori 1974; Mori 2005). The resulting stress and temperature fields are validated against an FE simulation in figures 3 and 4. Finally, since we only require the stress and temperature solutions on a set number of grid points, we further speed up our calculations by pre-computing and tabulating the stress and temperature at these points for a range of different times and then simply use linear interpolation to find the values at any specific time required in our calculations.

We note that Azarkhin & Barber (1988) have derived the solution for a similar problem in the form of a series solution, from which we could equally have formulated a solution to the problem discussed here. However, we find it more convenient to use the method presented here as it avoids the problem of determining the precise values of *x*, *y* and *t* where one can best make a smooth transition between the two forms of the series solution presented by Azarkhin & Barber (1988)—one for small values of (*x*^{2}+*y*^{2})/(4*kt*) and one for large. Such a smooth transition would be desirable to ensure that there are no sudden jumps in the size of the region that is over the threshold for plastic flow between time steps. Since our own method requires the use of numerical integration, probably it is slower than that of Azarkhin & Barber (1988), but this consideration is of less importance in our case, where results are to be calculated only once and then tabulated.

## 5. Derivation: the perturbing basis functions

As already discussed, eigenstrain will be inserted into each of the elements in the form of a uniform field across the element. The required basis functions thus consist of the set of stress response fields as a function of position ** x**=(

*x*,

*y*) owing to the insertion of a square of constant eigenstrain at the position of the

*pq*th grid element. We find that the most convenient way to construct these stress response functions is by their Fourier series, thus the derivation of the basis functions can be split into two parts: the construction of the response of the half-space to a given periodic eigenstrain field, and then the construction of the Fourier series for a square element of constant eigenstrain.

### (a) The harmonic response function

The stress response to a harmonic eigenstrain field in an infinite (three-dimensional) space is readily available in the literature (e.g. Mura 1991) and for a harmonic eigenstrain *ε*^{pl}(** x**) of amplitude and wavevector

**=(**

*ξ**ξ*

_{x},

*ξ*

_{y},

*ξ*

_{z}) defined by is given by where

*C*

_{ijkl}is the elastic stiffness tensor, and with

*ϵ*

_{ijk}being the usual Levi–Civita symbol.

To adapt this to our particular problem, there are two steps we need to take. Firstly, we should simplify the stress function to reflect the symmetry of the problem we are dealing with. We can accomplish this by setting *ξ*_{z}=0 to indicate that the stresses and strains in our problem are independent of the *z* coordinate and by setting to reflect the fact that (because of symmetry) we always expect the *z*-direction to be one of the principal strain axes. Inserting these conditions, together with Hooke’s Law for an isotropic medium, into the above, we get the stress field
where
and *λ* and *μ* are the usual Lamé’s constants. Note that we still have *σ*_{zz}≠*ν*(*σ*_{xx}+*σ*_{yy}), and thus these stress functions still admit a non-zero out-of-plane elastic strain.

The second modification we need to make is to correct these solutions so that they apply to the infinite half-space *y*<0 rather than to a fully infinite space. We can accomplish this by distributing cancelling point forces along *y*=0, just as was done in §4 for the thermoelastic solution, only here we must cancel both normal and shear tractions. It turns out that the convolution integrals required to do this are particularly simple to evaluate in this case. This is because we can carry them out in Fourier space, where our uncorrected stress field becomes simply a Dirac delta function (making the required inverse transform trivial). The final corrected stress function for a half-space, *σ*^{hs}, then becomes
(5.1)
where
and *B*(** ξ**) and

*C*(

**) are as before.**

*ξ*### (b) Creating and transforming the Fourier series

Having established the stress response to harmonic eigenstrain fields, we move on to calculate the Fourier series of a square element of constant eigenstrain and hence its own stress response. We wish to express the eigenstrain field in terms of a two-dimensional Fourier series
(5.2)
where the integers *m* and *n* index the possible frequencies along the *x*- and *y*-directions respectively. Obviously, since we will not be interested in the entire infinite half-space, but rather only those regions near the applied heat source, we will be able to truncate the region this Fourier series is constructed in to some finite region *L*_{x}×*L*_{y} (although this will necessarily involve some loss of accuracy, as discussed later) so that ** ξ**=(

*ξ*

_{x},

*ξ*

_{y})=(2

*πm*/

*L*

_{x},2

*πn*/

*L*

_{y}).

Since each basis function is composed from the eigenstrain in only a single square element (which we shall label the *pq*th element), and since this eigenstrain is constant across this element, the coefficients of this Fourier series become very simple to calculate
(5.3)
where the value of the (constant) eigenstrain found in the *pq*th element is represented by . The remaining integral over the element can be written
where (*x*_{0},*y*_{0}) are the coordinates of the start (the lower left corner as viewed looking down the *z*-axis) of the truncated region of the half-space, which we divide into *N*_{x}×*N*_{y} square elements. The integers *p* and *q*, respectively, label the *x* and *y* positions of the element into which the eigenstrain is inserted and can range from 1 to *N*_{x} and *N*_{y}. The integral is easily evaluated in a closed form, although care must be taken to treat the cases where *n* or *m* is zero separately.

To get the corresponding stress response, which can then be used as the perturbative stress in the solution matrix of our algorithm (equation (3.1)), we simply sum the stress responses from each of the frequency components
(5.4)
where the value of corresponding to *ε*^{pl}[*pq*] for any given wavevector is determined by equation (5.3) and *σ*^{hs} is the half-space stress response to an eigenstrain of given wavevector ** ξ** and magnitude as formulated in equation (5.1). Explicit forms for the expressions following these substitutions are given in the appendix (available in the electronic supplementary material), since they are somewhat lengthy.

The zeroth-order (*m*=*n*=0) Fourier coefficient is not used in the sum since it corresponds to an eigenstrain field that is constant everywhere and hence (since it satisfies the compatibility equations on its own without the need for a non-zero elastic strain) produces no stress. In practice, these stress responses will be calculated only once and then tabulated for subsequent use in the model.

We note that, as asserted in our earlier discussion, these response functions depend only linearly on the magnitude of the inserted eigenstrain.

### (c) Some practicalities

Although we have now outlined the main derivation of our basis functions, a few details remain that need to be addressed.

Since our model is symmetric across *x*=0, it is convenient to solve only in the quarter-space. In doing so, it should be noted that any eigenstrain inserted into the right-hand quarter-space will in reality be matched by a similar insertion in the left-hand quarter-space, and the resulting stress response must take account of this. Thus, although ultimately the iterative solution process embodied in equation (3.1) is carried out only over elements in the quarter-space (which make up what we shall call the ‘model grid’), the eigenstrain Fourier series used to construct the response functions needs to extend across the *half*-space in order to include eigenstrain contributions from both sides of the symmetry line. This and other processes (discussed below) which contribute to the size of the *L*_{x}×*L*_{y} region over which the Fourier series is constructed (which we refer to as the ‘Fourier grid’) are shown in figure 5.

As already discussed, we are concerned only with the behaviour of the solution in a truncated region of interest near the heat source (i.e. the model grid is finite), since far from the heat source, the stress and temperature fields become very small. Thus, the Fourier series of the eigenstrain, too, is constructed over only a finite region. This has important consequences, since the stress field resulting from this eigenstrain will itself have a Fourier component, and this component will have a periodic continuation outside of the truncated region of immediate interest. This has two effects upon the accuracy of the constructed stress field.

The first effect concerns additional unwanted constraints placed upon the possible form of the stress field. For there to be no far-field mechanical load on the half-space, the mean value of over the half-space must be everywhere finite (and similarly ). However, if a component of *σ*_{xx} is periodic, then, to satisfy the above, the integral of this component of *σ*_{xx} over one period must be zero. Thus, by truncating the region over which the Fourier series is taken and so giving it a finite period, we introduce an additional artificial condition upon the stress field. The smaller the truncation region, the more it will affect the accuracy of the predicted stress field. Note that if we extend the *L*_{x}×*L*_{y} Fourier grid into regions where eigenstrain will never be inserted, in order to minimize these inaccuracies, there is no strict need to do so in increments the size of the square elements used within the main region; however, it is much simpler if this is the case (so that the Fourier grid is always spanned by *N*_{x}×*N*_{y} square elements), and has been implicitly assumed in the formulation presented in this paper.

The second effect is due specifically to the periodicity along the *y*-axis. The stress response functions are by their nature inherently discontinuous, since they represent the response to a single element of uniform eigenstrain surrounded by nothing else. However, when added together to form the final solution, the resulting stress is expected to be only weakly discontinuous since the series of eigenstrains inserted in any given time step is expected to decay gradually away from the centre of the heat source. Periodicity in any component of the stress along the *y*-direction, however, will introduce a strong discontinuity into the stress field, since it effectively ‘wraps’ the surface (*y*=0) of the half-space, where eigenstrain is expected to be largest, on to the bottom of the truncated region, where eigenstrain is not expected at all. This also potentially affects the accuracy of the resulting stress field.

To mitigate the latter problem, we in fact choose to extend the eigenstrain field symmetrically to the other side of the half-space surface (even though *physically* there is no medium beyond the surface boundary), taking the Fourier series over this new larger region. This does not violate the validity of the stress fields since the eigenstrain *inside* the half-space is still of the desired form and the boundary itself remains traction free owing to the correction terms inserted into the stress field.

The minimum size of the region over which we can construct our eigenstrain Fourier series is thus twice as wide (to allow symmetrical insertion of eigenstrain into the left-hand quarter-space) and twice as high (to allow the mirroring of the eigenstrain field along the *y* direction discussed above) as the truncated region in which our final solution is actually employed. Although, in the interests of accuracy, choosing a region still larger than this is preferable as it minimizes the additional constraints on the form of the stress field imposed by periodicity. All these points have been summarized in figure 5.

The reason we cannot choose to make the Fourier grid (as parametrized by *N*_{x} and *N*_{y}) arbitrarily large is computational, since it affects how many frequencies need be included in the Fourier series and thus how readily the final stress function can be calculated. As written in equation (5.2), the eigenstrain Fourier series is in fact summed over all possible frequencies (i.e. ). However, since the resulting stress function will be sampled only at the centres of the square elements in the ‘model grid’, which we assume to be representative of the region as a whole, we in practice ignore those frequencies expected to describe fine scale detail on regions smaller than the grid elements. This in turn corresponds to excluding frequencies in excess of *m*=±*N*_{x}/2 and *n*=±*N*_{y}/2, respectively.

## 6. Results

In this section, we now present and discuss the predicted stresses that result when the model formulated in the last few sections is applied to the initial problem outlined in §2*a* and figure 1. We start though with a brief discussion concerning the material parameters used when running our model and the FE simulation that has been created to validate our results.

### (a) Material properties

The thermoelastic stress arising in the problem of figure 1 is dependent upon a number of material properties, and not all of these dependencies can be easily absorbed into dimensionless parameters. Thus, to facilitate the presentation of results, and also because we wish to validate these results against FE simulations, we assign illustrative values to these properties, given in table 1. These values are mostly representative of steel, although some stainless steels can have significantly lower diffusivity.

We further set *w*=0.01 m and s. In a truly realistic scenario, the yield stress would be a function of the material’s temperature (as indeed would many of the material properties listed in table 1). However, for simplicity, here we consider all material properties as temperature independent. As a result of this simplification, the values of *Q* and the yield stress are not in fact independent: increasing *Q* has exactly the same effect as reducing the yield stress, namely it increases the total extent of plastic flow and causes the initial onset of plasticity to move towards earlier times. For the current validation, we simply tune these parameters to ensure that there is a large plastic region created. Resulting stresses are then quoted in units of the von Mises yield stress, *σ*_{Y}.

### (b) Finite-element simulation

To provide us with a benchmark against which to validate our model’s results, we have constructed a simple FE simulation of the problem specified in figure 1 using the material properties stated in table 1. This simulation is constructed using the off-the-shelf FE package Abaqus (Dassault Systèmes 2007). The infinite quarter-space is represented by a 1×1 m plate, with symmetric boundary conditions along *x*=0 and traction-free boundary conditions on the remaining edges. The 10×10 mm region nearest the heat source is divided into 76×76 square elements (roughly 0.13 mm each in length), with larger elements being used further from the heat source for a total of 10 446 elements.

Using this set-up, a coupled thermal and stress simulation was carried out with linear, reduced integration point, plane-strain elements and a von Mises yield criterion under elastic–perfectly plastic conditions. In addition, heat was allowed to dissipate from the non-heated edges of the plate, to simulate the infinite extent of the idealized problem. The model and simulation are compared in the following section.

### (c) Main results

Using the material parameters discussed above, we ran two instances of our model: one using the Tresca-derived plastic flow rules of equation (2.2) and one using the standard Reuss flow rules (equation (2.3)). Stresses were calculated at the centre of each square element in a 76 × 76 grid covering the region of the quarter-space nearest the heat source and coincident with the grid used by the FE simulation, grid elements being of length 0.13158 mm to a side. The eigenstrain basis functions were calculated using *N*_{x}=*N*_{y}=304 (decreasing these values generally results in an overall decreased accuracy), and fixed time increments of 0.02 s were used up until *t*=2 s. The final state, , was achieved in just one additional increment.

The predicted residual stresses remaining in these two models after the heat source has been removed and the temperature has returned to ambient are shown in figures 6 (Tresca-based flow rules) and 7 (Reuss rules), plotted against the residual stresses found from the FE simulation described in §6*b*. Perhaps unsurprisingly, we see that the model based on the Reuss rules seems to perform better than that using Tresca-based plastic flow rules.

Specifically, we can see that, along the lines shown, the model based on Tresca flow rules fails to correctly reproduce the *σ*_{xx} component of stress. Other stress components are more accurately produced, in particular *σ*_{zz}. This suggests that while this set of flow rules produces roughly the correct amount of plastic strain in the out-of-plane direction, its distribution of the in-plane components of plastic strain does not provide a sufficiently good approximation for accurate prediction of residual stresses.

The model based on the Reuss flow rules, on the other hand, seems relatively successful at reproducing all the stress components. Although data are only shown in the figures for lines near the surface and the symmetry axis, we have checked the results along other lines and find a similarly good agreement with the FE simulation. It is interesting that although the Tresca *yield criterion* (subject to the normalization discussed in §2*b*(i)) provides a good approximation to the von Mises criterion used in the FE simulation, yet Tresca-based plastic flow rules do not seem to be able to substitute for the von Mises-based Reuss rules.

Lastly, we also show in figure 8 the set of grid elements that have experienced plastic flow (a total of 621 in our Reuss-based model), superimposing upon this the residual-stress contours. The extent of the plastic region closely matches that which would be predicted from an extrapolation of the purely elastic stress field derived in §4 (although of course residual stresses could not be predicted from such an extrapolation). In future work, we hope to investigate the ability of our model to converge and make predictions when the strength of the heat source is further increased to produce plastic regions that differ markedly in shape from purely elastic predictions. Other extensions that ought to prove amenable to future investigation are a temperature-dependent yield stress (since it only affects the right-hand side of equation (3.1)) and strain hardening (since our algorithm already acquires the plastic strain history of the material as the thermal load is applied).

## 7. Conclusions

We have implemented an adaptive semianalytical model for stresses arising from thermally induced plasticity by using basic eigenstrain solutions to perturb an initially perfectly elastic solution. The resulting model was applied to a simple test problem motivated by welding and was shown to agree well with FE simulations under a heat source strong enough to cause plastic flow in 621 of the model elements. It was shown that, while a Tresca yield criterion can prove sufficient for accurate results (at least with this test problem), the plastic flow rules cannot be similarly simplified. Further investigations are planned to determine whether the model is also able to function under the very high heating regimes which would give rise to plastic zones that differ significantly in shape from those predicted by elastic extrapolation.

The eigenstrain basis functions we have derived for this model are independent of the thermoelastic solutions, and so adapting the model to other thermal problems in the same geometry should just be a case of slotting a new elastic solution into our existing framework. There is also scope within the model for extensions to strain hardening and a temperature-dependent yield stress.

## Acknowledgements

J.M.B. would like to thank the Japan Society for the Promotion of Science (JSPS) for funding a research visit to Nagasaki University Mechanical Systems Engineering Department where part of this work was completed and also thanks all those at Nagasaki University for their hospitality.

## Footnotes

↵1 Note that we could have chosen to formulate this in terms of, for example, (

*σ*_{1}−*σ*_{2})^{4}instead, but the final result would have remained unchanged.↵2 See, for example, Barber (2002).

↵3 A full description of the FE simulations, including the input parameters, used for verification of our results is given in §6. Figure 2 is included for illustrative purposes.

- Received February 13, 2009.
- Accepted November 2, 2009.

- © 2009 The Royal Society