Temperature-Driven Instabilities in High-Pressure Vessel Flat Plates: A Thermal Buckling Study

1. Introduction

Designing pressure vessels for heating purposes presents significant challenges, especially in high-temperature environments and where corrosive substances are present [1]. The growing demand for advanced thermal processing systems underscores the need for robust heating surfaces that can withstand extreme conditions [2]. These systems typically involve heat transfer fluid on one side and processed materials on the other, necessitating materials with both corrosion resistance and mechanical strength at elevated temperatures (Figure 1). The presence of high-pressure fluids, combined with temperature gradients causing thermal expansion and resulting thermal stresses, further complicates the design. Additionally, ensuring corrosion resistance adds another layer of complexity to the configuration of heating surfaces [3].

In the conventional design of typical high-pressure vessels, which are commonly fabricated with curved bodies, the primary design consideration is often centered around the pressure exerted by the fluid contained within [4]. This assumption holds true for curved vessels where thermal expansion has minimal impact on the overall stress state, aligning well with ASTM recommendations that advocate for shell-type mesh in Finite Element Method (FEM) analysis of such vessels [5]. However, a critical question arises when high-pressure vessels incorporate flat components or are entirely composed of flat plates: can the same recommendation guarantee a reliable safety factor under these conditions?

This scenario presents a real-world industry problem where a high-pressure vessel, constructed from two parallel circular hollow plates interconnected by pins, contained a heating fluid under 50 psi pressure (see Figure 1). Initially, the system was assumed to maintain uniform temperature. Finite Element Method (FEM) modeling, employing a shell-type mesh as recommended by ASME, suggested a safety factor of 2. However, several weeks into operation, symmetrical patterns indicative of plastic deformations appeared on the heating plates (see Figure 2). Despite initial assessments using FEM modeling to ensure structural integrity, the design was compromised, resulting in permanent deformations of heating surfaces made from 304 stainless steel with a thickness of 0.5 inches. This begs the question: what was the root cause of this oversight? Which parameter was overlooked or underestimated? These were the inquiries driving engineers and researchers to seek answers.

Initial analyses initially suggested that friction between sliders and plates could be a potential cause of the observed deformations. However, this hypothesis was refuted by the absence of corresponding deformations in sliders with similar material properties. Instead, hypotheses surrounding localized temperature increases at contact points, possibly influenced by variations in material temperatures introduced by processed substances, have emerged as more plausible explanations. The irregular and symmetrical nature of the observed deformations, which do not align with the expected patterns from material conveyance mechanisms, further complicates our understanding of these failure modes.

Given that material design considerations must adhere to the linear elastic regime of stainless steel, we can address the problem by employing the superposition principle, dividing it into two distinct issues: one focusing on pressure and the other on temperature. Regarding pressure, Finite Element Method (FEM) modeling using a shell-type mesh has highlighted the critical role of pressure in design, providing assurance that the current design is safe for operation. However, a key question remains: does temperature contribute significantly to the design force or is its impact negligible?

Fortunately, the analysis of a hollow circular plate under thermal loading represents one of the few mechanical engineering problems with a viable analytical model [6]. This analytical model significantly enhances our understanding of the plate’s response to high temperatures, enabling us to explore the influence of temperature conditions on stress distribution within the plate. [7] have extensively studied the thermal deflections of axisymmetrically heated circular plates with fixed and simply supported edges. They approached the problem by considering the plate exposed to a temperature gradient as a 2D body. The resulting displacement and stress equations for a disk free of traction on the boundaries, subject to plane stress conditions due to thermal stresses, are as follows:

$${\sigma }_{rr}={\displaystyle \frac{E\alpha \left(T_a-T_b\right)}{2}}\left[\left(1-{\displaystyle \frac{b^2}{r^2}}\right){\displaystyle \frac{a^2}{a^2-b^2}}-{\displaystyle \frac{ln\left(\frac{r}{b}\right)}{ln\left(\frac{a}{b}\right)}}\right],$$

(1)

$${\sigma }_{\varphi \varphi }={\displaystyle \frac{E\alpha \left(T_a-T_b\right)}{2}}\left[\left(1+{\displaystyle \frac{b^2}{r^2}}\right){\displaystyle \frac{a^2}{a^2-b^2}}-{\displaystyle \frac{1+ln\left(\frac{r}{b}\right)}{ln\left(\frac{a}{b}\right)}}\right].$$

(2)

The parameters corresponding to Equations 1 and 2 are detailed in Table 1. It’s important to note that these equations apply specifically to the context of plane stress conditions within a thin plate, where the plate’s thickness is significantly smaller than its diameter—akin to the assumption underlying the shell-type mesh recommended by ASME for high-pressure vessels.

Expanding on Equations (3) and (4), it may initially seem that a uniform temperature ($T_a=T_b$) would not induce any stress within the plate. However, is this truly accurate? Can we continuously raise a plate’s temperature without inducing stress? The affirmative response is YES under ideal conditions of uniform temperature application and the absence of external perturbations, which would maintain the plate stress-free. Nevertheless, this ideal scenario is somewhat detached from reality, as the world rarely provides a disturbance-free environment. In essence, while the equilibrium equation may hold, the stability of this equilibrium under real-world conditions is a distinct matter—a concept known as the stability of equilibrium. In this context, it appears that temperature plays a critical role in inducing instability within high-pressure vessels, while pressure acts more as a disturbance force. This hypothesis forms the basis of the investigation that we aim to clarify in this article.

Another probable explanation for the observed deformations in the heating plates could be attributed to thermal buckling phenomena. Thermal buckling, characterized by symmetrical deformation patterns, is a credible scenario for several reasons. Studies have demonstrated that thermal buckling in circular plates can generate symmetric deformations under specific conditions [8]. Furthermore, thermal buckling can occur at stress levels lower than those typically captured by FEM simulations [9]. The disparity between FEM predictions and observed deformations underscores the credibility of this hypothesis. It is crucial to note that while FEM modeling predicts material behavior based on equilibrium equations, it may overlook the stability of these equations, particularly under thermal stress conditions. Therefore, an examination of the stability of equilibrium equations related to thermal stresses in circular plates could provide valuable insights into the root cause of these deformations.

This study seeks to deepen the understanding of mechanisms underlying thermal instabilities and deformations in high-temperature environments, providing insights into material behavior and critical design considerations essential for enhancing the reliability and performance of thermal processing systems.

This article presents the derivation of equilibrium equations governing thermal stress in circular annular plates (Section 2), followed by an evaluation of the stability of these equations (Section 3). Subsequently, it explores the factors that influence the thermal stability of these heating surfaces, along with an analysis of how design parameters impact thermal stability (Section 4). Finally, it provides a comprehensive discussion of critical design considerations (Section 5).

2. Mathematical Model of Circular Annular Plate

2.1. Kinetic and Kinematic Equations

Consider an annular plate of thickness h, inner radii b, and outer radii a, referred to as the polar coordinates (r,h,z), resting on an elastic foundation, as shown in Figure 3. In the present work, we assume that the modulus of elasticity E, thermal conductivity K, and the thermal expansion coefficient $\alpha$, Poisson’s ratio $\nu$ are considered to be constant across the thickness.

The displacement field for the thin plate ($h/a<0.05$) may be written as [10]

$$\begin{gathered} u(r,\theta ,z)=u_0(r,\theta )+z\varphi (r,\theta ), \\ v(r,\theta ,z)=v_0(r,\theta )+z\psi (r,\theta ), \\ w(r,\theta ,z)=w_0(r,\theta ). \end{gathered}$$

(3)

where $u_0(r,\theta )$, $v_0(r,\theta )$, and $w_0(r,\theta )$ represent the displacements on the middle sureface ($z=0$) along r, $\theta$, and z directions, respectively. In this work, the assumption of classical plate theory (CPT) is adapted and the effect of the transverse strains is ignored. Hence, based on the linear deformation assumptions in the polar coordinate, the out-of-plan strains (${\gamma }_{rz}$ and ${\gamma }_{\theta z}$) can be written as [10]:

$$\begin{gathered} {\gamma }_{rz}={\displaystyle \frac{\partial u}{\partial z}}+{\displaystyle \frac{\partial w_0}{\partial r}}=0, \\ {\gamma }_{\theta z}={\displaystyle \frac{\partial v}{\partial z}}+{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial w_0}{\partial \theta }}=0, \\ {\epsilon }_{zz}={\displaystyle \frac{\partial w_0}{\partial z}}=0. \end{gathered}$$

(4)

The last term of Equation (4) demonstrated that based on the assumption of CPT, the deflection does not vary through the thickness. In addition, the functions of $\varphi$ and $\psi$ can be described as follow [10]:

$$\begin{gathered} \varphi (r,\theta )=-{\displaystyle \frac{\partial w_0}{\partial r}}, \\ \psi (r,\theta )=-{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial w_0}{\partial \theta }}. \end{gathered}$$

(5)

Substituting Equation (5) into Equation (3) gives the displacement field equations for the thin circular plate as

$$\begin{gathered} u(r,\theta ,z)=u_0(r,\theta )-z{\displaystyle \frac{\partial w_0}{\partial r}}, \\ v(r,\theta ,z)=v_0(r,\theta )-{\displaystyle \frac{z}{r}}{\displaystyle \frac{\partial w_0}{\partial \theta }}, \\ w(r,\theta ,z)=w_0(r,\theta ). \end{gathered}$$

(6)

From Equation (6), the non-linear strain-displacement relations will be written as [11]:

$$\begin{gathered} {\epsilon }_{rr}={\displaystyle \frac{\partial u}{\partial r}}+{\displaystyle \frac{1}{2}}{\left({\displaystyle \frac{\partial w}{\partial r}}\right)}^2=\left({\displaystyle \frac{\partial u_0}{\partial r}}+{\displaystyle \frac{1}{2}}{\left({\displaystyle \frac{\partial w_0}{\partial r}}\right)}^2\right)-z\left({\displaystyle \frac{{\partial }^2{w}_0}{\partial r^2}}\right); \\ {\epsilon }_{\theta \theta }={\displaystyle \frac{u}{r}}+{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial v}{\partial \theta }}+{\displaystyle \frac{1}{2r^2}}{\left({\displaystyle \frac{\partial w}{\partial \theta }}\right)}^2=\left({\displaystyle \frac{u_0}{r}}+{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial v_0}{\partial \theta }}+{\displaystyle \frac{1}{2r^2}}{\left({\displaystyle \frac{\partial w_0}{\partial \theta }}\right)}^2\right)-z\left({\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial w_0}{\partial r}}+{\displaystyle \frac{1}{r^2}}{\displaystyle \frac{{\partial }^2{w}_0}{\partial {\theta }^2}}\right); \\ {\gamma }_{r\theta }=\left({\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial u}{\partial \theta }}+{\displaystyle \frac{\partial v}{\partial r}}-{\displaystyle \frac{v}{r}}+{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial w}{\partial \theta }}{\displaystyle \frac{\partial w}{\partial r}}\right), \\ =\left({\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial u_0}{\partial \theta }}-{\displaystyle \frac{z}{r}}{\displaystyle \frac{{\partial }^2{w}_0}{\partial \theta \partial r}}+{\displaystyle \frac{\partial v_0}{\partial r}}-{\displaystyle \frac{z}{r}}{\displaystyle \frac{{\partial }^2{w}_0}{\partial r\partial \theta }}+{\displaystyle \frac{z}{r^2}}{\displaystyle \frac{\partial w_0}{\partial \theta }}-{\displaystyle \frac{v_0}{r}}+{\displaystyle \frac{z}{r^2}}{\displaystyle \frac{\partial w_0}{\partial \theta }}+{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial w_0}{\partial \theta }}{\displaystyle \frac{\partial w_0}{\partial r}}\right), \\ =\left({\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial u_0}{\partial \theta }}+{\displaystyle \frac{\partial v_0}{\partial r}}-{\displaystyle \frac{v_0}{r}}+{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial w_0}{\partial \theta }}{\displaystyle \frac{\partial w_0}{\partial r}}\right)-2z\left({\displaystyle \frac{1}{r}}{\displaystyle \frac{{\partial }^2{w}_0}{\partial r\partial \theta }}-{\displaystyle \frac{1}{r^2}}{\displaystyle \frac{\partial w_0}{\partial \theta }}\right). \end{gathered}$$

(7)

Based on Equation (7), the strain components consist of two terms. Therefore, the stain components can be written as:

$$\begin{gathered} {\epsilon }_{rr}={\epsilon }_{rr}^0+z{\epsilon }_{rr}^1; \\ {\epsilon }_{\theta \theta }={\epsilon }_{\theta \theta }^0+z{\epsilon }_{\theta \theta }^1; \\ {\gamma }_{r\theta }={\gamma }_{r\theta }^0+z{\gamma }_{r\theta }^1, \end{gathered}$$

(8)

where ${\epsilon }_{rr}^0$, ${\epsilon }_{rr}^1$, ${\epsilon }_{\theta \theta }^0$, ${\epsilon }_{\theta \theta }^1$, ${\gamma }_{r\theta }^0$, and ${\gamma }_{r\theta }^1$ are defined by [12]:

$$\begin{gathered} {\epsilon }_{rr}^0={\displaystyle \frac{\partial u_0}{\partial r}}+{\displaystyle \frac{1}{2}}{\left({\displaystyle \frac{\partial w_0}{\partial r}}\right)}^2; \\ {\epsilon }_{rr}^1=-{\displaystyle \frac{{\partial }^2{w}_0}{\partial r^2}}; \\ {\epsilon }_{\theta \theta }^0={\displaystyle \frac{u_0}{r}}+{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial v_0}{\partial \theta }}+{\displaystyle \frac{1}{2r^2}}{\left({\displaystyle \frac{\partial w_0}{\partial \theta }}\right)}^2; \\ {\epsilon }_{\theta \theta }^1=-{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial w_0}{\partial r}}-{\displaystyle \frac{1}{r^2}}{\displaystyle \frac{{\partial }^2{w}_0}{\partial {\theta }^2}}; \\ {\gamma }_{r\theta }^0={\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial u_0}{\partial \theta }}+{\displaystyle \frac{\partial v_0}{\partial r}}-{\displaystyle \frac{v_0}{r}}+{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial w_0}{\partial \theta }}{\displaystyle \frac{\partial w_0}{\partial r}}; \\ {\gamma }_{r\theta }^1=-{\displaystyle \frac{2}{r}}{\displaystyle \frac{{\partial }^2{w}_0}{\partial r\partial \theta }}+{\displaystyle \frac{2}{r^2}}{\displaystyle \frac{\partial w_0}{\partial \theta }}. \end{gathered}$$

(9)

Within the framework of linear thermoelasticity, consistent with plane stress conditions, stress-strain relation are described as follow [13]:

$$\left\{\begin{array}{c}{\sigma }_{rr}\\ {\sigma }_{\theta \theta }\\ {\tau }_{r\theta }\end{array}\right\}={\displaystyle \frac{E}{1-{\nu }^2}}\left[\begin{array}{ccc}1& \nu & 0\\ \nu & 1& 0\\ 0& 0& {\displaystyle \frac{1-\nu }{2}}\end{array}\right]\left(\left\{\begin{array}{c}{\epsilon }_{rr}^0+z{\epsilon }_{rr}^1\\ {\epsilon }_{\theta \theta }^0+z{\epsilon }_{\theta \theta }^1\\ {\gamma }_{r\theta }^0+z{\gamma }_{r\theta }^1\end{array}\right\}-\Delta T\left\{\begin{array}{c}\alpha \\ \alpha \\ 0\end{array}\right\}\right),$$

(10)

where $\Delta T$ is the difference between the temperature distribution through the plate and the reference temperature.

Although the transverse strains (${\gamma }_{rz}$, ${\gamma }_{\theta z}$, and ${\epsilon }_{zz}$) are assumed to be identically zero, based on the kinematics assumption of the classical plate theory, their energy conjugates (${\tau }_{rz}$, ${\tau }_{\theta z}$, and ${\sigma }_{zz}$) in the reality are present to keep the plate in equilibrium. These stress components may be specified on the boundary. Thus, the transverse stresses must be accounted for in the boundary conditions and equilibrium of forces.

Stress resultants, based on the classical plate theory are related to the stress tensor components:

$$\begin{gathered} \left(N_{rr},N_{\theta \theta },N_{r\theta }\right)={\int }_{-{\displaystyle \frac{h}{2}}}^{{\displaystyle \frac{h}{2}}}\left({\sigma }_{rr},{\sigma }_{\theta \theta },{\tau }_{r\theta }\right)dz; \\ \left(M_{rr},M_{\theta \theta },M_{r\theta }\right)={\int }_{-{\displaystyle \frac{h}{2}}}^{{\displaystyle \frac{h}{2}}}z\left({\sigma }_{rr},{\sigma }_{\theta \theta },{\tau }_{r\theta }\right)dz; \\ \left(Q_r,Q_{\theta }\right)={\int }_{-{\displaystyle \frac{h}{2}}}^{{\displaystyle \frac{h}{2}}}\left({\tau }_{rz},{\tau }_{\theta z}\right)dz. \end{gathered}$$

(11)

Substituting Eqs. 7 and 8 into Equation (9) give the stress resultants in terms of the mid-plane displacement as

$$\left\{\begin{array}{c}{N}_{rr}\\ N_{\theta \theta }\\ N_{r\theta }\\ M_{rr}\\ M_{\theta \theta }\\ M_{r\theta }\end{array}\right\}={\displaystyle \frac{E}{1-{\nu }^2}}\left[\begin{array}{cccccc}h& \nu h& 0& 0& 0& 0\\ \nu h& h& 0& 0& 0& 0\\ 0& 0& \left({\displaystyle \frac{1-\nu }{2}}\right)h& 0& 0& 0\\ 0& 0& 0& {\displaystyle \frac{h^3}{12}}& {\displaystyle \frac{\nu h^3}{12}}& 0\\ 0& 0& 0& {\displaystyle \frac{\nu h^3}{12}}& {\displaystyle \frac{h^3}{12}}& 0\\ 0& 0& 0& 0& 0& {\displaystyle \frac{(1-\nu )h^3}{24}}\end{array}\right]\left\{\begin{array}{c}{\epsilon }_{rr}^0\\ {\epsilon }_{\theta \theta }^0\\ {\gamma }_{r\theta }^0\\ {\epsilon }_{rr}^1\\ {\epsilon }_{\theta \theta }^1\\ {\gamma }_{r\theta }^1\end{array}\right\}-\left\{\begin{array}{c}{N}^T\\ N^T\\ 0\\ M^T\\ M^T\\ 0\end{array}\right\},$$

(12)

where ${\epsilon }_{rr}^0$, ${\epsilon }_{rr}^1$, ${\epsilon }_{\theta \theta }^0$, ${\epsilon }_{\theta \theta }^1$, ${\gamma }_{r\theta }^0$, and ${\gamma }_{r\theta }^1$ are defined by Equation (9). $N^T$ and $M^T$ are the thermal force and moment resultants are calculated as

$$N^T={\int }_{{\displaystyle \frac{-h}{2}}}^{{\displaystyle \frac{h}{2}}}{\displaystyle \frac{\alpha E}{1-\nu }}\Delta Tdz,$$

(13)

$$M^T={\int }_{{\displaystyle \frac{-h}{2}}}^{{\displaystyle \frac{h}{2}}}{\displaystyle \frac{\alpha E}{1-\nu }}\Delta Tzdz.$$

(14)

If the temperature varies linearly within the thickness of the plate ($\Delta T=\left({\displaystyle \frac{T_b-T_t}{h}}\right)z+{\displaystyle \frac{T_b+T_t}{2}}$, in which $T_b$ and $T_t$ are the temperature of the bottom and top of the plate, respectively). Therefore, the thermal force and moment resultants can be calculated as:

$$N^T={\displaystyle \frac{\alpha Eh}{1-\nu }}\left({\displaystyle \frac{T_b+T_t}{2}}\right),$$

(15)

$$M^T={\displaystyle \frac{\alpha E{h}^2}{12(1-\nu )}}\left(T_b-T_t\right).$$

(16)

2.2. Equations of Motion

The equations of motion in classical plate theory are derived using the principle of virtual displacements. In these derivations, we incorporate thermal effects under the assumption that material properties remain independent of temperature, and the temperature T is known as a function of position ($\delta T=0$). Consequently, the temperature parameter only affects the formulation through constitutive equations (Equation (10)).

The dynamic version of the principle of virtual displacements or the minimum total potential energy, $\delta \Pi =\delta \left(U+V-K\right)$, is

$$0={\int }_0^t\left(\delta U+\delta V+\delta K\right)dt,$$

(17)

where U is the strain energy, V is the work done by external forces, and K is the kinetic energy. Suppose that $K_w$ and $K_g$ are the distributed linear and rotational springs of the foundation. The virtual strain energy is given by

$$\delta U={\int }_b^a{\int }_0^{2\pi }{\int }_{-{\displaystyle \frac{h}{2}}}^{{\displaystyle \frac{h}{2}}}\left({\sigma }_{rr}\delta {\epsilon }_{rr}+{\sigma }_{\theta \theta }\delta {\epsilon }_{\theta \theta }+{\tau }_{r\theta }\delta {\gamma }_{r\theta }\right)rdzd\theta dr$$

(18)

and the virtual work done by the external loads is given by

$$\delta V={\int }_b^a{\int }_0^{2\pi }\left(K_w{w}_0\delta w_0+K_g\left({\displaystyle \frac{\partial w_0}{\partial r}}\delta {\displaystyle \frac{\partial w_0}{\partial r}}+{\displaystyle \frac{1}{r^2}}{\displaystyle \frac{\partial w_0}{\partial \theta }}\delta {\displaystyle \frac{\partial w_0}{\partial \theta }}\right)\right)rd\theta dr.$$

(19)

To derive the equilibrium equations of the classical plate, we substitute Eqs. 15 and 16 into Equation (14), utilize Equation (12), integrate the displacement gradients by parts to accommodate virtual displacements, and then perform mathematical simplifications. This process yields the following equilibrium equations for the classical plate [14].

$$\begin{gathered} \delta u_0:{\displaystyle \frac{\partial N_{rr}}{\partial r}}+{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial N_{r\theta }}{\partial \theta }}+{\displaystyle \frac{1}{r}}\left(N_{rr}-N_{r\theta }\right)=0; \\ \delta v_0:{\displaystyle \frac{\partial N_{r\theta }}{\partial r}}+{\displaystyle \frac{2}{r}}{N}_{r\theta }+{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial N_{\theta \theta }}{\partial \theta }}=0; \\ \delta w_0:{\displaystyle \frac{{\partial }^2{M}_{rr}}{\partial r^2}}+{\displaystyle \frac{2}{r}}{\displaystyle \frac{\partial M_{rr}}{\partial r}}+{\displaystyle \frac{1}{r^2}}{\displaystyle \frac{{\partial }^2{M}_{\theta \theta }}{\partial {\theta }^2}}-{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial M_{\theta \theta }}{\partial r}}+{\displaystyle \frac{2}{r}}{\displaystyle \frac{{\partial }^2{M}_{r\theta }}{\partial r\partial \theta }}+{\displaystyle \frac{2}{r^2}}{\displaystyle \frac{\partial M_{r\theta }}{\partial \theta }} \\ +N_{rr}{\displaystyle \frac{{\partial }^2{w}_0}{\partial r^2}}+N_{\theta \theta }\left({\displaystyle \frac{1}{r^2}}{\displaystyle \frac{{\partial }^2{w}_0}{\partial r^2}}+{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial w_0}{\partial r}}\right) \\ +2N_{r\theta }\left({\displaystyle \frac{1}{r}}{\displaystyle \frac{{\partial }^2{w}_0}{\partial r\partial \theta }}-{\displaystyle \frac{1}{r^2}}{\displaystyle \frac{\partial w_0}{\partial \theta }}\right)-K_w{w}_0 \\ +K_g\left({\displaystyle \frac{{\partial }^2{w}_0}{\partial r^2}}+{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial w_0}{\partial r}}+{\displaystyle \frac{1}{r^2}}{\displaystyle \frac{{\partial }^2{w}_0}{\partial {\theta }^2}}\right)=0. \end{gathered}$$

(20)

3. Stability Analysis of Circular Annular Plate

For any type of boundary conditions, to produce the in-plane compressive thermal loads, immovability and periodicity conditions have to be satisfied. a superscript 0 indicates the pre-buckling conditions. Neglecting the lateral deflection of the plate in the pre-buckling state and solving the symmetrical type of the equilibrium equations yields

$$u_0^0(r,\theta )=0.$$

(21)

Therefore, by ignoring the non-linear terms of Equation (12), the pre-buckling forces are obtained

$$N_{rr}^0=N_{\theta \theta }^0=-{\displaystyle \frac{\alpha Eh}{1-\nu }}\left({\displaystyle \frac{T_b+T_t}{2}}\right),N_{r\theta }^0=0,$$

(22)

while the extra pre-buckling moments exist which are equal to

$$M_{rr}^0=M_{\theta \theta }^0=-{\displaystyle \frac{\alpha E{h}^2}{12(1-\nu )}}\left(T_b-T_t\right).M_{r\theta }^0=0.$$

(23)

Based on Equation (23), in the pre-buckling state, no deflection occurs in the plate due to thermal loading. In fact, the temperature gradient through the thickness of the plate is ignored. If an equilibrium state is assumed by displacement components $u_0^0$, $v_0^0$, and $w_0^0$, the displacement components of a neighbouring state of the stable equilibrium differ by $u_0^1$, $v_0^1$, and $w_0^1$ with respect to the equilibrium position. Thus, the total displacements of a neighbouring state are [15]

$$\begin{gathered} u_0=u_0^0+u_0^1, \\ {v}_0=v_0^0+v_0^1, \\ {w}_0=w_0^0+w_0^1. \end{gathered}$$

(24)

Similar to displacements, the stress resultants are divided into two terms representing the stable equilibrium and the neighbouring state. The stress resultants with superscript 1 are linear functions of displacement with superscript 1. Considering this, and using Eqs. 12 and 20, and performing proper simplifications, the stability equations become [16]

$$\begin{gathered} {\displaystyle \frac{\partial N_{rr}^1}{\partial r}}+{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial N_{r\theta }^1}{\partial \theta }}+{\displaystyle \frac{1}{r}}\left(N_{rr}^1-N_{r\theta }^1\right)=0; \\ {\displaystyle \frac{\partial N_{r\theta }^1}{\partial r}}+{\displaystyle \frac{2}{r}}{N}_{r\theta }^1+{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial N_{\theta \theta }^1}{\partial \theta }}=0; \\ {\displaystyle \frac{{\partial }^2{M}_{rr}^1}{\partial r^2}}+{\displaystyle \frac{2}{r}}{\displaystyle \frac{\partial M_{rr}^1}{\partial r}}+{\displaystyle \frac{1}{r^2}}{\displaystyle \frac{{\partial }^2{M}_{\theta \theta }^1}{\partial {\theta }^2}}-{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial M_{\theta \theta }^1}{\partial r}}+{\displaystyle \frac{2}{r}}{\displaystyle \frac{{\partial }^2{M}_{r\theta }^1}{\partial r\partial \theta }}+{\displaystyle \frac{2}{r^2}}{\displaystyle \frac{\partial M_{r\theta }^1}{\partial \theta }} \\ +N_{rr}^0{\displaystyle \frac{{\partial }^2{w}_0^1}{\partial r^2}}+N_{\theta \theta }^0\left({\displaystyle \frac{1}{r^2}}{\displaystyle \frac{{\partial }^2{w}_0^1}{\partial r^2}}+{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial w_0^1}{\partial r}}\right)+2N_{r\theta }^0\left({\displaystyle \frac{1}{r}}{\displaystyle \frac{{\partial }^2{w}_0^1}{\partial r\partial \theta }}-{\displaystyle \frac{1}{r^2}}{\displaystyle \frac{\partial w_0^1}{\partial \theta }}\right) \\ -K_w{w}_0^1+K_g\left({\displaystyle \frac{{\partial }^2{w}_0^1}{\partial r^2}}+{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial w_0^1}{\partial r}}+{\displaystyle \frac{1}{r^2}}{\displaystyle \frac{{\partial }^2{w}_0^1}{\partial {\theta }^2}}\right)=0. \end{gathered}$$

(25)

To obtain the stability equations based on the displacement field, Equation (12) has to be substituted into the above equations. Upon substitution, second and higher-order terms of incremental displacements may be omitted. The resulting equations are three stability equations based on the classical plate theory for a circular annular plate in contact with two parameters of elastic foundation

$$\begin{gathered} \left(u_{0,rr}^1+{\displaystyle \frac{1}{r}}{u}_{0,r}^1-{\displaystyle \frac{1}{r^2}}{u}_0^1-{\displaystyle \frac{1}{r^2}}{v}_{0,\theta }^1+{\displaystyle \frac{1}{r}}{v}_{0,r\theta }^1\right)+{\displaystyle \frac{1-\nu }{2}}\left({\displaystyle \frac{1}{r^2}}{u}_{0,\theta \theta }^1-{\displaystyle \frac{1}{r}}{v}_{0,r\theta }^1-{\displaystyle \frac{1}{r^2}}{v}_{0,\theta }^1\right)=0, \\ \left({\displaystyle \frac{1}{r^2}}{v}_{0,\theta \theta }^1+{\displaystyle \frac{1}{r}}{u}_{0,r\theta }^1+{\displaystyle \frac{1}{r^2}}{u}_{0,\theta }^1\right)+{\displaystyle \frac{1-\nu }{2}}\left(v_{0,rr}^1+{\displaystyle \frac{1}{r}}{v}_{0,r}^1-{\displaystyle \frac{1}{r^2}}{v}_0^1+{\displaystyle \frac{1}{r^2}}{u}_{0,\theta }^1-{\displaystyle \frac{1}{r}}{u}_{0,r\theta }^1\right)=0, \\ {\displaystyle \frac{E{h}^3}{12(1-{\nu }^2)}}\left(w_{0,rrrr}^1+{\displaystyle \frac{2}{r}}{w}_{0,rrr}^1-{\displaystyle \frac{1}{r^2}}{w}_{0,rr}^1+{\displaystyle \frac{1}{r^3}}{w}_{0,r}^1+{\displaystyle \frac{2}{r^2}}{w}_{0,rr\theta \theta }^1-{\displaystyle \frac{2}{r^3}}{w}_{0,r\theta \theta }^1+{\displaystyle \frac{4}{r^4}}{w}_{0,\theta \theta }^1+{\displaystyle \frac{1}{r^4}}{w}_{0,\theta \theta \theta \theta }^1\right) \\ +N^T\left(w_{0,rr}^1+{\displaystyle \frac{1}{r^2}}{w}_{0,\theta \theta }^1+{\displaystyle \frac{1}{r}}{w}_{0,r}^1\right)+K_w{w}_0^1-K_g\left(w_{0,rr}^1+{\displaystyle \frac{1}{r^2}}{w}_{0,\theta \theta }^1+{\displaystyle \frac{1}{r}}{w}_{0,r}^1\right)=0. \end{gathered}$$

(26)

Equation (26) shows that the stability equations are decoupled by themselves if the material properties of the plate are assumed to be constant through the thickness direction. The last relation of Equation (26) can be written as:

$$D_0\left({\nabla }^2\left({\nabla }^2{w}_0^1\right)\right)+\left(N^T-K_g\right)\left({\nabla }^2{w}_0^1\right)+K_w{w}_0^1=0$$

(27)

in which $D_0={\displaystyle \frac{E{h}^3}{12(1-{\nu }^2)}}$ is defined as a flexural rigidity of the plate. Also, the above equation can be written as:

$${\nabla }^4{w}_0^1+\left({\displaystyle \frac{N^T-K_g}{D_0}}\right){\nabla }^2{w}_0^1+\left({\displaystyle \frac{K_w}{D_0}}\right)w_0^1=0$$

(28)

For the sake of simplicity and generality, the following non-dimensional parameters are introduced

$$\begin{gathered} \overline{r}={\displaystyle \frac{r}{a}},\beta ={\displaystyle \frac{b}{a}},\delta ={\displaystyle \frac{h}{a}}, \\ {n}^T={\displaystyle \frac{N^T{a}^2}{D_0}},k_g={\displaystyle \frac{K_g{a}^2}{D_0}},k_w={\displaystyle \frac{K_w{a}^4}{D_0}}. \end{gathered}$$

(29)

By substituting Equation (29) into Eq 28, the stability equation of a circular annular plate in the form of non-dimensional can be written as:

$${\nabla }^4{\overline{w}}_0^1+\left(n^T-k_g\right){\nabla }^2{\overline{w}}_0^1+\left(k_w\right){\overline{w}}_0^1=0.$$

(30)

Hence, the thermal buckling condition can be driven by the equation below:

$$\left\{\left({\displaystyle \frac{{\partial }^2}{\partial {\overline{r}}^2}}+{\displaystyle \frac{1}{\overline{r}}}{\displaystyle \frac{\partial }{\partial \overline{r}}}+{\displaystyle \frac{1}{{\overline{r}}^2}}{\displaystyle \frac{{\partial }^2}{\partial {\theta }^2}}\right)\left({\displaystyle \frac{{\partial }^2{\overline{w}}_0^1}{\partial {\overline{r}}^2}}+{\displaystyle \frac{1}{\overline{r}}}{\displaystyle \frac{\partial {\overline{w}}_0^1}{\partial \overline{r}}}+{\displaystyle \frac{1}{{\overline{r}}^2}}{\displaystyle \frac{{\partial }^2{\overline{w}}_0^1}{\partial {\theta }^2}}+(n^T-k_g)\right)+k_w\right\}{\overline{w}}_0^1=0.$$

(31)

In the pressure vessels system under analysis, the pins between the plates act as an elastic foundation to stabilize the heating surface. Therefore, the values of $K_w$ and $K_g$ for the pin system can be simplified as:

$$K_w={\displaystyle \frac{A_p{E}_p}{L_p}}{N}_p,K_g={\displaystyle \frac{J_p{G}_p}{L_p}}{N}_p,$$

(32)

in which $A_p$ is the cross-section area of each pin, $J_p$ is the polar moment of inertia of each pin, $L_p$ is the length of pins, $N_p$ is the number of pins, and $E_p$ and $G_p$ are the elastic and shear modulus of the pins, respectively.

3.1. Solving the Stability Equation

For solving Equation (31) can be assumed that one suggested solution can be as [17]

$${\overline{w}}_0^1=W_n\left(\overline{r}\right)cos\left(n\theta \right),$$

(33)

where n is the number of nodal diameters. The value of $n=0$ indicates the symmetric buckled shape of the plate, while $n>0$ is associated with the asymmetric buckled shape. By substituting Equation (33) into Equation (31), it will be obtained:

$$\left({\displaystyle \frac{d^2}{d{\overline{r}}^2}}+{\displaystyle \frac{1}{\overline{r}}}{\displaystyle \frac{d}{d\overline{r}}}-{\displaystyle \frac{n^2}{{\overline{r}}^2}}\right)\left({\displaystyle \frac{d^2{W}_n\left(\overline{r}\right)}{d{\overline{r}}^2}}+{\displaystyle \frac{1}{r}}{\displaystyle \frac{d{W}_n\left(\overline{r}\right)}{d\overline{r}}}+\left(n^T-{\displaystyle \frac{n^2}{{\overline{r}}^2}}\right)\right)W_n\left(\overline{r}\right)cos\left(n\theta \right)=0.$$

(34)

A non-trivial solution for Equation (34) is [17]

$$W_n\left(\overline{r}\right)=C_1{J}_n\left(\sqrt{n^T}\overline{r}\right)+C_2{Y}_n\left(\sqrt{n^T}\overline{r}\right)+C_3\overline{r}+C_4\left\{\begin{array}{c}ln\left(\overline{r}\right)n=0\\ \overline{r^{-n}}n>0\end{array}\right\},$$

(35)

in which $J_n\left(\sqrt{n^T}\overline{r}\right)$ and $Y_n\left(\sqrt{n^T}\overline{r}\right)$ are the Bessel functions of the first and second kind, respectively, and $C_i$, $i=1,\cdots ,4$, are constants that are evaluated by applying the associated boundary conditions. The top form of Equation (35) ($n=0$) is associated with symmetric buckling, while the bottom one ($n>0$) is related to asymmetric buckling.

As the heating surfaces are clamped for both inner and outer edges, they exhibit the bifurcation type buckling for transverse thermal loading. Therefore, the boundary conditions are:

$$W_n\left(1\right)={\displaystyle \frac{d{W}_n}{d\overline{r}}}\left(1\right)=W_n\left(\beta \right)={\displaystyle \frac{d{W}_n}{d\overline{r}}}\left(\beta \right)=0.$$

(36)

By applying boundary conditions (Eq.36) to Equation (35), a system of homogeneous equations is formed. For the symmetric buckling $(n=0)$:

$$\left[\begin{array}{cccc}{J}_0\left(\sqrt{n^T}\right)& Y_0\left(\sqrt{n^T}\right)& 1& 0\\ -\sqrt{n^T}{J}_1\left(\sqrt{n^T}\right)& -\sqrt{n^T}{Y}_1\left(\sqrt{n^T}\right)& 1& 1\\ J_0\left(\beta \sqrt{n^T}\right)& Y_0\left(\beta \sqrt{n^T}\right)& 1& ln\left(\beta \right)\\ -\sqrt{n^T}{J}_1\left(\beta \sqrt{n^T}\right)& -\sqrt{n^T}{Y}_1\left(\beta \sqrt{n^T}\right)& 1& {\displaystyle \frac{1}{\beta }}\end{array}\right]\left[\begin{array}{c}{C}_1\\ C_2\\ C_3\\ C_4\end{array}\right]=0$$

(37)

and for asymmetric buckling ($n>0$):

$$\left[\begin{array}{cccc}{J}_n\left(\sqrt{n^T}\right)& Y_n\left(\sqrt{n^T}\right)& 1& 1\\ \sqrt{n^T}\left({\displaystyle \frac{J_{n-1}-J_{n+1}}{2}}\right)& \sqrt{n^T}\left({\displaystyle \frac{Y_{n-1}-Y_{n+1}}{2}}\right)& 1& -n\\ J_n\left(\beta \sqrt{n^T}\right)& Y_n\left(\beta \sqrt{n^T}\right)& 1& {\beta }^{-n}\\ \sqrt{n^T}\left({\displaystyle \frac{J_{n-1}-J_{n+1}}{2}}\right)& \sqrt{n^T}\left({\displaystyle \frac{Y_{n-1}-Y_{n+1}}{2}}\right)& 1& {\displaystyle \frac{-n}{{\beta }^{n+1}}}\end{array}\right]\left[\begin{array}{c}{C}_1\\ C_2\\ C_3\\ C_4\end{array}\right]=0.$$

(38)

To gain a non-trivial solution of Eqs. 37 and 38, the determinant of the coefficients matrices have to be set equal to zero, which yields a non-linear equation containing n and $n^T$. Now to find the non-dimensional critical buckling loads of the plate for every positive integer number n, the associated determinant equations have to be solved. This equation determines the stability criterion for the thermal buckling of circular annular plates. Finding the smallest positive root of the associated equations for each n and choosing the smallest between them, yields the associated critical value of $n^T$, which is called $n_{cr}^T$.

The amount of Eqs. 37 and 38 for $n=0,\cdots ,3$ are depicted in Figure 4, Figure 5, Figure 6 and Figure 7. The borders between the yellow and the blue in the figures indicate the zero amount of the determinant of Eqs. 37 and 38. From Figure 4, Figure 5, Figure 6 and Figure 7 show that the minimum amount of $n^T$, which is called $n_{cr}^T$, happened when $n=0$.

3.2. Case Study: A Real Heating Surface with Molten Salt

To explore the initial hypothesis regarding plastic deformation of a heating plate, a practical case study was conducted under controlled conditions using molten salt as a heat transfer fluid to ensure uniform radial temperature distribution [18]. The selection of 304 stainless steel (SS304) for the heating plate reflects a careful balance between corrosion resistance and cost-effectiveness in industrial applications [19,20]. This study exemplifies thermal stability analysis in a real-world heat transfer system, focusing on the behavior of a specific heating plate under operational conditions. By investigating thermal buckling phenomena, the analysis demonstrates that the heating plate operates within an unstable zone, where any perturbations could disrupt its equilibrium and potentially lead to plastic deformation or failure.

In this study, the heating plate is assumed to be clamped at both its inner and outer edges due to welding at these boundaries. The inner diameter ($b=34$ inches) and outer diameter ($a=96$ inches) of the plate configuration are chosen to represent typical dimensions encountered in industrial heating surface designs. The thermal properties of SS304, including its elastic modulus ($E=200$ GPa), Poisson’s ratio ($\nu =0.265$), and thermal expansion coefficient ($\alpha =17.3\times {10}^{-6}$${\mathrm{K}}^{-1}$), provide critical insights into material behavior under high-temperature conditions [21]. This case study serves as a practical demonstration of thermal instability in heating surface configurations, underscoring the importance of robust design considerations for ensuring structural integrity and performance in heat transfer systems.

To perform thermal stability analysis for the , we substitute $\beta =0.354167$ into Equations (37) and (38) and utilize MAPLE software to solve these equations, yielding the non-dimensional critical buckling parameter $\sqrt{n^T}$. This parameter is then used to calculate the critical temperature of buckling for the circular annular plate using the formula:

$$\Delta T_{cr}={\displaystyle \frac{h^2}{6\alpha \left(1+\nu \right)a^2}}{n}^T$$

(39)

The results of these calculations are summarized in Table 2.

4. Design Parameters Analysis

Designing plates for high-temperature applications presents significant challenges in selecting appropriate dimensions and materials. In this section, key design parameters essential for the design of pressure vessels operating at elevated temperatures are explored. Specifically, the scalability of the vessel, the utilization of alternative stainless steel compositions, and the impact of increased plate thickness are investigated. Each of these design parameters influences the critical thermal buckling temperature, providing valuable insights into optimizing plate design for high-temperature environments.

4.1. Scalability Assessment

The scalability of a heating system in certain applications could offer economic benefits. For instance, employing larger heating plates could enhance the overall heating process efficiency. Therefore, assessing the performance of larger heating plates under similar clamping conditions at both inner and outer edges, constrained by welding considerations, is advantageous. The material of choice is assumed to be SS304, with specified inner and outer diameters of $b=42$ inches and $a=144$ inches, respectively.

To assess the critical temperature for thermal buckling of these larger plates, the non-dimensional critical buckling parameter $\sqrt{n^T}$ is calculated using the value $\beta =0.29167$, substituted into Eqs. 37 and 38. The resulting values are computed using MAPLE software to determine the critical buckling temperature difference $\Delta T_{cr}$ based on Equation (39). The summarized outcomes of these calculations are presented in Table 3 and Table 4 for plate thicknesses of 0.375 inches and 0.5 inches, providing insights into the critical buckling behavior under varying thicknesses for the planned larger heating plates. These findings contribute to the ongoing design considerations aimed at optimizing the thermal stability and performance of the heating surfaces.

By comparing the critical thermal buckling temperatures of the initial case study problem, where the inner and outer diameters were 34 and 96 inches (as summarized in Table 3), with the current scaled-up heating plate featuring inner and outer diameters of 42 and 144 inches (as summarized in Table 3), it can be concluded that increasing the outer diameter may decrease the critical temperature for thermal buckling of the plate. This finding is valuable for engineers seeking to balance economic efficiency with ensuring the stability of heating surfaces.

4.2. Effect of Thermal Expansion Coefficient ($\alpha$)

Given the presence of molten salts and the associated high risk of corrosion, the selection of materials for heating surfaces is considerably restricted. Consequently, our options for heating plates are primarily limited to austenitic stainless steel due to its robust performance under such demanding conditions. Austenitic stainless steels are renowned for their ability to withstand high temperatures while retaining excellent corrosion resistance.

To facilitate informed material selection for future heating plate designs, Table 5 summarizes some of the most commonly used austenitic stainless steels along with their respective critical buckling temperatures. This information serves as a valuable reference for identifying suitable materials that can effectively withstand the thermal and corrosive challenges posed by molten salt environments, supporting the development of durable and reliable heating surfaces for industrial applications.

As is shown in Table 5, using the SS310 instead of SS304 will increase the critical temperature of buckling by 20%. Therefore, it can be suggested as a solution for enhancing the critical temperature difference for the buckling of heating plates.

4.3. Effect of Thickness of Plate (h)

The thickness of the plate not only impacts the cost and weight of the system but also plays a critical role in influencing the thermal buckling behavior of heating surfaces. In this section, we explore the effect of thickness on thermal buckling. As shown in Figure 8, increasing the plate thickness enhances its resistance to buckling, as evidenced by the relationship between thickness and critical buckling temperature difference illustrated in the figure.

Based on Figure 8, it is evident that thicker plates exhibit higher critical buckling temperatures, indicating improved buckling resistance with increased thickness. However, the relationship between thickness and critical thermal buckling temperature is nonlinear. For instance, increasing the thickness from 1 inch to 2 inches raises the critical thermal buckling temperature by 78°C, whereas increasing the thickness from 4 inches to 5 inches results in a 256 °C increase in critical temperature. Therefore, one effective approach to mitigate thermal buckling in heating surfaces could be to increase plate thickness. This strategy enhances the structural integrity of heating surfaces, reducing the risk of buckling under high-temperature conditions. Engineers can leverage this insight when designing high-temperature flanges to prevent mechanical instability.

4.4. Effect of the inner and outer diameter ratio ($\beta$)

Another variable that significantly influences the critical thermal buckling temperature is the ratio of the inner diameter to the outer diameter of the plate. This ratio, depicted in Figure 9, demonstrates the relationship between the critical temperature and the diameter proportions for plates with thicknesses of 0.375 inches and 0.5 inches.

In Figure 9, we observe that the effect of the diameter ratio on the critical temperature is most pronounced when the ratio exceeds 0.8. Below this threshold, the influence of the diameter ratio on the critical temperature appears to be minimal or negligible. This observation underscores the importance of considering the diameter ratio when designing heating surfaces to mitigate thermal buckling, particularly for ratios exceeding 0.8, where changes in the diameter proportions can notably impact the critical buckling temperature. Understanding this relationship is crucial for optimizing the design parameters of heating surfaces to ensure stability and performance under high-temperature conditions. This insight can also be valuable for engineering applications, particularly in designing flanges for high temperatures to mitigate the risk of permanent deformation and prevent seal failure under elevated temperature conditions.

5. Conclusions

The aim of this article is to demonstrate how elevated temperatures in pressure vessels composed of flat plates can induce an unstable condition, where any perturbation could lead to structural failure. This study highlights that, even when adhering to ASME guidelines for employing shell-type mesh in FEM simulations of pressure vessels, the critical influence of temperature on destabilizing equilibrium is not adequately captured by the ASME-recommended FEM modeling.

In conclusion, this work has investigated the potential for thermal buckling in pressure vessels at elevated temperatures, yielding several key findings:

• The critical buckling temperature of the heating surface increases with the clamped edge of the plate. This suggests that fixing the boundaries of the plate raises the threshold temperature at which buckling occurs. Analogously, a column with clamped ends exhibits a higher critical force before buckling.
• Increasing the thickness of the plate enhances its resistance to thermal buckling, reducing the likelihood of buckling as thickness increases.
• Holding the outer diameter of the plate constant while increasing the inner diameter results in a higher critical buckling temperature.
• Choosing materials with lower thermal expansion coefficients can improve the thermal buckling resistance of heating surfaces.
• Clamping the plates together, effectively creating an elastic foundation through additional pinning and fixing between the plates, increases the critical buckling temperature.

These findings offer crucial insights for optimizing the design and material selection of heating surfaces, thereby reducing the risk of thermal buckling and enhancing the reliability and performance of high-temperature processing systems. Further research and experimentation in these areas can lead to more effective strategies for designing robust heating surfaces capable of withstanding extreme thermal conditions.

Engineers and designers should prioritize understanding thermal buckling phenomena when designing high-temperature pressure vessels. Based on the findings of this study, it is recommended to explore alternative mesh types, such as 3D solid mesh, for FEM modeling of pressure vessels composed of flat plates operating at high temperatures. This approach will help ensure accurate safety predictions and enhance the reliability of pressure vessel designs under extreme thermal conditions.