Using the penalty method

Consider an elastic-plastic analysis of the thick-wall tube shown. Compute stress distribution under plane strain conditions using
(a) axisymmetric 8-node elements,
(b) 3D elements with the symmetry condition enforced by the penalty method.
Compare the results.

$E=2\times10^5\text{ MPa},$ $\nu=0.3,$ $\sigma_Y=400\text{ MPa}.$ Prandtl–Reuss–von Mises elastic–perfectly plastic model.

Plane strain condition.

Internal pressure $p=287.47\text{ MPa}.$

• Axisymmetric solution, TUBE7P1
• Penalty solution, TUBE56P1

Denote by $x$ the radial coordinate, $r=1\text{ m}$ the inner radius, $r_Y$ the radius of plastic zone, and $R=2\text{ m}$ the outer radius of the tube, respectively. The elastic solution valid in the domain $x\ge r_Y$ satisfying $\sigma_r(R)=0$ at the boundary takes the form $$\begin{align} \sigma_r &= C\left[1-\left(\frac{R}{x}\right)^2\right] \\ \sigma_t &= C\left[1+\left(\frac{R}{x}\right)^2\right] \\ \sigma_o &= \nu(\sigma_r+\sigma_t) \end{align} \tag{1}\label{1}$$ where $C$ is a constant to be determined from the yield condition $\sigma_e(r_y)=\sigma_Y.$ Substituting \eqref{1} into von Mises' function $\sigma_e$ with $x=r_Y,$ we get $$C(r_y)=\sigma_Y\left[(1-2\nu)^2+3\left(\frac{R}{r_Y}\right)^4\right]^{-\frac{1}{2}}. \tag{2}\label{2}$$

A closed-form solution in the plastic zone $x<r_Y$ can only be obtained for Tresca's condition or if $\nu=1/2.$ With von Mises's criterion and $\nu=0.3,$ we must resort to some simplifying assumptions.

The state of plane strain enforces $$\varepsilon_o=\varepsilon_o^e+\varepsilon_o^p =\frac{1}{E}\left[\sigma_o-\nu(\sigma_r+\sigma_t)\right]+\varepsilon_o^p=0. \tag{3}\label{3}$$

At the elastic-plastic interface we have $\varepsilon_o^p=0,$ which gives $$\sigma_o=\nu(\sigma_r+\sigma_t)\quad\text{at}\quad x=r_Y. \tag{4}\label{4}$$ If the loading approaches its plastic limit, the stress components cease to change, therefore $\dot{\boldsymbol\varepsilon}^e=\mathbf{0}$ and \eqref{3} requires that $\dot\varepsilon_o=\dot\varepsilon_o^p=0.$ Using the Prandtl-Reuss equations $$\dot{\boldsymbol\varepsilon}^p=\frac{3}{2}\frac{\dot\varepsilon_p}{\sigma_Y}\mathbf{S} \tag{5}\label{5}$$ and calculating the axial component of deviatoric stress tensor $\mathbf{S}$ we derive $$S_o=\sigma_o-\textstyle\frac{1}{3}(\sigma_r+\sigma_t+\sigma_o)=0. \tag{6}$$

Hence $$\sigma_o=\textstyle\frac{1}{2}(\sigma_r+\sigma_t) \tag{7}\label{7}$$ for pronounced yielding. Equations \eqref{4} and \eqref{7} suggest that the axial stress is bounded in the plastic zone by $$\nu(\sigma_r+\sigma_t)<\sigma_o<\textstyle\frac{1}{2}(\sigma_r+\sigma_t) \tag{8}$$ when $\sigma_o=\nu(\sigma_r+\sigma_t)$ at $x=r_Y$ whereas the upper bound is approached in a limit as $\varepsilon_p\to\infty.$

Now, an approximation to the yielding function $\sigma_e$ is made by assuming $$\sigma_o=\textstyle\frac{1}{2}(\sigma_r+\sigma_t)\quad\text{for}\quad x<r_Y. \tag{9}\label{9}$$

Setting $\sigma_e=\sigma_Y$ we obtain the yield condition in the form $$\sigma_t-\sigma_r=\frac{2}{\sqrt{3}}\sigma_Y\quad\text{for}\quad x<r_Y. \tag{10}\label{10}$$

In view of \eqref{4} the maximum deviation should be expected in the vicinity of the elastic-plastic boundary. The error can easily be estimated by inserting elastic solution \eqref{1} and \eqref{2} into the expression $$\sigma_t-\sigma_r=2\sigma_Y\left[\left(1-2\nu\right)^2\left(\frac{R}{r_Y}\right)^4+3\right]^{-\frac{1}{2}}\quad\text{at}\quad x=r_Y. \tag{11}\label{11}$$

Clearly, if $\nu=0.5$ the right-hand sides of \eqref{10} and \eqref{11} are coincident. In this example $\nu=0.3,$ $R/r_Y=2/1.5,$ thus $$\sigma_t-\sigma_r=\frac{2}{1.87}\sigma_Y\quad\text{at}\quad x=r_Y. \tag{12}$$ Comparison with the factor $2/\sqrt{3}$ shows that \eqref{10} is a good approximation to the yield condition at any point in the plastic zone.

Solving the equilibrium equations $$x\frac{\mathrm{d}\sigma_r}{\mathrm{d}x}+\sigma_r-\sigma_t=0 \tag{13}$$ together with \eqref{10} and the boundary condition $\sigma_r(r)=-p$ yields $$\begin{align} \sigma_r &= -p+\displaystyle\frac{2}{\sqrt{3}}\sigma_Y\displaystyle\ln\frac{x}{r}, \\ \sigma_t &= -p+\displaystyle\frac{2}{\sqrt{3}}\sigma_Y\left(1+\displaystyle\ln\frac{x}{r}\right). \end{align} \tag{14}\label{14}$$ The relationship between $p$ and $r_Y$ can be established by equating \eqref{1}$_1$ and \eqref{14}$_1$ at $x=r_Y,$ from which $$\frac{p}{\sigma_Y}=\left[\left(\frac{R}{r_Y}\right)^2-1\right]\left[\left(1-2\nu\right)^2+3\left(\frac{R}{r_Y}\right)^4\right]^{-\frac{1}{2}}+\frac{2}{\sqrt{3}}\ln\frac{r_Y}{r}. \tag{15}$$ In this example we choose $r_Y=1.5\text{ m},$ thus $p=287.47\text{ MPa}.$ The estimate of $\sigma_o$ given by \eqref{9} is satisfactory for the approximation of the yielding function but it would have been very inaccurate for the calculation of $\sigma_o(x).$ Fortunately, a second-order correction can now be introduced by applying the Prandtl-Reuss equations. We use \eqref{9} and \eqref{14} as a first approximation to compute the deviatoric stress $$\begin{align} S_r &= -\sigma_Y/\sqrt{3}, \\ S_t &= \sigma_Y/\sqrt{3}, \\ S_o &= 0. \end{align} \tag{16}$$ Since the deviator is independent of plastic strain the Prandtl-Reuss equations \eqref{5} may easily be integrated $$\begin{align} \varepsilon_r^p &= -\varepsilon_p\sqrt{3}/2, \\ \varepsilon_t^p &= \varepsilon_p\sqrt{3}/2, \\ \varepsilon_o^p &= 0. \end{align} \tag{17}$$ By virtue of \eqref{3} the total strain components are computed as $$\begin{align} \varepsilon_r &= \displaystyle\frac{1-\nu^2}{E}\left[\sigma_r-\frac{\nu}{1-\nu}\sigma_t\right]+\displaystyle\frac{\sqrt{3}}{2}\varepsilon_p \\ \varepsilon_t &= \displaystyle\frac{1-\nu^2}{E}\left[\sigma_t-\frac{\nu}{1-\nu}\sigma_r\right]+\displaystyle\frac{\sqrt{3}}{2}\varepsilon_p \end{align} \tag{18}$$ and substituted into compatibility condition $$x\frac{\mathrm{d}\varepsilon_t}{\mathrm{d}x}+\varepsilon_t-\varepsilon_r=0. \tag{19}$$ With the aid of \eqref{14} and after some manipulations the differential equation for cumulated plastic strain $\varepsilon_p$ is obtained $$-x\frac{\mathrm{d}\varepsilon_p}{\mathrm{d}x}=2\varepsilon_p+\frac{8\left(1-\nu^2\right)}{3E}\sigma_Y. \tag{20}\label{20}$$ Solving \eqref{20} with the boundary condition $\varepsilon_p(r_Y)=0$ we arrive at $$\varepsilon_p=\frac{4(1-\nu^2)}{3E}\sigma_Y\left[\left(\frac{r_Y}{x}\right)^2-1\right]. \tag{21}$$ Once the cumulated plastic strain has been calculated we receive a better estimate of $\varepsilon_o^p.$ Integration of the Prandtl-Reuss equations \eqref{5} leads to $$\varepsilon_o^p = \frac{1}{\sigma_Y}\int_0^{\varepsilon_p}\left[\sigma_o-\textstyle\frac{1}{2}(\sigma_r+\sigma_t)\right]\mathrm{d}\varepsilon_p \simeq \frac{\varepsilon_p}{\sigma_Y}\left[\sigma_o-\textstyle\frac{1}{2}\left(\sigma_r+\sigma_t\right)\right]. \tag{22}$$ Finally, $\sigma_o$ is resolved from \eqref{3} as $$\sigma_o=\frac{1}{2}(\sigma_r+\sigma_t)\frac{2\nu\sigma_Y+\varepsilon_p}{\sigma_Y+E\varepsilon_p}. \tag{23}\label{23}$$ Note that $\sigma_o=\nu(\sigma_r+\sigma_t)$ if $\varepsilon_p=0$ and $\sigma_o=\frac{1}{2}(\sigma_r+\sigma_t)$ as $\varepsilon_p\to\infty.$ The results of numerical analysis and the analytical solution described by equations $\eqref{1},$ $\eqref{2},$ \eqref{14} and \eqref{23} are shown below.