Analytical solution of thermal stress in a bi-metallic cylinder

As a side project I have looked to solve using Python an analytical solution for one-dimensional transient heat transfer in a thick cylinder with cladding on the inside surface, wherein the fluid in contact with the inner surface (cladding) undergoes an arbitrary temperature variation. The solution consists of estimating the time-varying temperature distribution in the thickness of the cylinder and then determine the elastic thermal stresses from knowledge of the temperature distribution. The cylinder in question has an internal radius of 2000mm and external radius of 2207.5 mm. The cladding thickness on the inner surface is 7.5 mm, with the remainder of the cylinder being of 200mm thickness. The cladding possesses different thermal and mechanical properties to the bulk of the cylinder thickness. I have looked to verify the accuracy of the analytical solution using Abaqus FE software.

Solving analytically for the temperature distribution in both the cladding and the cylinder is relatively trivial and matches closely the numerical solution derived from FE analysis 1.

1.Transient Temperature Solution

Hooke’s law is used to determine the elastic through-wall stresses at any time considering the distribution of temperature through-thickness. For instance, the axial stress at each radial position, r, in the thickness of the cylinder is calculated according to the 2, where rv is the internal radius, re is the external radius, and the integral is usually solved numerically using trapezoidal/simpson’s rule. It is made a bit more difficult than the case of a single material since the Young’s modulus (E) and coefficient of thermal expansion (cte) vary spatially as a function of the radial position; therefore it cannot be assumed constant and must be considered in the integration. At the boundary between the cladding and cylinder, there is a step change in the material properties due to the material discontinuity.

2.Axial Stress, Hooke’s law

Considering the case of t=3000s during the thermal transient. Assuming I want to analytically solve for the stress distribution at 22 points through-thickness, 11 in the cladding and 11 in the cylinder, I have the following arrays which define the radial position (r), Young’s modulus (E), thermal expansion coefficient (cte), and temperature distribution through-wall (temp) at this instant in time (see example code below). However, when solving for the axial stress at this time, the through-wall stress distribution looks reasonable for the bulk cylinder but is clearly incorrect in the cladding and does not resemble the numerical solution obtained from FE analysis 3. I have tried using a much finer discretisation through-thickness (400 points) and still have the same issue. I suspect it may be due to the discontinuity in the integrand.

3.Incorrect through-wall stress

Any suggestions as to what may be causing the error here and how I might resolve it?

Many thanks for your time.

def Ir(r,cte,T):
    return cte*T*r

from numpy import trapz, array

# radial position (internal radius = 2000, external radius = 2207.5, interface = 2007.5 all in mm)
r = array([2000,2000.75,2001.5,2002.25,2003,2003.75,2004.5,2005.25,2006,2006.75,2007.5,
           2007.5,2027.5,2047.5,2067.5,2087.5,2107.5,2127.5,2147.5,2167.5,2187.5,2207.5],dtype=float)

# Young's modulus of cladding and cylinder (MPa)
E  = array([197000,197000,197000,197000,197000,197000,197000,197000,197000,197000,197000,
           204000,204000,204000,204000,204000,204000,204000,204000,204000,204000,204000],dtype=float)

# Coefficient of thermal expansion of cladding and cylinder (degC^-1)
cte = array([1.64e-5,1.64e-5,1.64e-5,1.64e-5,1.64e-5,1.64e-5,1.64e-5,1.64e-5,1.64e-5,1.64e-5,1.64e-5,
1.09e-5,1.09e-5,1.09e-5,1.09e-5,1.09e-5,1.09e-5,1.09e-5,1.09e-5,1.09e-5,1.09e-5,1.09e-5],dtype=float)

# Temperature in cladding and cylinder (degC)
temp = array([23.13,24.78,26.445,28.11,29.76,31.43,33.09,34.75,36.41,38.07,39.73,
              39.73,57.31,73.91,89.21,102.97,114.97,125.03,133.0,138.77,142.26,143.43],dtype=float)

# Poisson's ratio (assumed to be 0.3 for both cladding and cylinder)
v = 0.3

# Solve numerically for the integral
I1 = trapz(Ir(r,cte,temp),x=r)

ax_stress = ((E)/(1.0-v))*((2.0/(r[-1]**2-r[0]**2))*I1-temp*cte)