Elastic–plastic analysis of the plane strain under combined thermal and pressure loads with a new technique in the finite element method

Akeel Ali Wannas; Auday Shaker Hadi; Noor H. Hamza

doi:10.1515/eng-2022-0049

Article Open Access

Elastic–plastic analysis of the plane strain under combined thermal and pressure loads with a new technique in the finite element method

Akeel Ali Wannas , Auday Shaker Hadi and Noor H. Hamza

Published/Copyright: July 22, 2022

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information

From the journal Open Engineering Volume 12 Issue 1

Abstract

This research aims to find a new way to address thermal loads within the framework of the elastic–plastic relationship, especially when the loads, such as thermal and pressure ones, are combined. While the residual force method was employed to attempt to find the convergence of a nonlinear solution, it was unable to do so. For this purpose, a mathematical relationship was derived to address thermal loads and add them to the hierarchy of the specified elements method as a subroutine. The findings of the developed program were verified by comparing the numerical results with those of the analytical solution of a thick-walled cylinder loaded with heat load and internal pressure; the results proved the correctness and accuracy of the method used. The method offers a way to redevelop old programs that are unable to effectively address the thermal load in elastic–plastic relations without changing said programs significantly but only by adding the subroutine and some very simple modifications. The solution technique provided in this article can be utilized in many related cases, such as plane stress, axisymmetric solid, and three-dimensional stress analyses.

Keywords: thermal load; residual force; thick cylinder; elastic–plastic analysis; plane strain

1 Introduction

Most practical engineering problems are difficult to solve theoretically. And while experimental methods can give accurate results, they are often expensive and time-consuming. The problem becomes more complicated with loads such as pressure, point load, and gravity with temperature effects (which are usually accompanied by initial stress and strain, alteration of material properties, and creep). In general, the literature on the finite element method in the solution and analysis of elastic–plastic relationships can be accessed in textbooks. However, current concerns are focused on improving the performance of the software by adding variables to the constitutive law [1]. The more the variables, the more complex the mathematical and numerical solutions.

With the advent of computers in the 1960s, the finite element method gained real impetus [2], especially within the study of elastic relationships of materials. Many books solving the linear state, including thermal loads, have been written on this subject [3]; on the topic of the residual force method for solving nonlinear problems, such as elastic–plastic and elastic-viscoelastic analysis, Owen and his co-authors [4,5,6] as well as others [7] have written several books and conducted much research.

The residual force method was initially used to solve nonlinear stress analysis in the elastic–plastic relationships by the finite element method. However, this method was ineffective in solving nonlinear relationships with the presence of thermal loads, especially when combined with other loads, such as pressure ones. We note that in their book, Hinton and his co-authors [3] dealt with the programming of the finite element method in linear relations and took into account the thermal loads in it. But in their next book, published in 1980 [6] and based on the previous one – which dealt with the programming of the finite element method in plasticity – they avoided mentioning the thermal loads, suggesting that no solution for this problem had been found until that book. During that period, researchers also began developing the Constitutive Law and its programming in the finite element method, which in turn solved most of the non-linear problems – if not all of them – including thermal loads.

The main questions for this study are: Why does the residual force method not work in nonlinear relationships with thermal loads? And if we can answer this question, how can this problem be solved?

2 Thermal stress in elasticity

In thick-walled cylinders, when there is a heat source either inside or outside, it will cause the heat to flow to the other side. This will result in a temperature difference across the wall of the cylinder, meaning that each layer will expand differently from the one next to it due to the presence of thermal expansion. If we assume isotropic material with a thermal expansion α coefficient, then the thermal strains are given by ref. [8]

(1) { ε t } = { α } Δ T ,

where Δ T is a temperature change between assumed layers and ε t is the thermal strain.

Now, it is possible to find the nodal force on an element { f t } due to thermal stress, if we assume free expansion is inhibited

(2) { f t } = ∫ Ω [ B ] T [ D ] { ε t } dΩ ,

where [ B ] T is the strain-displacement matrix, [ D ] is the stress-strain matrix, and d Ω is the domain of interest. Now, we can combine the force matrix of one element with all the others to form a global force matrix.

To calculate the stresses{σ} in a thermally loaded system, the thermal strains {ϵ ^t} must be subtracted from the sum of the strains {ϵ} after solving the set of arrays and extracting the nodal displacements { u }

(3) { ϵ } = [ B ] { u } ,

(4) { σ } = [ D ] ( { ϵ } − { ϵ t } ) ,

(5) { σ } = ( { σ ′ } − { σ t } ) ,

where { σ ′ } represents the stress caused by the thermal force, and { σ t } the stress caused by the thermal strain.

It is possible to find the relationship by calculating the specific force that caused the true stress

(6) { f } = ∫ Ω [ B ] T ( { σ ′ } − { σ t } ) dΩ .

This indicates that there are two parts of thermal expansion: one that expands without causing thermal stress, and another that prevents expansion, turning it into thermal stresses.

3 Numerical solution processes for the residual force method

Here, we produce a system formula of large class nonlinear simultaneous equations when solving the elastic–plastic formula by the finite element method

(7) H φ + f = 0 ,

where φ is the basic unknowns (displacements) vector, f is the applied ‘loads’ vector, and H is the assembled “stiffness” matrix.

If the elements of the matrix H depend on the values of the unknown vector φ , the problem is likely to be nonlinear. In this case, the direct solution of the equation system is not possible; therefore, the iterative solution mode must be used. The current research focuses on using the residual force method to solve the nonlinear simultaneous equations strain; Armen mentions the advantages of this technique in his research [9].

During any step of the iterative solving process, equation (7) will not be satisfied unless a convergence process occurs. It is possible to assume a residual force system [6], so that

(8) ψ = H φ + f ≠ 0 ,

where the residual force ψ value can be interpreted as the magnitude of moving away from equilibrium.

Since H is a function of φ , at any stage of the iteration ψ = ψ ( φ ) .

If the true solution to the problem exists at φ r + ∆ φ r , the Newton–Raphson approximation for the general terms to the vector of residual force ψ r solvable at φ r is

(9) ψ i r = − ∑ j = 1 N Δ φ j r ∂ ψ i ∂ φ j r ,

Here, the superscript r denotes the number of iterations to the r th approximation of the true solution, and N is the total number of variables. Substituting for ψ i from equation (8), the whole expression for all the residual force components can be written in matrix form as

(10) ψ ( φ r ) = − J ( φ r ) ∆ φ r ,

where a term of the Jacobian matrix J is

(11) J i j = ∂ ψ i ∂ φ j r = h i j r + ∑ k = 1 m ∂ h i k ∂ φ j φ k r ,

and h ij is the public term of matrix H . Nonsymmetrical terms in the Jacobian matrix will appear in the last term of equation (11). If these nonsymmetric terms are neglected to keep up symmetry, then substitution of equation (11) in equation (10) results in

(12) H ( φ r ) · ∆ φ r = − ψ ( φ r ) ,

or since

∆ φ r = φ r + 1 − φ r

Equation (12) can be reduced, when using equation (8), to

(13) H ( φ r ) · φ r + 1 + f = 0 .

It is necessary to keep the asymmetric elements in order to achieve a better convergence associated with the Newton–Raphson method. The explicit form of the nonlinear terminology in equation (11) will depend on the way in which the coefficients of the stiffness matrix h ij depend on the unknown φ. It is possible to assemble the Jacobian matrix terms (11) into a general expression:

(14) J ( φ ) = H ( φ ) + H ′ ( φ ) .

The process of Newton–Raphson can be written using equations (10) and (14) with the following formula:

∆ φ r = − [ J ( φ r ) ] − 1 · ψ ( φ r ) ,

(15) ∆ φ r = − [ H ( φ r ) + H ′ ( φ r ) ] − 1 · ψ ( φ r ) .

The above formula allows us to correct the unknown vector φ , which we can obtain from the residual force vector ψ for any iteration. Therefore, an iterative procedure must be followed with the unknown vector φ , which is corrected at each stage of the iteration according to formula (15) until convergence is achieved.

4 Expressions of the finite element

The basic form of the solution can be obtained by using the principle of virtual work. If we take into account the solid, the external force f , the internal stresses σ , and the distribution of force per unit volume b , all of them form an equilibrium field when we subject an arbitrary virtual displacement pattern δ d * which results in internal displacements δ u * and incompatible strains δ ϵ * . The virtual working principle will require the following formula:

(16) ∫ Ω ( δ ε ⁎ T σ − δ u ⁎ T b ) dΩ − σ d ⁎ T f = 0 .

Then, the normal procedure of discretizing the finite element leads to the following expressions for the strains and displacements within any element

(17) δ u * = N δ d * , δ ϵ * = B δ d * ,

where B and N are the elastic strain matrix and shape functions, respectively. Then, the process of element assembly gives

(18) ∫ Ω δ d ⁎ T ( B T σ − N T b ) dΩ − δ d ⁎ T f = 0 .

Equation (16) will not generally be satisfied at any stage of the calculation, and

(19) ψ = ∫ Ω B T σ dΩ − f + ∫ Ω N T b dΩ = 0 ,

before calculating the residual force, the load matrix must be recalculated from the results of the first iteration (stress, strain, and displacement according to linear relationships) using equation (6). The new load matrix replaces the old one for the rest of the iterations until equilibrium is achieved. As for strains and displacements, it remains the same

(20) f = ∫ Ω B T ( σ ′ − σ t )dΩ = ∫ Ω B T σ dΩ .

The stiffness of a material is constantly changing with the stress/strain relationship for an elastic–plastic situation. To estimate the material tangential stiffness matrix K T at any step [6,10], the form of incremental (19) must be employed. Thus, within an increase of load, we achieve

(21) Δ ψ = ∫ Ω B T Δ σ dΩ − Δ f + ∫ Ω N T Δ b dΩ ,

(22) d σ = D ep d ϵ ,

where D ep is the elastic–plastic matrix (the associated flow rule is employed to find it)

(23) Δ ψ = K T d − Δ f + ∫ Ω N T Δ b dΩ ,

where

(24) K T = ∫ Ω B T D ep B dΩ .

5 Program structure

All the mathematical relationships mentioned in the above paragraphs are found in the literature on the subject, except for relationship (20), which is considered the main focus of our research. With this relationship, a new force matrix was created to replace the original one, which in turn contained unreal values that worked, upon solution, only to increase the displacements (resulting from thermal expansion) without affecting the stress values. This step comes directly after the linear solution in the first iteration, and we preserve the values of the displacements and strain. Then, we can complete the solution as instructed in the literature of this field (as if it was not a thermal load).

The structure of the program in general depends on references [6,11], which were considered when writing the program for our research. Some modifications were introduced to suit the proposed method, such as the one where a subroutine was added to calculate the new load matrix according to equation (20) and replaced it with the original matrix. Figure 1 provides an explanation and description of the main program structure and the function of the sub-programs. The program was written in a programming language (Fortran 90).

Figure 1

Program structure for elastic–plastic applications.

Several terms and symbols appear in Figure 1, the details of which we did not address or define, as this document does not allow for listing such a large amount of information and equations. Therefore, we focus solely on presenting new information in order not to confuse the reader, and also because prior research can be reviewed through references [6,10,12]. The section in the dashed rectangle in Figure 1 represents a new addition to the software architecture of this type, which does not necessarily demand rewriting the existing programs that are unable to handle thermal loads; it suffices to simply add a new subroutine to upgrade them, which is what has been done.

6 Numerical example

Verifying the results of the proposed technique requires comparing them with the results of experiments or those of an analytical solution to the same issue. Here, we decided on the second option due to its ease and lack of cost.

6.1 Analytical solution

To verify the results of the method used to calculate thermal stresses, deformations, and elastic–plastic expansion, we compare the results of the finite element program with those of a theoretically solved state. We follow the example of reference [13] which deals with solving the problem of a thick-walled cylinder loaded with internal pressure and thermal load.

A thick-walled cylinder (plain strain) made of isotropic material has the yield stress Y , the modulus of elasticity E , the inner radius a , the outer radius b , a position of the elastic–plastic expansion c , a temperature of the inner surface T a , and the temperature of the outer surface T b , while the distribution of temperature across the wall of the cylinder T satisfies the Laplace equation ∇ 2 T = 0 , the solution of which may be written as the stress equations

In elastic region c ≤ r ≤ b
(25) σ r = − 1 2 Y − β E ln b a c 2 r 2 − c 2 b 2 − β E ln b r ln b a ,

(26) σ θ = 1 2 Y − β E ln b a c 2 r 2 − c 2 b 2 + β E 1 − ln b r ln b a .
In plastic region a ≤ r ≤ c

(27) σ r = − 1 2 Y − β E ln b a 1 − c 2 b 2 + ln c 2 r 2 − β E ln b r ln b a ,

(28) σ θ = 1 2 Y − β E ln b a 1 + c 2 b 2 − ln c 2 r 2 + β E 1 − ln b r ln b a ,

(29) σ θ − σ r = Y = 2 k ,

where equation (29) is Tresca’s yield criterion for a non-hardening material and Y is the yield stress; it is equal to twice the shear stress k. The shear stress is 0.577 (for mild steel) of the uniaxial yield stress (for mild steel 24 dN/mm²). However, we cannot compare the axial stress with the analytical solution because the initial conditions cannot be recreated; it is not possible to simulate the axial force resulting from the internal pressure. Luckily, this will not impact the remaining results of the stresses, because they are not affected by the axial force.

The equation below can be used to find the pressure required when the elastic–plastic expansion c is assumed:

(30) p = 1 2 Y − β E ln b a 1 − c 2 b 2 + ln c 2 a 2 + β E ,

where β = α ( T a − T b ) 2 ( 1 − ν ) , ν is the Poisson’s ratio, and α the coefficient of thermal expansion.

Table 1 presents the values for cylinder geometry and material properties.

Table 1

The values of cylinder geometry and material properties

a (mm)	b (mm)		E (dN/mm²)
100	200		21,000
Y (dN/mm²)	υ	H ′	α ( m / ( m ° C )
27.696	0.3	0.0	0.00001

6.2 Finite element case

Figure 2 represents a numerical example of a quarter of a cylinder. The geometry data for this example shown in Table 2 were used in the finite element program to calculate the stresses, elastic–plastic expansion, and deformations to be compared with the results of the analytical solution of the same example.

Figure 2

Mesh employed in the elastic–plastic analysis of a thick cylinder subjected to thermal and internal pressure loads under plane strain conditions.

To choose the minimum number of nodes and elements that give the least possible errors, convergence was studied. Here, we obtained fewer numbers than in Table 2 but increased them to achieve as smooth curves of results as possible.

Table 2

Data geometry for the numerical example

No.	Title	The description
1	Number of nodes	341
2	Number of elements	100
3	Element type	The 8-node quadrilateral element
4	Temperature distribution	T = T b + ( T a − T b ) ln ( b / r ) / ln ⁡ ( b / a )
5	Uniaxial yield stress	24.0 dN / ( m m 2 )
6	Yield criteria	Tresca yield criteria
7	Strain hardening parameter H′	0.0
8	Initial condition	All horizontal nodes parallel to the x-axis that start with node 1 and end with node 321 have free motion in the direction of the x-axis only
8	Initial condition	All vertical nodes parallel to the y-axis that start with node 21 and end with node 341 have free motion in the direction of the y-axis only

7 Results and discussion

To ensure that the program results are satisfactory, we must test the program and the method used in the research in many conditions. For example, the temperature change between inside and outside of the cylinder (whether the inside is larger than the outside, or the inside is smaller than the outside; this will be elaborated on later) and the change of pressure and the extent of its influence on the program outputs (such as hoop stress, radial stress, bore strain, and elastic–plastic expansion).

Figures 3 and 4 show the changes in pressure and bore strain, respectively, during the elastic–plastic expansion when considering three cases of the temperature distribution. The first occurs when the inside and outside temperatures are the same; the second when the internal surface temperature is 40°C higher than the outer surface; in the third case, the temperature of the inner surface is 40° lower than that of the outer surface. These results corroborate those of the analytical and numerical solutions.

Figure 3

Variation of internal pressure during the elastic–plastic expansion.

Figure 4

Variation of dimensionless radial displacement during the elastic–plastic expansion.

Figures 5–8 show the results between hoop and radial stress through the cylinder wall with elastic–plastic expansion for three cases where the internal surface temperature is higher or lower than the outer surface by 40°C. Again, the results match those of the analytical and numerical solutions.

Figure 5

Distribution of hoop stresses in an elastic–plastic cylinder subjected to internal pressure when β = 0.25.

$Figure 6 Distribution of radial stresses in an elastic–plastic cylinder subjected to internal pressure when β = 0 . 25 \hspace{0em}\beta =0.\hspace{-0.33em}25 .$

Figure 6

Distribution of radial stresses in an elastic–plastic cylinder subjected to internal pressure when β = 0 . 25 .

$Figure 7 Distribution of hoop stresses in an elastic–plastic cylinder subjected to internal pressure when β = − 0.25 \beta =-0.25 .$

Figure 7

Distribution of hoop stresses in an elastic–plastic cylinder subjected to internal pressure when β = − 0.25 .

$Figure 8 Distribution of radial stresses in an elastic–plastic cylinder subjected to internal pressure when β = − 0.25 \beta =-0.25 .$

Figure 8

Distribution of radial stresses in an elastic–plastic cylinder subjected to internal pressure when β = − 0.25 .

Any changes in the properties of the material due to different temperatures were not considered, as the properties of steel do not alter at a temperature of 40°C; the theoretical equations by which the results were compared are also free of such effects.

8 Conclusion

According to the above results, the method used in this research is effective in finding the values of stress, strain, and elastic–plastic expansion, even with the combination of all types of loads with convection. In solving for convergence, the tangential stiffness matrix method was used; however, it is also possible to employ the initial stiffness matrix or any other method based on the residual forces. One can also use the same method in all cases to calculate stress, strains, and other factors, such as plane stress, axisymmetric solid, or three-dimensional analysis.

Conflict of interest: The authors declare that there is no conflict of interest regarding the publication of this article.

References

[1] Lehmann T. The constitutive law in thermoplasticity an introduction. Vienna: Springer; 1984.10.1007/978-3-7091-2636-3Search in Google Scholar

[2] Hinton E, Irons B. Least squares smoothing of experimental data using finite elements. Strain. 1968;4(3):24–7. 10.1111/j.1475-1305.1968.tb01368.x.Search in Google Scholar

[3] Hinton E, Owen DRJ. Finite element programming, computational mathematics and applications. London: Academic Press; 1977.Search in Google Scholar

[4] Griffiths JR, Owen DRJ. An elastic–plastic stress analysis for a notched bar in plane strain bending. J Mech Phys Solids. 1971;19(6):419–31. 10.1016/0022-5096(71)90009-3.Search in Google Scholar

[5] Owen DRJ, Holbeche J, Zienkiewicz OC. Elastic–plastic analysis of fibre-reinforced materials. Fibre Sci Technol. 1969;1(3):185–207. 10.1016/0015-0568(69)90017-7.Search in Google Scholar

[6] Owen DRJ, Hinton E. Finite in elements in plasticity: theory and practice. UK: Pineridge Press Limited; 1980.Search in Google Scholar

[7] Aggarwal SK, Nayak GC. Elasto-plastic analysis as a basis for design of cylindrical pressure vessels with different end closures. Int J Press Vessel Pip. 1982;10(4):271–96. 10.1016/0308-0161(82)90036-9.Search in Google Scholar

[8] Smith IM, Griffiths DV, Margetts L. Programming the finite element method. UK: John Wiley & Sons Inc.; 2015.10.1002/9781119189237Search in Google Scholar

[9] Armen H. Assumptions, models, and computational methods for plasticity. Comput Struct. 1979;10(1–2):161–74. 10.1016/0045-7949(79)90084-1.Search in Google Scholar

[10] Zienkiewicz OC, Valliappan S, King IP. Elasto‐plastic solutions of engineering problems ‘initial stress’, finite element approach. Int J Numer Methods Eng. 1969;1(1):75–100. 10.1002/nme.1620010107.Search in Google Scholar

[11] Kalali AT, Hassani B, Hadidi-Moud S. Elastic–plastic analysis of pressure vessels and rotating disks made of functionally graded materials using the isogeometric approach. J Theor Appl Mech. 2016;54(1):113–25. 10.15632/jtam-pl.54.1.113.Search in Google Scholar

[12] Yamada Y, Yoshimura N, Sakurai T. Plastic stress-strain matrix and its application for the solution of elastic–plastic problems by the finite element method. Int J Mech Sci. 1968;10(5):343–54. 10.1016/0020-7403(68)90001-5.Search in Google Scholar

[13] Charabarty J. Theroy of plasticity. 3rd ed. Amsterdam: Elsevier Butterworth-Heinemann; 2006.Search in Google Scholar

Received: 2022-04-07

Revised: 2022-06-19

Accepted: 2022-06-29

Published Online: 2022-07-22

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

https://doi.org/10.1515/eng-2022-0049

Keywords for this article

thermal load; residual force; thick cylinder; elastic–plastic analysis; plane strain

Creative Commons

BY 4.0