# Introduction The main objective of this report is to analyze the convergence of the [ASME linearized stresses](https://www.ramsay-maunder.co.uk/knowledge-base/technical-notes/stress-linearisation-for-practising-engineers/) (membrane and membrane plus bending) in a pressurized infinitely-long circular pipe. This problem was chosen because 1. it is a typical problem tackled analytically in undergraduate courses in mechanical engineering (at least the displacements and stresses, not the ASME linearization), 2. it can be easily solved using standard finite-element tools ([@fig:quarter-caeplex]) so it serves as a benchmark problem, and 3. there is a practical interest to compare stress analysis results from general piping systems to the canonical case of the infinite pipe. ![A quarter of an infinite pipe solved with the cloud-based tool CAEplex. The project can be seen (and solved) online with just a regular web browser (even a mobile one) at ](quarter-caeplex.png){#fig:quarter-caeplex width=70%} In this work, the term convergence is taken as the study of how the linearized stresses computed with a finite-element method change with respect to the number of elements through the pipe thickness. Briefly, this technical report * analytically obtains the ASME-linearized (membrane and membrane plus bending) stresses for an infinite pipe subject to an uniform internal pressure, * serves as a benchmarking test, comparing the results obtained by the [free and open-source](https://en.wikipedia.org/wiki/Free_and_open-source_software) finite-element tool [Fino](](fino.svg){width=25%} # Problem {#sec:problem} Let us consider an infinite pipe (i.e. a cylinder) of internal radius $a$ and external radius $b$ with uniform mechanical properties---Young modulus $E$ and Poisson’s ratio $\nu$---subject to an internal uniform pressure $p$. We would like to compute the linearized membrane and membrane plus bending stresses according to ASME VIII Section 5 [@asme8div2, annex 5-A]. To obtain a numerical solution with the finite-element method, let us fix the numerical values of the parameters as follows \begin{align*} dnl b &= \SI{esyscmd([[echo "INCLUDE pipe-linearize/problem.was\nPRINT b" | wasora - | xargs]])}{\milli\meter} \\ dnl a &= \SI{esyscmd([[echo "INCLUDE pipe-linearize/problem.was\nPRINT a" | wasora - | xargs]])}{\milli\meter} \\ dnl E &= \SI{esyscmd([[echo "INCLUDE pipe-linearize/problem.was\nPRINT 1e-3*E" | wasora - | xargs]])}{\giga\pascal} \\ dnl \nu &= esyscmd([[echo "INCLUDE pipe-linearize/problem.was\nPRINT %.2f nu" | wasora - | xargs]]) \\ dnl p &= \SI{esyscmd([[echo "INCLUDE pipe-linearize/problem.was\nPRINT p" | wasora - | xargs]])}{\mega\pascal} \\ b &= esyscmd([[echo "INCLUDE pipe-linearize/problem.was\nPRINT b" | wasora - | xargs]])~\text{mm} \\ a &= esyscmd([[echo "INCLUDE pipe-linearize/problem.was\nPRINT a" | wasora - | xargs]])~\text{mm} \\ E &= esyscmd([[echo "INCLUDE pipe-linearize/problem.was\nPRINT 1e-3*E" | wasora - | xargs]])~\text{GPa} \\ \nu &= esyscmd([[echo "INCLUDE pipe-linearize/problem.was\nPRINT %.2f nu" | wasora - | xargs]]) \\ \end{align*} which correspond to a esyscmd([[grep inch pipe-linearize/problem.was | cut -c2-]]) carbon steel pipe according to ASME B36-10M [@asmeb36]. Again, for the sake of numerical results, we fix the internal pressure to $$ p = esyscmd([[echo "INCLUDE pipe-linearize/problem.was\nPRINT p" | wasora - | xargs]])~\text{MPa} $$ that is a typical pressure found in nuclear power plants (so ASME III applies). # Reference analytical solution {#sec:analytical} Given the cylindrical symmetry of the problem, there is only one independent variable, namely the radial coordinate $r$. Moreover, there are only two displacement fields that need to be considered: the axial $u_a(r)$ and the radial $u_r(r)$. The former is identically zero due to the fact that the cylinder is infinite ([plane strain condition](https://www.ramsay-maunder.co.uk/downloads/NBR03.pdf)). The equilibrium of momentum along the radial direction $r$ is $$ \frac{d\sigma_r}{dr} + \frac{\sigma_r(r) - \sigma_\theta(r)}{r} = 0 $${#eq:equilibrium} Defining the strains as \begin{align*} \epsilon_r &= \frac{du_r}{dr} \\ \epsilon_\theta &= \frac{u}{r} \end{align*} and the stresses as \begin{align*} \sigma_r &= \frac{E}{(1-\nu)(1-2\nu)} \cdot \Big[ (1-\nu) \cdot \epsilon_r + \nu \cdot \epsilon_\theta \Big] \\ \sigma_\theta &= \frac{E}{(1-\nu)(1-2\nu)} \cdot \Big[ \nu \cdot \epsilon_r + (1-\nu) \cdot \epsilon_r \Big] \\ \end{align*} then the differential equation [-@eq:equilibrium] can be casted in terms of the axial displacement $u_r(r)$ as $$ \frac{d^2 u}{dr^2} + \frac{1}{r} \cdot \frac{du}{dr} - \frac{u}{r^2} = 0 $$ that has the general solution $$ u(r) = c_1 \cdot r + \frac{c_2}{r} $$ For the boundary conditions of the particular problem that the radial stress should be equal to the negative of the internal pressure^[Exercise for the reader: why the negative of the pressure and not the pressure itself?] at $r=a$ and null at $r=b$, namely $$ \begin{cases} \sigma_r(a) &= -p\\ \sigma_r(b) &= 0 \end{cases} $$ the axial displacement has the particular solution $$ u_r(r) = p \cdot \frac{1+\nu}{E} \cdot \frac{a^2}{b^2-a^2} \cdot \left[ 1-2\nu + \frac{b^2}{r^2} \right]\cdot r $$ and the stresses the radial, tangential (hoop) and axial (longitudinal) stresses are $$ \sigma_r(r) = \frac{p \cdot a^2}{b^2-a^2} \cdot \left( 1 - \frac{b^2}{r^2}\right) $$ {#eq:sigmar} $$ \sigma_\theta(r) =\frac{p \cdot a^2}{b^2-a^2} \cdot \left( 1 + \frac{b^2}{r^2}\right) $$ {#eq:sigmatheta} $$ \sigma_l(r) = 2\nu \cdot \frac{p \cdot a^2}{b^2-a^2} $${#eq:sigmal} ## Stress linearization According to the ASME VIII code [@asme8div2], a _stress classification line_ should be chosen to compute the linearized stresses. Theses SCLs should be normal to iso-stress lines, which in this case may be any line radial from $r=a$ to $r=b$. The ASME code defines that each element of the _membrane_ stress tensor $\text{M}$ is \begin{equation*} \text{M}_{ij} = \frac{1}{b-a} \cdot \int_{a}^{b} \sigma_{ij}(r) \, dr \end{equation*} \noindent for $i=x,y,z$ and $j=x,y,z$. Having computed the nine values (of which only six are different) of such tensor, different kind of scalar membrane stresses can be obtained depending on how a single scalar stress “intensity” is computed out of the nine tensor elements. In particular we take into account the Tresca,^[Even though recent revisions of the ASME VIII code dropped the Tresca criterion in favor of Von Mises, the former is still the main way of computing the stress “intensity” in ASME III.] Von Mises and the three principal membrane stresses. As all the shear stresses are zero, $\sigma_l$, $\sigma_r$ and $\sigma_\theta$ happen to also be the principal stresses with \begin{align*} \sigma_1 &= \sigma_\theta \\ \sigma_2 &= \sigma_l \\ \sigma_3 &= \sigma_r \end{align*} So the membrane stresses to be considered are \begin{align*} \text{M}_1 &= \frac{1}{b-a} \cdot \int_{a}^{b} \sigma_\theta(r) \, dr \\ \text{M}_2 &= \frac{1}{b-a} \cdot \int_{a}^{b} \sigma_l(r) \, dr \\ \text{M}_3 &= \frac{1}{b-a} \cdot \int_{a}^{b} \sigma_r(r) \, dr \\ \text{M}_\text{tresca} &= \max\Big\{ \left| \text{M}_1 - \text{M}_2 \right|, \left| \text{M}_2 - \text{M}_3 \right|, \left| \text{M}_3 - \text{M}_1 \right| \Big\} \\ \text{M}_\text{vonmises} &= \sqrt{ \frac{(\text{M}_1-\text{M}_2)^2 + (\text{M}_2-\text{M}_3)^2 + (\text{M}_3-\text{M}_1)^2}{2}} \end{align*} In particular, replacing the stress distribution of the infinite pipe given by [@eq:sigmal; @eq:sigmar; @eq:sigmatheta] and carrying out the integrals, the principal membrane stresses are \begin{align*} \text{M}_1 &= {{a\,p}\over{b-a}} \\ \text{M}_2 &= {{2\,a^2\,\nu\,p}\over{b^2-a^2}} \\ \text{M}_3 &= {{a^2\,\left(2\,b-{{b^2+a^2}\over{a}}\right)\,p}\over{\left(b-a\right)\,\left(b^2-a^2\right)}} \end{align*} \medskip The other linearized stress, namely the _membrane plus bending_ stress tensor \text{MB}---again according to ASME VIII Annex 5-A---is $$ \text{MB}_{ij} = \text{M}_{ij} \pm \frac{6}{(b-a)^2} \cdot \int_{a}^{b} \sigma_{ij}(r) \cdot \left( \frac{a+b}{2}-r\right) \, dr $${#eq:mb} where the sign should be taken such that the resulting combined scalar stress (i.e. Tresca, Von Mises, etc.) represents the worst-case scenario. As before, the membrane plus bending scalar stresses are \begin{align*} \text{MB}_1 &= \text{M}_{1} \pm \frac{6}{(b-a)^2} \cdot \int_{a}^{b} \sigma_\theta(r) \cdot \left( \frac{a+b}{2}-r\right) \, dr \\ \text{MB}_2 &= \text{M}_{2} \pm \frac{6}{(b-a)^2} \cdot \int_{a}^{b} \sigma_l(r) \cdot \left( \frac{a+b}{2}-r\right) \, dr \\ \text{MB}_3 &= \text{M}_{3} \pm \frac{6}{(b-a)^2} \cdot \int_{a}^{b} \sigma_r(r) \cdot \left( \frac{a+b}{2}-r\right) \, dr \\ \text{MB}_\text{tresca} &= \max\Big\{ \left| \text{MB}_1 - \text{MB}_2 \right|, \left| \text{MB}_2 - \text{MB}_3 \right|, \left| \text{MB}_3 - \text{MB}_1 \right| \Big\} \\ \text{MB}_\text{vonmises} &= \sqrt{ \frac{(\text{MB}_1-\text{MB}_2)^2 + (\text{MB}_2-\text{MB}_3)^2 + (\text{MB}_3-\text{MB}_1)^2}{2}} \end{align*} For the infinite pipe case, $\text{MB}_1$ is maximum and $\text{MB}_2$ is minimum if we take the positive sign in [@eq:mb]. The integral in the expression for $\text{MB}_2$ above vanishes. The principal membrane plus bending stresses are then \begin{align*} \text{MB}_1 &= \displaystyle{{6\,a^2\,\left({{b^3+\left(2\,a\,\log a+a\right)\,b^2-a^2\,b}\over{2\,a}}-{{2\,b^2\,\log b+b^2}\over{2}}\right)\,p}\over{\left(b-a\right)^2\,\left(b^2-a^2\right)}}+{{a\,p}\over{b-a}} \\ \text{MB}_2 &= \displaystyle{{2\,a^2\,\nu\,p}\over{b^2-a^2}} \\ \text{MB}_3 &= \displaystyle{{6\,a^2\,\left({{2\,b^2\,\log b+b^2+2\,a\,b}\over{2}}-{{b^3+\left(2\,a\,\log a+a\right)\,b^2+a^2\,b}\over{2\,a}}\right)\,p}\over{\left(b-a\right)^2\,\left(b^2-a^2\right)}}+{{a^2\,\left(2\,b-{{b^2+a^2}\over{a}}\right)\,p}\over{\left(b-a\right)\,\left(b^2-a^2\right)}} \end{align*} For the particular parameters fixed in [@sec:problem], the numerical values for the theoretical linearized membrane and membrane plus bending scalar stresses are divert(-1) \begin{align*} \text{M}_\text{tresca} &= \SI{esyscmd(grep -A1 ev\(M_tresca pipe-linearize/analytical.txt | tail -n 1 | tr -d ' ' | awk -F\) '{printf("%.4f", $2) }')}{\mega\pascal} \\ \text{M}_\text{vonmises} &= \SI{esyscmd(grep -A1 ev\(M_vonmises pipe-linearize/analytical.txt | tail -n 1 | tr -d ' ' | awk -F\) '{printf("%.4f", $2) }')}{\mega\pascal} \\ \text{M}_1 &= \SI{esyscmd(grep -A1 ev\(M_1 pipe-linearize/analytical.txt | tail -n 1 | tr -d ' ' | awk -F\) '{printf("%.4f", $2) }')}{\mega\pascal} \\ \text{M}_2 &= \SI{esyscmd(grep -A1 ev\(M_2 pipe-linearize/analytical.txt | tail -n 1 | tr -d ' ' | awk -F\) '{printf("%.4f", $2) }')}{\mega\pascal} \\ \text{M}_3 &= \SI{esyscmd(grep -A1 ev\(M_3 pipe-linearize/analytical.txt | tail -n 1 | tr -d ' ' | awk -F\) '{printf("%.4f", $2) }')}{\mega\pascal} \\ \end{align*} and \begin{align*} \text{MB}_\text{tresca} &= \SI{esyscmd(grep -A1 ev\(MB_tresca pipe-linearize/analytical.txt | tail -n 1 | tr -d ' ' | awk -F\) '{printf("%.4f", $2) }')}{\mega\pascal} \\ \text{MB}_\text{vonmises} &= \SI{esyscmd(grep -A1 ev\(MB_vonmises pipe-linearize/analytical.txt | tail -n 1 | tr -d ' ' | awk -F\) '{printf("%.4f", $2) }')}{\mega\pascal} \\ \text{MB}_1 &= \SI{esyscmd(grep -A1 ev\(MB_1 pipe-linearize/analytical.txt | tail -n 1 | tr -d ' ' | awk -F\) '{printf("%.4f", $2) }')}{\mega\pascal} \\ \text{MB}_2 &= \SI{esyscmd(grep -A1 ev\(MB_2 pipe-linearize/analytical.txt | tail -n 1 | tr -d ' ' | awk -F\) '{printf("%.4f", $2) }')}{\mega\pascal} \\ \text{MB}_3 &= \SI{esyscmd(grep -A1 ev\(MB_3 pipe-linearize/analytical.txt | tail -n 1 | tr -d ' ' | awk -F\) '{printf("%.4f", $2) }')}{\mega\pascal} \\ \end{align*} respectively. divert(0) \begin{align*} \text{M}_\text{tresca} &= esyscmd(grep -A1 ev\(M_tresca pipe-linearize/analytical.txt | tail -n 1 | tr -d ' ' | awk -F\) '{printf("%.4f", $2) }')~\text{MPa} \\ \text{M}_\text{vonmises} &= esyscmd(grep -A1 ev\(M_vonmises pipe-linearize/analytical.txt | tail -n 1 | tr -d ' ' | awk -F\) '{printf("%.4f", $2) }')~\text{MPa} \\ \text{M}_1 &= esyscmd(grep -A1 ev\(M_1 pipe-linearize/analytical.txt | tail -n 1 | tr -d ' ' | awk -F\) '{printf("%.4f", $2) }')~\text{MPa} \\ \text{M}_2 &= esyscmd(grep -A1 ev\(M_2 pipe-linearize/analytical.txt | tail -n 1 | tr -d ' ' | awk -F\) '{printf("%.4f", $2) }')~\text{MPa} \\ \text{M}_3 &= esyscmd(grep -A1 ev\(M_3 pipe-linearize/analytical.txt | tail -n 1 | tr -d ' ' | awk -F\) '{printf("%.4f", $2) }')~\text{MPa} \\ \end{align*} and \begin{align*} \text{MB}_\text{tresca} &= esyscmd(grep -A1 ev\(MB_tresca pipe-linearize/analytical.txt | tail -n 1 | tr -d ' ' | awk -F\) '{printf("%.4f", $2) }')~\text{MPa} \\ \text{MB}_\text{vonmises} &= esyscmd(grep -A1 ev\(MB_vonmises pipe-linearize/analytical.txt | tail -n 1 | tr -d ' ' | awk -F\) '{printf("%.4f", $2) }')~\text{MPa} \\ \text{MB}_1 &= esyscmd(grep -A1 ev\(MB_1 pipe-linearize/analytical.txt | tail -n 1 | tr -d ' ' | awk -F\) '{printf("%.4f", $2) }')~\text{MPa} \\ \text{MB}_2 &= esyscmd(grep -A1 ev\(MB_2 pipe-linearize/analytical.txt | tail -n 1 | tr -d ' ' | awk -F\) '{printf("%.4f", $2) }')~\text{MPa} \\ \text{MB}_3 &= esyscmd(grep -A1 ev\(MB_3 pipe-linearize/analytical.txt | tail -n 1 | tr -d ' ' | awk -F\) '{printf("%.4f", $2) }')~\text{MPa} \\ \end{align*} respectively. See [@sec:maxima] and [@sec:wasora] for the input files used to compute both the integrals and the numerical values, either using [Maxima](http://maxima.sourceforge.net/) (symbolic integration) or [wasora](https://www.seamplex.com/wasora) [@wasoradesc] (numerical integration). # Parametric finite-element analysis solution In principle, the problem can be solved as a bi-dimensional plane-strain case. Even more, it is actually a one-dimensional problem with radial dependence only. Nevertheless, we solve a full three-dimensional pipe because we are interested in using the same set of methods and tools we would use in a real-case piping system. We do not even use any half or quarter symmetry but the full 360º annuli. In order to study the convergence of the linearized stresses computed using the stresses obtained in a finite-element solution, the [free and open source](https://en.wikipedia.org/wiki/Free_and_open-source_software) tool [Fino]( # Multi-freedom boundary condition for radial-only displacements {#sec:radial} To set the condition that in the external surfaces all displacement ought to be radially-symmetrical around the barycenter of the surface, a multi-freedom Dirichlet boundary condition needs to be set. That is to say, it is not the value of the actual displacements that is fixed but a certain relationship between the three displacements $u$, $v$ and $w$ and the spatial coordinates $x$, $y$ and $z$. In effect, let $\mathbf{x_0}=[\pm\ell/2,0,0]$ be the location of the barycenter of the face where the condition is to be applied. To force all displacements to be only radial, we require the vectors $\mathbf{u}=[u,v,w]$ and $(\mathbf{x}-\mathbf{x_0})$ to have the same direction. Using the scalar product notation: \begin{align*} \mathbf{u} \cdot \left(\mathbf{x}-\mathbf{x}_0\right) &= \| \mathbf{u} \| \cdot \| \left(\mathbf{x}-\mathbf{x}_0\right)\| \\ u \cdot (x-x_0) + v \cdot (y-y_0) + w \cdot (z-z_0) &= \sqrt{u^2 + v^2 + w^2} \cdot \sqrt{(x-x_0)^2 + (y-y_0)^2 + (z-z_0)^2} \end{align*} and taking into account that * $u=0$ because of the null axial displacement condition, * $\mathbf{x_0}=[\pm\ell/2,0,0]$, and * $x=\pm \ell/2$, \noindent then by squaring both members of the equation, we can simplify this expression to $$ (v y + w z)^2 = (v^2 + w^2) \cdot (y^2 + z^2) $$ Further algebraic manipulation leads to \begin{align*} (vy)^2 + 2 \cdot v y \cdot w z + (wz)^2 &= (vy)^2 + (vz)^2 + (wy)^2 + (wz)^2 \\ 2 \cdot v y \cdot w z &= (vz)^2 + (wy)^2 \\ 0 & = (vz - wy)^2 \end{align*} from which [@eq:radial] follows. \newpage
# Host, codes, scripts and input files {#sec:codes} The complete set of scripts and input files needed to run the cases discussed in this report can be found at :::: {text-center} :::: \medskip The results shown in this report correspond to revision ```{style=terminal} esyscmd(cd pipe-linearize; git log -1) ``` \noindent and were obtained using the free-as-in-free-speech^[And not just free as in “free beer”] tools [Gmsh](http://gmsh.info/) and [Fino](https://www.seamplex.com) using the following versions: ```{style=terminal} include(pipe-linearize/versions.txt) ``` These codes were executed on a host with the following features: ```{style=terminal} include(pipe-linearize/cpu.txt) ``` The complete case can be run by executing the following script: ```{.bash style=bash caption=run.sh} include(pipe-linearize/run.sh) ``` ## Analytical solutions with Maxima {#sec:maxima} The reference solutions presented in [@sec:analytical] where computed from the theoretical stresses by symbolic integration using [Maxima](http://maxima.sourceforge.net/). ```{.maxima style=maxima caption=problem.max} include(pipe-linearize/problem.max) ``` ```{.maxima style=maxima caption=analytical.max} include(pipe-linearize/analytical.max) ``` The results are illustrated in the following terminal mimic: ```{style=terminal} $ maxima -b analytical.max include(pipe-linearize/analytical.txt)dnl $ ``` ## Theoretical solutions with wasora {#sec:wasora} Instead of symbolic integration, adaptive numerical quadrature can be performed easily using [wasora](https://www.seamplex.com/wasora) [@wasoradesc]. First, we put all the problem parameters in a file `problem.was`: ```{.wasora style=wasora caption=problem.was} include(pipe-linearize/problem.was) ``` And then include this file from the main one `analytical.was`: ```{.wasora style=wasora caption=analytical.was} include(pipe-linearize/analytical.was) ``` ```{style=terminal} $ wasora analytical.was include(pipe-linearize/analytical.ppl)dnl $ ``` ## Geometry and mesh generation with Gmsh {#sec:gmsh} ```{.c style=gmsh caption=pipe.geo.m4} include(pipe-linearize/pipe.geo.m4) ``` ## Parametric finite-element analysis with Fino {#sec:fino} ```{.fino style=fino caption=pipe.fin} include(pipe-linearize/pipe.fin) ``` \newpage
# References