Open access

The Use of Matlab in Advanced Design of Bonded and Welded Joints

Written By

Paolo Ferro

Submitted: November 15th, 2010 Published: September 9th, 2011

DOI: 10.5772/21601

Chapter metrics overview

6,227 Chapter Downloads

View Full Metrics

1. Introduction

From a mathematical viewpoint, welding can be considered as a transient boundary problem in which the thermal input varies in space and time. Thermal coefficients involved with heat lost by convection and radiation, are usually temperature-dependent; displacement constraint conditions are also imposed on the geometry. Temperature and stress development in the material are steered by the well know differential equations of heat exchange and elasto-plastic equilibrium. Two complications in the models are due to the microstructure transformations (which influence the mechanical behaviour of the joint), and to the filler metal (which influences also the chemical composition of the parent metal and the temperature distribution). It should be noted that coupled phenomena are involved because the latent heat of phase transformation influences the temperature distribution due to the welding source. Moreover, the constitutive relations which connect stresses to strains are both temperature and phase dependent.

The development of new models for joint planning is of great importance in the industrial and research fields. The prediction of residual stresses, temperature distributions, phase transformations, asymptotic stress fields near the weld toe or near the interface between the matrix and the adhesive in bonded joints, may be of fundamental importance for a good joining plan and operation. In this chapter, some models developed for bonded and welded joints and solved by means of Matlab program, are presented. In the first part, models for temperture ditributions and phase transformations diagrams are considered with particular attention, referred both to conventional and innovative welding precesses, such as laser and friction stir welding. In the second part, mechanical models are described with particolar attention put on residual stress calculation and advanced joints planning methodologies.

Only analytical or semi-analyitical models are taken into account due to their efficiency compared to Finite Element models. As a matter of fact, by using analytical models for temperature distributions prediction it is possible to optimize the process parameters such as power source, welding speed and pre-heating temperature, with low effort in terms of time and cost. Such models offer also the possibility to predict the fusion zone (FZ) and heat affected zone (HAZ) extension and, finally, to perform a parametric sudy of welding. In this work, the Rosenthl solution (Rosenthal, 1941; Rosenthal and Shamerber, 1938) of the welding thermal problem, will be described. Moving point source, linear source or combinations of the last two are used to reproduce the fusion zone shape of the joint, leading to good results in a very short time. Referring to the residual stresses calcualtion, the ‘’rod model“ by Cañas et al. will be presented. In this model, the equations referred to the problem under consideration, are written in matrix notation and advantageously used in an efficient algorithm solved with Matlab. As refers to phase transformations, the possibility to calculate the Time-Temperature Transformation (TTT) and Continuous Cooling Transformation (CCT) diagrams of steels with use of the Kirkaldy model is also proved.

In the matter of advanced joint planning, it was shonwn that better predictions of static and fatigue resitance are possible if the intensity of the asymptotic stress fields is taken into account. Dealing with the fatigue strength of welded joints, such local approach models the weld toe region as a sharp, zero radius, V-shaped notch. Under these conditions, the intensity of asymptotic stress distribuions, obeying Williams’ solution, are quantified by means of the notch stress intensity factors (NSIFs). The formulation of this method is completed analytically and the resulting set of ordinary differential equations is solved numerically by means of Matlab. In this chapter it is described with a particular attention put on bonded joints.


2. Thermo-metallurgical analysis

The thermal history induced by a welding process can be calculated by solving the fundamental equation of heat transfer (1):

ρ C p T ˙ = k T i j L i j ( T ) p i j E1

where ρ is the material density, C p is the specific heat capacity of the material (weighted according to proportions of various phases), k is the thermal conductivity of material (weighted according to proportions of various phases), T is the temperature, T ˙ is the temperature rate (Newton’s notation is used for the time derivative of a function), L ij (T) is the latent heat (at temperature T) of the i→j transformation, p ij is the phase proportion of i-th phase which is transformed into j-th phase in the time unit, and

= i 1 x 1 + i 2 x 2 + i 3 x 3 E2

is the 3D gradient vector operator. The heat transfer boundary conditions of the problem are

q = k T E3

where q is the heat flux at the boundary which, in the welding process, consists of a prescribed function of the time and space (heat source), convective and radiative heat loss, and zero flux in a symmetry plane.

The Rosenthal solution of Eq. (1) is described both for fusion and friction stir welding with some guidelines put on its practical use in welding process plan.

2.1. Temperature distribution and cooling rate in fusion welding

The analytical solution of Eq. (1) was given by Rosenthal (1941), who considered a point source moving on a semi-infinite plate under steady-state conditions, with temperature-independent material properties, at convective and radiative heat loss and phase transformations neglected. In a reference system linked to the source, this solution is given by the following equation:

T T 0 = Q 2 π k e λ v ξ e λ v R R E4

where T 0 is the reference temperature, R=(ξ 2 + y 2 + z 2 ) 1/2 is the radial distance of a point of the plate from the source axis, v is the welding speed, t is the time, ξ = x – vt is a moving co-ordinate, λ = 1/(2α) (where α is the diffusivity), and Q is the effective thermal power absorbed by the material. In the case of a line source in a plate of thickness H, the relation (3) becomes:

T T 0 = Q 2 π k e λ v ξ K 0 ( λ v r ) H E5

where K 0 is the modified Bessel function of the second kind and zero order. In the case of arc welding, the effective thermal power equals to

Q = η V I E6

where V is the arc voltage, I is the current intensity, and η is the arc efficiency. Eq. (4) can be easily solved by means of Matlab. Figure 1 shows a 3D representation of Eq. (4) and a comparison between analytical and Finite Element (FE) solution (Ferro et al., 2002). It can be observed that the asymptotic solution at the centre of the source gives non-realistic results (asymptotic temperature distribution). This means that Eqs. (3) and (4) are valid only when referred to the points distant from the heat source axis.

Figure 1.

Thermal analytical solution (Rosenthal, 1941) (a); comparison between FE and analytical thermal solution (b). Material: AA-5083-O, welding technology: GMAW, voltage = 23.4 V, current = 170 A, welding speed: 11 mm/s

The time derivative of Eq. (3) gives an estimation of the cooling rate (at the point of its maximum value (y = z = 0)) and ξ< 0):

T t = 2 π k Q v ( T T 0 ) 2 E7

Equation (6) shows the strong (squared) dependence of the pre-heating (T-T 0) on cooling rate compared to the others process parameters such as welding speed and power. It means that pre-heating is the most efficient variable, in the case of steel welding, which can be modified in order to obtain sound welds. Knowing the critical cooling rate of the steel, it is possible to estimate, by using Eq. (6), the pre-heating value needed to avoid a martensitic microstructure in the weld bead.

2.2. Temperature distribution in Friction Stir Welding

The main difficulty in the formulation of any model for friction stir welding (FSW) is due to the high coupling between thermal and mechanical phenomena. Thus, in Ferro et al. (2010), the solution of the formulated equation was obtained by a numerical routine written in Matlab code under the simplification of isothermal condition at the matrix/tool interface. The formulation of heat flow is based on Rosethal’s solution, while the heat generation is described as a surface flux, which depends on the variation of the shear yield stress with temperature.

2.2.1. Governing equations

As proposed by P. Vilaca et al. (2007), the heat source can be considered to be concentrated at the mid-thickness of the plate, simulating the typical location of the nugget centre and travelling with a constant linear velocity (v) (fig. 2)(Eq. 3).

Figure 2.

Point Heat Source located at the nugget centre.

In this case, in Eq. (3), Q is the total heat generation due to frictional and plastic dissipation. Eq. (3) gives good results in case of thick plates. For medium thick plates, good results can be obtained by using eq. (7) (Ferro et al. (2010)):

T T 0 = Q 2 π k e v ξ / 2 α ( i = i = + 1 r i e v 2 α r i ) E8

where r i = ξ 2 + y 2 + ( z 2 i H ) 2 and H is the thickness of the plates. Finally, if thin plates has to be modelled, the best results can be found by using linear heat source (Eq. (4)) instead of point-source.

The heat flux in FSW is primarily generated by the friction and the deformation process. However, the mechanical loads applied by the pin tool to the workpiece result in a yielded region only in the immediate vicinity of the former while most of the workpiece remains unyielded. Thus, idealizing the localized yielded region as being coincident with the tool surface and treating the workpiece as a rigid matter, the heat generation originated from both frictional and plastic dissipation, can be modelled via surface flux boundary condition at tool/matrix interface (Schmidt et al., 2008, Perivilli et al., 2008). It can be found (Ferro et al., 2010) that the heat generated is expressed by the formula:

Q w o r k p i e c e = 2 3 π η ω τ 0 ( 1 T * T M ) [ ( R s h 3 R p 3 ) ( 1 + tan β ) + R p 3 + 3 R p 2 H p ] E9

where β is the tool shoulder cone angle, H p is the tool pin height, R sh and R p are the shoulder and pin radius respectively (Fig. 3), T M is the melting temperature, τ0 is a fitting material parameter, η is the thermal efficiency of the process, ω is the angular velocity of the tool and T* is the temperature at the tool/matrix interface. Thus, in the case of thick plates, the temperature field induced by friction stir welding under steady-state conditions can be described by eq. (9):

T = T 0 + 2 3 π η ω τ 0 ( 1 T * T M ) [ ( R s h 3 R p 3 ) ( 1 + tan β ) + R p 3 + 3 R p 2 H p ] 2 π k e v ξ / 2 α e v 2 α r r E10

Figure 3.

Schematic representation of the tool geometry.

(Ferro et al., 2010). The unknown parameters in Eqs. (8) and (9) are: η and T*, showing the coupled thermo-mechanical characteristic of the problem under consideration. However, Eq. (9) can be solved by using a simple numerical routine, as one described in the flow chart of Fig. 5 where ΔT is the temperature increment and T* trial is an interface trial temperature value. This routine was written in Matlab code. Finally, the thermal efficiency ( η ) is calculated by a reverse analysis as in any other analytical or phenomenological numerical model of welding process. Because of the temperature singualrity of the Rosenthal’s solution, the validity of Eq. (9) is limited to a zone sufficentely far from the source centre. Moreover, in order to obtain a numerical solution of Eq. (9), in the proposed model, T* refers to the temperature reached at a distance R* close to R sh (Fig. 4). Good results were obtained with R*=R sh for aluminum alloys. The modelling procedure, although not perfect, is belived to be a reasonable approach. Finally, since the Rosenthal solution depends on the thickness of the plates (Eqs. (3, 4,7 )), the correct formulation has to be used according to the analysed plates thickness.

Figure 4.

Schematic representation of the tool/matrix interface temperature and R* parameter

Figure 5.

Flow chart of the numerical routine written in Matlab code

The model was validated by comparison with different experimental data found in literature (Ferro et al., 2010). In Fig. 6 it can be observed that the thermal history calculated by the model is in good agreement with that measured experimentally.

Figure 6.

Variation of transient temperature at different locations of the thermocouples for rotational speed of 240 rpm, welding speed 3.32 mm/s (η=0.5, R*=R sh, 0=51.56 MPa) (lines: Semi-analytical solution, symbols: test data (Chao et al. (2003)))

2.3. TTT and CCT diagrams calculation

The Time-Temperature Transformation (TTT) and Continuous Cooling Transformation (CCT) diagrams are useful tools in thermomechanical processing of steels. Such diagrams depend on a so great number of variables that it is impossible to produce enough experimental diagrams for generalised use. For this reason, significant work has been undertaken to develop models that can calculate TTT and CCT diagrams for steels.

Starting from Kirkaldy’s model, the general formulation of TTT diagrams is described by the relation (10):

τ ( X , T ) = F ( C , M n , S i , N i , C r , M o , N ) Δ T n exp ( Q e f f / R T ) S ( X ) E11

where τ is the time needed to transform X volume fraction of austenite, T is the temperature, F is a function of steel composition (expressed in wt%), N is the prior austenite grain size (ASTM number), ΔT is the undercooling, Qeff is the effective activation energy for diffusion, and the exponent n is an empirical constant, determined by the effective diffusion mechanism; n = 2 for volume and n = 3 for boundary diffusion. S(X) is the reaction rate term, which approximates the sigmoidal effect of phase transformation. In the work of Victor Li at al. (1998), S(X) is expressed as follows:

S ( X ) = 0 X d X X 0.4 ( 1 X ) ( 1 X ) 0.4 X E12

In a TTT diagram, the location of the ‘nose’ of each C curve correlates to the maximum reaction rate. The exact locations of the nose are jointly determined by the values of n and Q. From Eq. (10), at the nose temperature the denominator has to be maximum, thus:

d d T ( Δ T n exp ( Q e f f / R T ) ) = 0 E13

which lead to the relationship:

Q e f f = n R T N o s e 2 Δ T E14

where TNose is the temperature at the nose position. Eq. (13) can be calibrated by using experimental data taken from literature and getting good estimations of the optimal value of Qeff. It was found in Victor et al. (1998) that the Qeff values have a median value of 27500 kcal/mol °C. For simplicity, this value is used for each diffusion-controlled phase transformation. The kinetic coefficients of alloying elements in Eq. (10) are then determined by calibrating such equation with TTT diagrams in the open literature. Under isothermal conditions, the ferrite transformation can be represented by:

τ F = exp ( 1.00 + 6.31 C + 1.78 M n + 0.31 S i + 1.12 N i + 2.70 C r + 4.06 M o ) 2 0.41 N ( A e 3 T ) 3 exp ( 27500 / R T ) S ( X ) E15

the pearlite transformation is represented by

τ P = exp ( 4.25 + 4.12 C + 4.36 M n + 0.44 S i + 1.71 N i + 3.33 C r + 5.19 M o ) 2 0.32 N ( A e 1 T ) 3 exp ( 27500 / R T ) S ( X ) E16

and the bainite transformation under isothermal condition is represented by

τ B = exp ( 10.23 + 10.18 C + 0.85 M n + 0.55 N i + 0.90 C r + 0.36 M o ) 2 0.29 N ( B S T ) 2 exp ( 27500 / R T ) S ( X ) E17


B s ( ° C ) = 637 58 C 35 M n 15 N i 34 C r 41 M o E18

While, the martensite start temperature can be expressed by the following equation:

M s ( ° C ) = 539 423 C 30.4 M n 17.7 N i 12.1 C r 7.5 M o + 10 C o 7.5 S i E19

Once the TTT diagram is calculated, it is possible to transform it into a CCT diagram using the well–established additivity rule:

0 t d t τ T T T ( X , T ( t ) ) = 0 t Δ T n exp ( Q / R T ) F ( C , M n , S i , N i , C r , M o , G ) S ( X ) d t = 1 E20

where τTTT(X,T(t)) represents the isothermal transformation time for X at temperature T, and t is the total non-isothermal transformation time. Fig. 7 shows an example of TTT diagram computed with Matlab, starting from Eqs. (14-18).

Figure 7.

Calculated TTT diagrams; steel composition (wt%): C 0.37, Mn 0.77, Si 0.15, Ni 0.04, Cr 0.98, Mo 0.21; ASTM grain size number: 7; A= Austenite, F = Ferrite, P = Pearlite, B = Bainite, Ms = martensite start temperature.

The above described model is limited to carbon and low-alloy steels. However, efforts are made in literature in order to develop models for general steels, including medium to high alloy types, tool steels, 13%Cr steels etc.

2.4. Practical use of thermo-metallurgical models for welding

An efficient recognition of the thermal field induced by a welding process, offers the possibility to make a parametric study of welding in order to check the influence of process, material and geometrical parameters on the temperature distribution within the plates. Fig. 8 shows an example of influence of the welding speed on the temperature distribution during the welding. It is clear that the lower the velocity the wider the HAZ.

Figure 8.

Isotherms [°C]: at I = 170 A, V = 23.4 Volt, H = 6.6 mm, T0 = 0 °C and a) v = 5 mm/s b) v = 11 mm/s [Material: AA-5083-O, welding technology: GMAW]

By comparing the solidus temperature isotherm width with that of FZ, it is also possible to estimate the thermal efficiency of the welding process. Finally, by using the CCT diagram and Eq. (6), it is possible to evaluate the pre-heating temperature needed to avoid a martensitic microstructure in the weld bead.


3. Mechanical analysis

3.1. Residual stresses calculation on butt-welded joints

High thermal gradients occurred during a welding process make thermal stresses in welded plates, that develop, after cooling, a state of permanent stress, generally defined as residual stresses. The magnitude of this residual stress state is so great that the mechanical behaviour of the welded joint, and in particular: fatigue, fracture, instability strength and stress corrosion, may be compromised. Thus the knowledge and, above all, the prediction of residual stresses is very important for a correct choice of parameters of the welding process and a good planning of welded joints. Unfortunately, the evaluation of thermal and residual stresses is not an easy task because of the complexity of the phenomena. However, in the case of simple shapes like butt or edge welded plates, some simplifications are possible. In such instances, one can calculate the residual stress field induced by a welding operation by using simple equations which can be solved with an iterative procedure with high time efficiency. The use of such analytical models were already been proposed by several authors. In particular, Goff (1979) simplified the problem assuming temperature-independent materials properties and a linear temperature distribution in the transversal direction of the plate. Using the singular Rosenthal solution of the thermal field, Tall (1964) suggested a step by step trial and error method in which at any temperature increment one had to determinate the equilibrium of thermal stress seen as the summation of temperature and equilibrium stresses; thermal stresses due to the temperature increment were summed to stresses calculate at the previous time step. Agapakis and Masubuchi (1984) developed the previous work by Tall, solving the equilibrium and stress-strain consititutive relations by an iterative procedure. Finally, Cañas et al (1996), proposed a ‘rod model’ in which the previous equations were written in matrix form and advantageously used in an efficient algorithm.

In what follows a model for residual stresses calcualtion in butt-welded joints is described. The equations which solve the problem are written in matrix form like in Cañas‘s work (1996). The thermal field induced by welding is given by Eq. (4) in which some fundamental assumptions are made:

  1. the plate is infinitely large and very thin;

  2. Eq. (4) describes a line source and then no temperature gradient exists through the thickness of the plate;

  3. steady-state conditions;

  4. for the stress calculation it is assumed that at time t, each longitudinal section is a part of an infinitely long plate subject to the same temperature distribution over its entire length (T=T(x,t)) (Fig. 9).

Moreover, plates without lateral constraints are considered so that longitudinal stresses are much greater than the transversal ones which are for this reason neglected. In order to write the equations in matrix form the welded plates are divided into n bars as shown in Fig. 9.

Figure 9.

Schematic representation of welded plates (plate thickness: h)

According to the configuration shown in Fig. 9, the conditions of equilibrium of forces and moments, written in matricial form, are:

C T N y = 0 E21

where C and Ny represent the following matrixes:

C = [ 1 d 1 ... ... 1 d i ... ... 1 d n ] N y = [ N y 1 ... N y i ... N y n ] E22

Nyi is the axial force of the bar i and n is the total number of bars employed to represent the welded plates. Assuming an ideal elasto-plastic non-holonomic material behaviour the constitutive equations are:

q y = q y e + q y t + q y p + Δ q y p E23


  • qy represents the elongation vector in y direction (qi y is the elongation of the bar i in y direction);

  • q e i y = L E ( T i ) h ( N y i b )
  • q t i y = L α ( T i ) ( T i T 0 )
  • qpi y is the accumulated inelastic elongation of the bar i during the previous time increments;

  • Δqpi y is the change of inelastic elongation of the bar i in y direction during the current time increment;

The compatibility equations are:

q y = C u E24

where u is the displacement vector associated to the degrees of freedom δ and θ to which the equilibrium conditions are applied (fig. 10).

Figure 10.

Scheme of welded general plate deformation.

3.1.1. Field equations reduction

Eq. (22) can be written as:

q y = A N y + A t ( T T 0 ) + q p y + Δ q p y E25

where A is a diagonal matrix, each element of it representing the flexibility coefficient of the bar i (L/(E(Ti)hb), with E = Young’s modulus), At is another diagonal matrix where each element of it is Lα(Ti) and T is the temperature vector.

Now by using Eqs. (22) and (23), the vector Ny turns out to be:

N y = K ( C u A t ( T T 0 ) q p y Δ q p y ) E26

where K = [A]-1. From Eqs. (24) and (20):

u = [ C T K C ] 1 C T K ( A t ( T T 0 ) + q y p + Δ q p y ) E27

Now, the only unknown terms are the current plastic elongations (Δqpi y) but a convergence procedure may be used to calculate them. At each time step, the temperature vector (T-T0) can be first calculated by Eq. (4) and thus the temperature dependent material characteristics (E(Ti) etc). The matrixes, K, A and At (qp y is known from the previous time step) may be then determined. Initially assuming that no plastic elongation exists, Eqs. (24) and (25) can be used for a first approximation of Ny. Imposing

| N i y | N p i for i  = 1 n E28

where Npi is the yield force of the bar i, a first approximation of Δqp y can be obtained by using the previous value of u and Eqs. (22) and (23). In this way a second approximation of Ny can be obtained by means of Eqs. (24), (25) and (26) and this procedure can be repeated until convergence is reached for the current time step. A program for the automatic solution of such iterative procedure can be easily written in Matlab code. The input data requested are: the net energy (Q) of the heat source, the welding speed (v), the liquidus and reference temperature, plates dimensions (L, 2B, h), the total number of bars (n). Figure 11 shows an example of the computed residual stresses in a butt welding joint and the comparison with Finite Element and experimental results.

Figure 11.

Residual stresses along the centre line transverse to welding direction: analytical solution (Material: AA-5083-O. Geometrical parameters: h = 0.66 cm, 2B = 36 cm, L = 25 cm. Welding procedure: GMAW; Welding voltage: 23.4 V; Welding current: 170 A; Arc efficiency: 0.64; Welding speed: 11 mm s-1; Filler wire diameter: 1.2 mm; Number of passes: 1; Shielding gas: Argon; Shielding gas flow rate: 0.21 s-1)

Figure 12.

Influence of welding speed (v) on residual stresses

The analytical model is very cost-effective in its computer implementation. Therefore, it allows for a series of parametric analyses to be performed and the relative importance of various parameters to be easily investigated. Fig. 12 shows, for example, the influence of welding speed on the residual stress field in welded plate without lateral constraints. Because of the heat input reduction, a decrement of plastic zone when welding speed increases is found.

3.2. Asymptotic stress distributions in welded and bonded joints

It is well know that fatigue resistance in mecahical components is controlled above all by the singular stress fields which arise near the gemetric discontinuities such as the wleld toe in welded joints (Livieri at al. (2005)) (Fig. 13) or the interface between the substrate and the adhesive in bonded joints (Lazzarin et al. (2002)).

Figure 13.

Welded joint geometry and weld toe.

In several works (Lazzarin at al., 1998; Livieri et al, 2005) the weld toe region is modelled as a sharp, zero radius, V-shaped notch and the intensity of asymptotic stress distributions obeying Williams’ solution (Williams, 1952) are quantified by means of the Notch Stress Intensity Factors (NSIFs). When the constancy of the angle included between weld flanks and main plates is assured and the angle is large enough to make mode II contribution non-singular, mode I NSIF can be directly used to summarised the fatigue strength of welded joints having very different geometry (Livieri et al, 2005). Furthermore, the NSIFs parameters can be used also for the evaluation of thermal fatigue resistance of different components (Ferro et al., 2006; Ferro et al., 2009).

The same approach is used for the advanced planning of bonded joints (Lazzarin et al., 2002). Adhesively bonded joints inevitably present high stress concentration zones, due to the different elastic properties of the connected materials. While in a homogeneous material, linear elastic stress distributions need the presence of a V-shaped corner to become singular (Williams, 1952), in the bi-material problems stress singularity arises, as well known, also in absence of any geometrical discontinuity (Bogy, 1968).

Strength evaluation needs both the order and the intensity of the stress singularity to be quantified in terms of joint geometry, material elastic properties and applied load. Several researchers (Gradin, (1982); Adams, et. al (1987); Groth, (1988); Hattory et al., (1988); Hattory, (1991); Reedy, (1990)) have used H (or other symbols, as K or Q) as a “generalised” stress intensity factor suitable as a failure criterion for bonded joints made of dissimilar materials. In the work of Hattory et al. (1991), for example, the stress field parameter (denoted Qxy) was determined on the basis of the singular distribution of the shear stress σxy present at the interface of Large Scale Integrated (LSI) electronic circuit devices subjected to thermal stresses. The critical value of Qxy (able to provoke delamination in the components) was plotted against the order of the singularity γ. Different values of γ were obtained by using various configurations of epoxy/Fe-Ni blocks bonded together. Previously, Qxy had already been used by Gradin (1982), by introducing his static criterion for brittle edge-bonded bi-material bodies.

More recently, a H-based approach has also been used by Lefrebvre and Dillard (1999) to predict the fatigue crack initiation in epoxy-aluminium wedge specimens, in a manner similar to the use of the Notch Stress Intensity Factors in welded structures (Lazzarin at al., 1998).

In view of the use of stress fields and stress intensity factors to predict the fatigue life at crack initiation, it is important to have the complete and correct description of the stress field very near to the apex. A method for the evaluation of the singular stress field in bonded joints of different geometry is presented and solved with Matlab; the stress distributions are represented by a two terms stress expansion, under the hypothesis that both first and second terms are in the variable separable form and therefore each term can be represented by a radial component with unknown exponent (eigenvalue) and an angular function also unknown.

The resulting analytical formulation of the stress distributions can be given as:

σ i j ( r , ϑ ) = H 0 r s f i j ( 0 ) ( ϑ ) + H 1 r t f i j ( 1 ) ( ϑ ) E29

The method is based on the numerical solution of the ordinary differential equation (ODE) system that, under the hypothesis of plane strain state, derives from the equilibrium and compatibility equations of a bonded joint or, more generally, of a bi-material body.

The capability of the formulation to account for the actual elastic properties of the substrates, allows us to obtain the accurate description of the stress field even in the case of joints made of materials with comparable elastic properties. It is worth noting that the “stress function approach” (where the formulation is completed analytically and the resulting set of equations is solved numerically) has already been used successfully by many researchers, mainly engaged with in-plane crack and notch problems in materials obeying a power-hardening law (Lazzarin et al., 2001).

3.2.1. Analytical frame

Let us consider the problem of the elastic equilibrium in a bi-material joint, in presence of a V-shaped corner with an opening angle ( ϑ 1 + ϑ 2 ) as shown in Figure 14. Both materials are thought of as homogeneous and isotropic and subjected to plane strain conditions. Under linear elastic hypothesis, strains components are:

ε i j = 1 + ν m E m σ i j ν m E m σ k k δ i j E30

where subscript m =1, 2 denotes the material, δ i j is the Kroneker delta, and summation convention is used for repeated indexes.

In writing the problem of the elastic equilibrium we can now consider separately the two materials and later find the solution by applying the boundary conditions at the traction free surfaces and at the interface. It is possible therefore to omit, from now on and until not differently evidenced, the material subscript m, being the equations valid for the both substrates.

Figure 14.

Schematic view of the singular zone showing the cartesian and polar coordinate systems.

By assuming a polar coordinate system, in absence of body forces, the equilibrium conditions can be written as:

σ r r r + 1 r σ r ϑ ϑ + σ r r σ ϑ ϑ r = 0 E31
1 r σ ϑ ϑ ϑ + σ r ϑ r + 2 r σ r ϑ = 0 E32

The compatibility equations between strains and displacements are:

ε r r = U r r E33
ε ϑ ϑ = U r r + 1 r U ϑ ϑ E34
ε r ϑ = 1 2 ( 1 r U r ϑ + U ϑ r U ϑ r ) E35

where U r and U ϑ are the displacement components.

According to the direct approach, first suggested by Ponte Castanêda (1985) and then used also in (Yuan et al., (1994); Lazzarin et al. (2001)), a variable separable two term expansion is used for the stresses:

σ i j ( r , ϑ ) = r s f i j ( 0 ) ( ϑ ) + r t f i j ( 1 ) ( ϑ ) E36

where the exponent s is to be thought of as negative to give a stress field singular and, moreover, it is stated for hypothesis that s < t.

It should be noted that in Eq. (34) the generalised stress intensity factors H0 and H1, related to the first and second order component of the stress distribution respectively, are, for the time being, included in the angular stress distribution functions f i j which are always defined within a constant value. Such a value is to be later determined by means of FE analyses.

Substitution of Eq. (34) into (28) gives the following expressions for strains:

ε i j ( r , ϑ ) = r s ε i j ( 0 ) ( ϑ ) + r t ε i j ( 1 ) ( ϑ ) E37


ε i j ( 0 ) ( ϑ ) = 1 + ν E f i j ( 0 ) ( ϑ ) ν E f k k ( 0 ) ( ϑ ) δ i j E38
ε i j ( 1 ) ( ϑ ) = 1 + ν E f i j ( 1 ) ( ϑ ) ν E f k k ( 1 ) ( ϑ ) δ i j E39

The relevant displacement components are:

U i ( r , ϑ ) = r s + 1 U i ( 0 ) ( ϑ ) + r t + 1 U i ( 1 ) ( ϑ ) E40

As it is known, the exponents depend on the combination of the material elastic properties.

In the close neighbourhood of the singularity point (r which tends towards zero), the first term of the stress distribution becomes dominant. Let us consider, therefore, in the stress and strain expansions only the leading-order term, s being the relevant exponent.

Substitution of Eqs. (34), (35) and (38) into Eqs. (29) to (33), together with the plane strain conditions

f z ( 0 ) ( ϑ ) = ν [ f r r ( 0 ) ( ϑ ) + f ϑ ϑ ( 0 ) ( ϑ ) ] E41

gives the following system:

( s + 1 ) f r r ( 0 ) ( ϑ ) + f r ϑ , ϑ ( 0 ) ( ϑ ) f ϑ ϑ ( 0 ) ( ϑ ) = 0 f ϑ ϑ , ϑ ( 0 ) ( ϑ ) + ( s + 2 ) f r ϑ ( 0 ) ( ϑ ) = 0 ( s + 1 ) U r ( 0 ) ( ϑ ) 1 + ν E f r r ( 0 ) ( ϑ ) + ν E [ f r r ( 0 ) ( ϑ ) ( 1 + ν ) + f ϑ ϑ ( 0 ) ( ϑ ) ( 1 + ν ) ] = 0 U r ( 0 ) ( ϑ ) + U ϑ , ϑ ( 0 ) ( ϑ ) 1 + ν E f ϑ ϑ ( 0 ) ( ϑ ) + ν E [ f r r ( 0 ) ( ϑ ) ( 1 + ν ) + f ϑ ϑ ( 0 ) ( ϑ ) ( 1 + ν ) ] = 0 1 2 [ U r , ϑ ( 0 ) ( ϑ ) + s U ϑ ( 0 ) ( ϑ ) ] 1 + ν E f r ϑ ( 0 ) ( ϑ ) = 0 E42

where f , ϑ and U , ϑ mean f ϑ and U ϑ , respectively.

The boundary conditions, at the traction free surface and at the interface between the two materials, are as follows:

f ϑ ϑ m a t .1 ( γ 1 ) = f ϑ ϑ m a t .2 ( + γ 2 ) = 0 f r ϑ m a t .1 ( γ 1 ) = f r ϑ m a t .2 ( + γ 2 ) = 0 f ϑ ϑ m a t .1 ( 0 ) = f ϑ ϑ m a t .2 ( 0 ) f r ϑ m a t .1 ( 0 ) = f r ϑ m a t .2 ( 0 ) U ϑ m a t .1 ( 0 ) = U ϑ m a t .2 ( 0 ) U r m a t .1 ( 0 ) = U r m a t .2 ( 0 ) E43

By using the third equation of the system 40 to eliminate the U r ( 0 ) ( ϑ ) , the problem gives four differential equations for the material 1 and four differential equations for the material 2. The angular stress and displacement components f ϑ ϑ , f r ϑ , f r r and U ϑ are the eigenfunctions of the problem, which will be given to within a constant value, as the problem is homogeneous.

A PC-based application has been developed for the solution of the system, by using the MATLAB® program and, particularly, its ODE45® routine.

According to the procedure illustrated in the flowchart of Figure 15, the solution begins with the introduction of an arbitrary value for U ϑ ( γ 1 ) (for example, 104÷106) together with a couple of s and f r r ( γ 1 ) guess values and proceeds in iterative form till up:

a b s { f ϑ ϑ m a t .2 ( γ 2 ) } + a b s { f r ϑ m a t .2 ( γ 2 ) } t o l e r a n c e E44

the tolerance value being usually set equal to 1015, while the angular stress components are close to the unity.

Figure 15.

Procedure for the solution of the ODE system (40).

To obtain the solution, the MATLAB® FMIN® function can, automatically, slightly modify the guessed values and restart the procedure. In cases where the convergence is not reached within a certain number of iterations, substantially different guessed values have to be chosen and the solution to be restarted.

As the stress expansion is in a variable separable form, and due also to the hypothesis of linear elastic behaviour, if we consider, in the stress and strain expansions, only the second term and the associated t order exponent, the system solution provides the related eigenfunctions.

Once the solution of the ODE system (40) is completed, the eigenfunctions, that is the angular stress and displacement distribution functions, and also the singularity strength for both the leading and second order term and for material 1 and 2 are available.

For the complete description of the stress field, according to Eq. (27), the generalised stress intensity factors H0 and H1 have still to be determined and this can be done through FE analyses, by minimising, for a generic ϑ * direction (arbitrarily chosen on the steel side), the following error function:

E . F . ( H 0 , H 1 ) = r , min r , max a b s . [ H 0 r s f i j ( 0 ) ( ϑ * ) + H 1 r t f i j ( 1 ) ( ϑ * ) σ i j ( F E ) ] E45

being the rmin and rmax values usually set to 104 mm and one tenth of the substrate thickness, respectively. It is important to note that, kept constant the mesh refinement, the adhesive and substrate thickness varied from model to model, being 0.5-1.0 mm and 2.0-4.0 mm, respectively, their more typical value ranges. Fig. 16 shows the asymptotic stress distribution near the interface between the adhesive and the substrate in bonded joints obtained with Matlab.

Figure 16.

Asymptotic stress distribution near the interface between the adhesive and the substrate in bonded joints


4. Conclusion

Different thermal and mechanical models solved with Matlab were presented which are very useful for a good planning of welded or bonded joints. About the thermal and residual stress fields predictions, some limitations are due to the geometry and material properties, but, compared to the finite element models, they are more user friendly and more efficient in terms of computational time; thus they can be used both for a rapid check of the thermal and stress filed induced by welding and for a parametric study of the process.

The ‘stress field approach’ for fatigue resistance prediction of welded and bonded joints was also presented. In particular, an analytical method for the description of the singular stress distributions on bonded joints of different geometry has been developed. The stress distributions near the singularity have been assumed to be represented by a power series expansion in the variable separable form as usually done in the elastoplastic analyses; under the further hypothesis of plane strain state, applying the equilibrium and compatibility conditions results in a ordinary differential equation system, which has been numerically solved by using a “shooting” technique.


  1. 1. Rosenthal D. 1941Mathematical theory of heat distribution during welding and cutting, Weld. J. Res. Supp., 20, 220s.
  2. 2. Rosenthal D. Shamerber H. 1938Thermal study of Arc Welding. Experimental Verification of Theoretical Formulas. The Welding Journal Res. Suppl., 17 4 2 8
  3. 3. Ferro P. Tiziani A. Bonollo F. 2002Analisi numerica e teorica del campo di tensioni residue indotte dal processo di saldatura in piastre in lega leggera. Proceedings of XXXI AIAS Conference, Parma, Italy
  4. 4. Ferro P. Bonollo F. 2010A Semianalytical Thermal Model for Friction Stir Welding. Metallurgical and Material transaction A, 41A 440 449
  5. 5. Vilaca P. Quintino L. dos Santos. J. F. Zettler R. Sheikhi S. 2007Quality assessment of friction stir welding joints via an analytical thermal model, iSTIR’. Materials Science and Engineering A, 445-446 501 508
  6. 6. Victor Li. M. Niebuhr V. Meekisho L. Atteridge G. 1998A Computational Model for the Prediction of Steel Hardenability. Metallurgical and Materials Transactions B, 29B 661 672
  7. 7. Schmidt H. B. Hattel J. H. 2008Thermal modelling of friction stir welding’. Scripta Materialia, 58 332 337
  8. 8. Perivilli S. Peddieson J. Cui J. 2008Simplified Two-Dimensional Analytical Model for Friction Stir Welding heat Transfer. Journal of Heat Transfer, 130 1 9
  9. 9. Goff R. F. D. Porter 1979A Simplified Analysis of the Residual Longitudinal Stresses and Strains Due to the Gas-Cutting and Welding of Thin Steel Plate. Int. J. Mech. Sci. 21 287 300
  10. 10. Chao Y. J. Qi X. Tang W. 2003Heat Transfer in Friction Stir welding- experimental and Numerical Studies. Transaction of ASME, 125 138 145
  11. 11. Tall L. 1964Residual Stresses in Welded plates- a Theoretical Study. Weld. J. Res. Supp. 43s.
  12. 12. Agapakis J. E. Masubuchi K. 1984Analytical Modeling of Thermal Stress Relieving in Stainless and High Strength Steel Weldments”, Weld J. Res. Suopp. 36, 187s.
  13. 13. Cañas J. Picòn R. París F. Blazquez A. Marín J. C. 1996A Simplified Numerical Analysis of Residual Stresses in Aluminium Welded Plates, Computers & Structures, 58 1 59 69
  14. 14. Cañas J. Picòn R. París F. Del Río J. I. 1996A One-Dimensional Model for the Prediction of Residual Stress and its Relief in Welded Plates, Int. J. Mech. Sci. 38 735 751
  15. 15. Lazzarin P. Quaresimin M. Ferro P. 2002A two terms stress function approach to evaluate stress distributions in bonded joints of different geometry, Journal of Strain Analysis for Engineering Design, 37 5 385 398
  16. 16. Livieri P. Lazzarin P. 2005Fatigue strength of steel and aluminium welded joints based on generalised stress intensity factors and local energy values, International Journal of Fracture, 133 247 276
  17. 17. Williams M. L. 1952Stress singularities resulting from various boundary conditions in angular corners of plates in exstension. Journal of Applied Mechanics, 19 526 528
  18. 18. Bogy D. B. 1968Edge-bonded dissimilar orthogonal elastic wedges under normal and shear loading. J. Appl. Mechanics, 35 460 466
  19. 19. Gradin P. A. 1982A fracture criterion for edge-bonded bi-material bodies. J. Composite Mater., 16 448 456
  20. 20. Adams R. D. Harris J. A. 1987The influence of local geometry on the strength of adhesive joints. Int. J. Adhes. Adhes., 7 69 80
  21. 21. Groth H. L. 1988Stress singularities and fracture at interface corners in bonded joints, Int. J. Adhes. Adhes.,, 8 107 113
  22. 22. Hattory T. Sakata S. Hatsuda T. Murakami G. 1988A stress singularity parameter approach for evaluating adhesives strength. JSME Int. J., Series I, 31 718 723
  23. 23. Hattory T. 1991A stress singularity parameter approach for evaluating the adhesives strength of single lap joint. JSME. Int. J., Series I, 34 326 331
  24. 24. Reedy Jr E. D. 1990Intensity of stress singularity at the interface-corner between a bonded elastic and rigid layer. Engng Fracture Mechanics, 36 575 583
  25. 25. Lefebvre D. R. Dillard D. A. 1999A stress singularity approach for the prediction of fatigue crack initiation in adhesive bonds. PART 1: Theory. J. Adhesion, 70 119 138
  26. 26. Lazzarin P. Tovo R. 1998A Notch Stress Intensity Factor Approach to the Stress Analysis of Welds. Fatigue Fracture Engng Mater. Structs, 21 1089 1104
  27. 27. Lazzarin P. Zambardi R. Livieri P. 2001Plastic Notch Stress Intensity Factors for Large V-Shaped Notches under Mixed Load Conditions. Int. J. Fracture, 107 361 377
  28. 28. Ponte Castañeda. P. 1985Asymptotic fields in steady crack growth with linear strain-hardening. J. Mech. Phys. Solids, 35 227 268
  29. 29. Yuan H. Lin G. 1994Analysis of elastoplastic sharp notches. Int. J. Fracture,, 67 187 216
  30. 30. Ferro P. Petrone N. 2009Asymptotic Thermal and Residual Stress Distributions due to Transient thermal Loads. Fatigue and Fracture of Engineering Materials and Structures. 32 936 948
  31. 31. Ferro P. Berto F. Lazzarin P. 2006Generalized stress intensity factors due to steady and transient thermal loads with applications to welded joints. Fatigue Fract. Engng. Mater. Struct, 29 440 453

Written By

Paolo Ferro

Submitted: November 15th, 2010 Published: September 9th, 2011