Mar APR May 5 2013 2014 2015
1 capture 5 Apr 14 - 5 Apr 14

Close Help

  
Using the penalty method

Problem description: Consider an elasto-plastic analysis of the thick-wall tube shown. Compute stress distribution under plane strain conditions using i) the axisymmetric 8-node elements and ii) the 3D elements with the symmetry condition enforced by the penalty method. Compare the results.

Mesh: Use the following element types:
ITE=7--see appendix A.3, ten quadrilaterals
ITE=56--see appendix B.5, ten hexahedra.

Material properties: $E=2\times 10^5$ MPa, $\nu=0.3$, $\sigma_Y=400$ MPa,
Prandtl-Reuss-von Mises elastic-perfectly plastic model.

Support: Plane strain condition.

Loading: The internal pressure p=287.47 MPa.

Solution: Denote by x the radial coordinate, r=1m the inner radius, r_Y the radius of plastic zone, and R=2m 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{array}{rcl} \sigma_r&=&C[1-\left(\frac{R}{x}\right)^2... ...}{x}\right)^2]\\ \sigma_o&=&\nu(\sigma_r+\sigma_t) \end{array}$$ (V.1)

C is a constant to be determined from the yield condition $\sigma_e(r_y)=\sigma_Y$. Substituting (V.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}}$$ (V.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

$$\epsilon_o=\epsilon_o^e+\epsilon_o^p= \frac{1}{E}[\sigma_o-\nu(\sigma_r+\sigma_t)]+\epsilon_o^p=0$$ (V.3)

At the elastic-plastic interface we have $\epsilon_o^p=0$, which gives

$$\sigma_o=\nu(\sigma_r+\sigma_t)~~\mbox{at}~~x=r_Y$$ (V.4)

If the loading approaches its plastic limit, the stress components cease to change, therefore $\mbox{\boldmath$\dot\epsilon$ }^e={\bf0}$ and (V.3) requires that $\dot\epsilon_o=\dot\epsilon_o^p=0$. Using the Prandtl-Reuss equations

$$\mbox{\boldmath$\dot\epsilon$ }^p=\frac{3}{2}\frac{\dot\epsilon_p}{\sigma_Y}{\bf S}$$ (V.5)

and calculating the axial component of deviatoric stress tensor ${\bf S}$ we derive

$$S_o=\sigma_o-\textstyle\frac{1}{3}(\sigma_r+\sigma_t+\sigma_o)=0$$ (V.6)

Hence

$$\sigma_o=\textstyle\frac{1}{2}(\sigma_r+\sigma_t)~~ \mbox{for pronounced yielding}$$ (V.7)

Equations (V.4) and (V.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)$$ (V.8)

when $\sigma_o=\nu(\sigma_r+\sigma_t)$ at x=r_Y whereas the upper bound is approached in a limit as $\epsilon_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)~~\mbox{for}~~x<r_Y$$ (V.9)

Setting $\sigma_e=\sigma_Y$ we obtain the yield condition in the form

$$\sigma_t-\sigma_r=\frac{2}{\sqrt{3}}\sigma_Y~~\mbox{for}~~x<r_Y$$ (V.10)

In view of (V.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 (V.1), (V.2) into the expression

$$\sigma_t-\sigma_r=2\sigma_Y \left[(1-2\nu)^2\left(\frac{R}{r_Y}\right)^4+3\right]^{-\frac{1}{2}} ~~\mbox{at}~~x=r_Y$$ (V.11)

Clearly, if $\nu=0.5$ the right-hand sides of (V.10) and (V.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~~\mbox{at}~~x=r_Y$$ (V.12)

Comparison with the factor $2/\sqrt{3}$ shows that (V.10) is a good approximation to the yield condition at any point in the plastic zone.

Solving the equilibrium equations

$$x\frac{\,\mbox{\rm d}\sigma_r}{\,\mbox{\rm d}x} +\sigma_r-\sigma_t=0$$ (V.13)

together with (V.10) and the boundary condition $\sigma_r(r)=-p$ yields

$$\begin{array}{rcl} \sigma_r &=& -p+\displaystyle\frac{2}{\sq... ...2}{\sqrt{3}}\sigma_Y(1+\displaystyle\ln\frac{x}{r}) \end{array}$$ (V.14)

The relationship between p and r_Y can be established by equating (V.1)₁ and (V.14)₁ at x=r_Y, from which

$$\frac{p}{\sigma_Y}= \left[\left(\frac{R}{r_Y}\right)^2-1\rig... ...t)^4\right]^{-\frac{1}{2}} +\frac{2}{\sqrt{3}}\ln\frac{r_Y}{r}$$ (V.15)

In this example we choose r_Y=1.5m, thus $p=287.47\,$MPa.

The estimate of $\sigma_o$ given by (V.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 (V.9) and (V.14) as a first approximation to compute the deviatoric stress

$$\begin{array}{rcr} S_r &=& -\sigma_Y/\sqrt{3}\\ S_t &=& \sigma_Y/\sqrt{3}\\ S_o &=&0 \end{array}$$ (V.16)

Since the deviator is independent of plastic strain the Prandtl-Reuss equations (V.5) may easily be integrated

$$\begin{array}{rcr} \epsilon_r^p &=& -\epsilon_p\sqrt{3}/2\\ ... ...n_t^p &=& \epsilon_p\sqrt{3}/2\\ \epsilon_o^p &=&0 \end{array}$$ (V.17)

By virtue of (V.3) the total strain components are computed as

$$\begin{array}{rcl} \epsilon_r &=&\displaystyle\frac{1-\nu^2}... ...\sigma_r]+\displaystyle\frac{\sqrt{3}}{2}\epsilon_p \end{array}$$ (V.18)

and substituted into compatibility condition

$$x\frac{\,\mbox{\rm d}\epsilon_t}{\,\mbox{\rm d}x}+\epsilon_t-\epsilon_r=0$$ (V.19)

With the aid of (V.14) and after some manipulations the differential equation for cumulated plastic strain $\epsilon_p$ is obtained

$$-x\frac{\,\mbox{\rm d}\epsilon_p}{\,\mbox{\rm d}x}=2\epsilon_p+\frac{8(1-\nu^2)}{3E}\sigma_Y$$ (V.20)

Solving (V.20) with the boundary condition $\epsilon_p(r_Y)=0$ we arrive at

$$\epsilon_p=\frac{4(1-\nu^2)}{3E}\sigma_Y \left[\left(\frac{r_Y}{x}\right)^2-1\right]$$ (V.21)

Once the cumulated plastic strain has been calculated we receive a better estimate of $\epsilon_o^p$. Integration of the Prandtl-Reuss equations (V.5) leads to

$$\epsilon_o^p=\frac{1}{\sigma_Y}\int_0^{\epsilon_p} [\sigma_o... ...}{\sigma_Y}[\sigma_o-\textstyle\frac{1}{2}(\sigma_r+\sigma_t)]$$ (V.22)

Finally, $\sigma_o$ is resolved from (V.3) as

$$\sigma_o=\frac{1}{2}(\sigma_r+\sigma_t) \frac{2\nu\sigma_Y+\epsilon_p}{\sigma_Y + E \epsilon_p}$$ (V.23)

Note that $\sigma_o=\nu(\sigma_r+\sigma_t)$ if $\epsilon_p=0$ and $\sigma_o=\frac{1}{2}(\sigma_r+\sigma_t)$ as $\epsilon_p\to\infty$.

The results of numerical analysis and the analytical solution described by eqns. (V.1),(V.2), (V.14) and (V.23) are shown in figure.