more changes from angus suggestions

nafems4.md

@@ -24,7 +24,7 @@ An infinitely-long pressurised thick pipe as taught in college and an isometric
 
 Like the pendulums above, we will be swinging back and forth between a case study about fatigue analysis in piping systems of a nuclear power plant and more generic and even romantic topics related to finite elements and computational mechanics. These latter regressions will not remain just as abstract theoretical ideas. Not only will they be directly applicable to the development of the main case, but they will also apply to a great deal of other engineering problems tackled with the finite element method (and its cousins, [@sec:formulations]).
 
-Finite elements are like magic to me. I mean, I can follow the whole derivation of the equations, from the strong, weak and variational formulations of the equilibrium equations for the mechanical problem (or the energy conservation for heat transfer) down to the [algebraic multigrid](https://en.wikipedia.org/wiki/Multigrid_method) [preconditioner](https://en.wikipedia.org/wiki/Preconditioner) for the inversion of the [stiffness matrix](https://en.wikipedia.org/wiki/Stiffness_matrix) passing through [Sobolev spaces](https://en.wikipedia.org/wiki/Sobolev_space) and the [grid generation](https://en.wikipedia.org/wiki/Mesh_generation). Then I can sit down and program all these steps into a computer, including the shape functions and its derivatives, the assembly of the discretised stiffness matrix ([@sec:building]), the numerical solution of the system of equations ([@sec:solving]) and the computation of the gradient of the solution ([@sec:stress-computation). Yet, the fact that all these a-priori unconnected steps give rise to pretty pictures that resemble reality is still astonishing to me.
+Finite elements are like magic to me. I mean, I can follow the whole derivation of the equations, from the strong, weak and variational formulations of the equilibrium equations for the mechanical problem (or the energy conservation for heat transfer) down to the [algebraic multigrid](https://en.wikipedia.org/wiki/Multigrid_method) [preconditioner](https://en.wikipedia.org/wiki/Preconditioner) for the inversion of the [stiffness matrix](https://en.wikipedia.org/wiki/Stiffness_matrix) passing through [Sobolev spaces](https://en.wikipedia.org/wiki/Sobolev_space) and the [grid generation](https://en.wikipedia.org/wiki/Mesh_generation). Then I can sit down and program all these steps into a computer, including the [shape functions](https://en.wikipedia.org/wiki/Finite_element_method_in_structural_mechanics#Interpolation_or_shape_functions) and its derivatives, the assembly of the discretised stiffness matrix ([@sec:building]), the numerical solution of the system of equations ([@sec:solving]) and the computation of the gradient of the solution ([@sec:stress-computation). Yet, the fact that all these a-priori unconnected steps give rise to pretty pictures that resemble reality is still astonishing to me.
 
 
 Again, take all this information as coming from a fellow that has already taken such a journey from college’s pencil and paper to real engineering cases involving complex numerical calculations. And developing, in the meantime, both an actual working finite-element [back-end](https://www.seamplex.com/fino) and [front-end](https://www.caeplex.com) from scratch.^[The dean of the engineering school I attended used to say “It is not the same to read than to write manuals, and we should aim at writing them.”]
@@ -75,13 +75,21 @@ An important part of the analysis that almost always applies to nuclear power pl
 
 ## Fatigue {#sec:fatigue}
 
-Mechanical systems can fail due to a wide variety of reasons. The effect known as fatigue can create,  migrate and grow microscopic cracks at the atomic level, called dislocations. Once these cracks reach a critical size, then the material fails catastrophically even under stresses much lower than the tensile strength limits. There are not complete mechanistic models from first principles which can be used in general situations, and those that exist are very complex and hard to use. Instead, using an experimental approach very much like the [Hooke’s Law](https://en.wikipedia.org/wiki/Hooke's_law) experiment, the stress amplitude of a periodic cycle can be related to the number of cycles where failure by fatigue is expected to occur. For each material, this dependence can be computed using normalised tests and a family of “fatigue curves” like the one depicted in\ [@fig:SN] (also called $S$-$N$ curve) for different temperatures can be obtained.
+Mechanical systems can fail due to a wide variety of reasons. The effect known as fatigue can create,  migrate and grow microscopic cracks at the atomic level, called dislocations. Once these cracks reach a critical size, then the material fails catastrophically even under stresses much lower than the tensile strength limits. There are not complete mechanistic models from first principles which can be used in general situations, and those that exist are very complex and hard to use. There are two main ways to approach practical fatigue assessment problems using experimental data very much like the [Hooke’s Law](https://en.wikipedia.org/wiki/Hooke's_law) experiment:
+
+ 1. stress life, or
+ 2. strain life
+
+The first one is suitable for cases where the loads are nowhere near the yield stress of the material. In those cases where plastic deformation is expected to occur, strain-life methods ought to be employed.
+
+For the case study, as the loads come principally from operational loads, the ASME\ stress-life approach should be used. The stress amplitude of a periodic cycle can be related to the number of cycles where failure by fatigue is expected to occur. For each material, this dependence can be computed using normalised tests and a family of “fatigue curves” like the one depicted in\ [@fig:SN] (also called $S$-$N$ curve) for different temperatures can be obtained.
 
 
 ![A fatigue or $S$-$N$ curve for two steels.](SN.svg){#fig:SN width=95%}
 
 It should be noted that the fatigue curves are obtained in a particular load case, namely purely-periodic and one-dimensional, which cannot be directly generalised to other three-dimensional cases. Also, any real-life case will be subject to a mixture of complex cycles given by a stress time history and not to pure periodic conditions. The application of the curve data implies a set of simplifications and assumptions that are translated into different possible “rules” for composing real-life cycles. There also exist two safety factors which increase the stress amplitude and reduce the number of cycles respectively. All these intermediate steps render the analysis of fatigue into a conservative computation scheme ([@sec:kinds]). Therefore, when a fatigue analysis performed using the fatigue curve method arrives at the conclusion that “fatigue is expected to occur after ten thousand cycles” what it actually means is “we are sure fatigue will not occur before ten thousand cycles, yet it may not occur before one hundred thousand or even more.”  
 
+The usage of $S$-$N$ curves imply that the
 
 # Solid mechanics, or what we are taught at college
 
@@ -694,32 +702,63 @@ then the sums (in plural because there are three eigenvalues) of their eigenvalu
 
 \medskip
 
-The moral of this fable is that if we have a case that is the combination of two other simpler cases (say one with only surface loads and one with only volumetric loads), in general we cannot just add the stresses and expect to obtain the correct principal stresses. We have to solve the full case again (both the surface and the volumetric loads at the same time). In case we are stubborn enough and still want to stick to solving the cases separately, there is a trick we can resort to. Instead of summing stresses, what we can do is to sum displacements which are fully linear. So we might pre-deform case B with the results from case A and then expect to obtain the correct stresses for the combined case. However, this scheme is actually far more complex than just solving the combined case in a single run and it also needs an educated guess and/or trial and error to know at what time the pre-deformation should be applied to take into account the seismic load.
+The moral of this fable is that if we have a case that is the combination of two other simpler cases (say one with only surface loads and one with only volumetric loads), in general we cannot just add the principal stresses (or Von Mises) of two cases and expect to obtain a correct answer. We have to solve the full case again (both the surface and the volumetric loads at the same time). In case we are stubborn enough and still want to stick to solving the cases separately, there is a trick we can resort to. Instead of summing principal stresses, what we can do is to sum either displacements or the individual stress components, which are fully linear. So we might pre-deform (or pre-stress) case B with the results from case A and then expect the FEM program to obtain the correct stresses for the combined case. However, this scheme is actually far more complex than just solving the combined case in a single run and it also needs an educated guess and/or trial and error to know at what time the pre-deformation or pre-stressing should be applied to take into account the seismic load.
 
 
 ## ASME stress linearisation (not linearity!)
 
-After discussing linearity, let us now dig into linearisation. The name is similar but these two animals are very different beasts. We said in\ [@sec:case] that the ASME Boiler and Pressure Vessel Code was born long before modern finite-elements methods were developed and of course being massively available for general engineering analysis (democratised?). Yet the code provides a comprehensive, sound and, more importantly, a widely and commonly-accepted body of knowledge as for the regulatory authorities to require its enforcement to nuclear plant owners. One of the main issues of the ASME code refers to what is known as “membrane” and “bending” stresses, which are defined in [section\ VIII](https://en.wikipedia.org/wiki/ASME_Boiler_and_Pressure_Vessel_Code#ASME_BPVC_Section_VIII_-_Rules_for_Construction_of_Pressure_Vessels) (although widely used in other sections, particularly section\ III) annex 5-A. Briefly, they give the zeroth-order (membrane) and first-order (bending) moments of the stress distribution along a so-called Stress Classification Line or SCL, which should be chosen depending on the type of problem under analysis.
+After discussing linearity, let us now dig into linearisation. The name is similar but these two animals are very different beasts. We said in\ [@sec:case] that the ASME Boiler and Pressure Vessel Code was born long before modern finite-elements methods were developed and of course being massively available for general engineering analysis (democratised?). Yet the code provides a comprehensive, sound and, more importantly, a widely and commonly-accepted body of knowledge as for the regulatory authorities to require its enforcement to nuclear plant owners. One of the main issues of the ASME code refers to what is known as “membrane” and “bending” stresses, which are defined in [section\ VIII](https://en.wikipedia.org/wiki/ASME_Boiler_and_Pressure_Vessel_Code#ASME_BPVC_Section_VIII_-_Rules_for_Construction_of_Pressure_Vessels) (although widely used in other sections, particularly section\ III) annex 5-A. Briefly, they give the zeroth-order (membrane) and first-order (bending) [moments](https://en.wikipedia.org/wiki/Moment_(mathematics)) of the stress distribution along a so-called Stress Classification Line or SCL, which should be chosen depending on the type of problem under analysis.
 
-The computation of these membrane and bending stresses are called [“stress linearisation”](https://www.ramsay-maunder.co.uk/knowledge-base/technical-notes/stress-linearisation-for-practising-engineers/) because (I am guessing) it is like computing the Taylor expansion of a real stress distribution along a line, and retaining the first two terms. That is to say, to obtain a linear approximation.
+The computation of these membrane and bending stresses are called [“stress linearisation”](https://www.ramsay-maunder.co.uk/knowledge-base/technical-notes/stress-linearisation-for-practising-engineers/) because (I am guessing) it is like computing the [Taylor expansion](https://en.wikipedia.org/wiki/Taylor_series) (or for the case, in [Legendre polynomials](https://en.wikipedia.org/wiki/Legendre_polynomials)) of a real stress distribution along a line, and retaining the first two terms. That is to say, to obtain a linear approximation. More (optional) mathematical details below.
 
 dnl **figures**
 
 So what about the SCLs? Well, the ASME standard says that they are lines that go through a wall of the pipe (or vessel or pump, which is what the ASME code is for) from the inside to the outside and ought to be normal to the iso-stress curves. Stop. Picture yourself a stress field, draw the iso-stress curves (those would be the lines that have the same colour in your picture) and then imagine a set of lines that travel in a perpendicular direction to them. Finally, choose the one that seems the prettiest (which most of the time is the one that seems the easiest). There you go! You have an SCL. But there is a catch. So far, we have referred to a generic concept of “stress.” Which of the several stresses out there should you picture? One of the three normals, the three shear, Von\ Mises, Tresca? Well, actually you will have to imagine tensors instead of scalars. And there might not be such a thing as “iso-stress” curves, let alone normal directions. So pick any radial straight line through the pipe wall at a location that seems relevant and now you are done. In our case study, there will be a few different locations around the material interfaces where high stresses due to differential thermal expansion are expected to occur.
 
-If you cannot wait to know, the expression for computing the $i$-$j$-th element of the membrane tensor is
+Now the optional (but recommended) mathematical details. Accoding to ASME\ VIII Div.\ 2 Annex\ 5-A, the expression for computing the $i$-$j$-th element of the membrane tensor is
 
 $$
 \text{M}_{ij} = \frac{1}{t} \cdot \int_0^t \sigma_{ij}(t^\prime) \, dt^\prime
 $$
-where $t$ is the length and $t^\prime$ is the parametrisation variable of the SCL. The other linearised stress, namely the _membrane plus bending_ stress tensor \text{MB} again according to ASME VIII Annex 5-A is
+where $t$ is the length and $t^\prime$ is the parametrisation variable of the SCL. The other linearised stress, namely the _bending_ stress tensor $\text{B}$ is
 
 $$
-\text{MB}_{ij} = \text{M}_{ij} \pm \frac{6}{t^2} \cdot \int_0^t \sigma_{ij}(t^\prime) \cdot \left( \frac{t}{2}-t^\prime\right) \, dt^\prime
+\text{B}_{ij} = \frac{6}{t^2} \cdot \int_0^t \sigma_{ij}(t^\prime) \cdot \left( \frac{t}{2}-t^\prime\right) \, dt^\prime
 $$
-where the sign should be taken such that the resulting combined scalar stress (i.e. Tresca, Von Mises, Principal\ 1, etc.) represents the worst-case scenario.
 
-No need to know or even understand these two integrals which for sure are not introduced to students in regular college courses. But it would be good to, as linearisation is a cornerstone subject for any serious mechanical analysis of pressurised components following the code. So get a copy of ASME sec.\ VIII div.\ 2 annex\ 5-A and then search online for for “stress linearisation” (or “linearization”).
+For the fatigue assessment, it is the sum\ $\text{MB}$ which is taken as a measure of the stress
+
+$$
+\text{MB}_{ij} = \text{M}_{ij} \pm \text{B}_{ij}
+$$
+where the sign should be taken such that the resulting stress intensity (i.e. Tresca, Von Mises, Principal\ 1, etc.) represents the worst-case scenario. Older versions of ASME\ VIII preferred the [Tresca criterion](https://en.wikipedia.org/wiki/Henri_Tresca), while newer versions switched to [Von\ Mises](https://en.wikipedia.org/wiki/Von_Mises_yield_criterion). ASME\ III sticked to Tresca.
+
+In any case, $\text{MB}_{ij}$ should be taken as a measure of the primary stress at the internal surface of the pipe. In effect, let us assume we have a stress distribution\ $\sigma(t^\prime)$ along a certain line of length $t$ such that\ $\sigma(0)$ and $\sigma(t)$ are the stresses at the internal and external surfaces, respectively. The primary stress distribution, i.e. the stress that is not self balancing will be of the simple linear form
+
+$$
+\sigma(t^\prime) = a\cdot t^\prime + b
+$$
+because any other higher-order polynomial term would self-balance.
+The membrane and bending stresses, according to ASME would then be
+
+$$
+\text{M} = \frac{1}{t} \int_{0}^{t} \big[ a\cdot t^\prime + b \big] \, dt^\prime = \frac{1}{2} \cdot a \cdot t + b
+$$
+$$
+\text{B} = \frac{6}{t^2} \int_{0}^{t} \big[ a\cdot t^\prime + b \big] \cdot \left( \frac{t}{2} - t^\prime \right) \, dt^\prime = -\frac{1}{2} \cdot a \cdot t
+$$
+
+We can see that
+$$
+\sigma(0) = b = \text{M} + \text{B}
+$$
+and
+
+$$
+\sigma(t) = a \cdot t + b = \text{M} - \text{B}
+$$
+
+No need to know or even understand these integrals which for sure are not introduced to students in regular college courses. But it would be good to, as linearisation is a cornerstone subject for any serious mechanical analysis of pressurised components following the code. So get a copy of ASME sec.\ VIII div.\ 2 annex\ 5-A and then search online for for “stress linearisation” (or “linearization”).
 
 
 
@@ -773,7 +812,7 @@ You can get both the exponential nature of each added bullet and how easily we c
 Error in the computation of the linearised stresses vs. CPU time needed to solve the problem using FEM
 :::::
 
-An additional note should be added. The FEM solution, which not only gives the nodal displacements but also a method to interpolate these values inside the elements, does not fully satisfy the original equilibrium equations at every point (i.e. the strong formulation). It is an approximation to the solution of the [weak formulation](https://en.wikipedia.org/wiki/Weak_formulation) that is close (measured in the vector space spanned by the shape functions) to the real solution. Mechanically, this means that the FEM solution leads only to nodal equilibrium but not pointwise equilibrium.
+An additional note should be added. The FEM solution, which not only gives the nodal displacements but also a method to interpolate these values inside the elements, does not fully satisfy the original equilibrium equations at every point (i.e. the strong formulation). It is an approximation to the solution of the [weak formulation](https://en.wikipedia.org/wiki/Weak_formulation) that is close (measured in the vector space spanned by the [shape functions](https://www.quora.com/What-is-a-shape-function-in-FEM) to the real solution. Mechanically, this means that the FEM solution leads only to nodal equilibrium but not pointwise equilibrium.
 
 ## Elements, nodes and CPU {#sec:elements-nodes}
 
@@ -808,9 +847,12 @@ The linear FEM problem leads of course of a system of\ $NG$ linear equations, ca
 
 $$K \cdot \mathbf{u} = \mathbf{b}$${#eq:kub}
 
-The objective of the solver is to find the vector\ $\mathbf{u}$ of nodal displacements that satisfy the momentum equilibrium. Luckily (well, not purely by chance but by design) the stiffness matrix is almost empty. It is called a [sparse matrix](https://en.wikipedia.org/wiki/Sparse_matrix) because most of its elements are zero. If it was fully filled, then a problem with just 100k nodes would need more than 700\ Gb of RAM just to store the matrix elements, rendering FEM as practically impossible. And even though the stiffness matrix is sparse, its inverse is not so we cannot solve the elastic problem by “inverting” the matrix. Particular methods to represent and more importantly to solve linear systems involving these kind of matrices have been developed, which are the methods used by finite-element (and the other finite-cousins) programs. These methods are [iterative](https://en.wikipedia.org/wiki/Iterative_method) in nature, meaning that at each step of the iteration the numerical answer improves, i.e. $\left|K \cdot \mathbf{u} - \mathbf{b}\right| \rightarrow 0$. Even though in particular for [Krylov-subspace](https://en.wikipedia.org/wiki/Krylov_subspace) methods, $K \cdot \mathbf{u} - \mathbf{b} = \mathbf{0}$ after\ $NG$ steps, a few dozen of iterations are usually enough to assume that the problem is effectively solved (i.e. the residual is less than a certain threshold).
+The objective of the solver is to find the vector\ $\mathbf{u}$ of nodal displacements that satisfy the momentum equilibrium. Luckily (well, not purely by chance but by design) the stiffness matrix is almost empty. It is called a [sparse matrix](https://en.wikipedia.org/wiki/Sparse_matrix) because most of its elements are zero. If it was fully filled, then a problem with just 100k nodes would need more than 700\ Gb of RAM just to store the matrix elements, rendering FEM as practically impossible. And even though the stiffness matrix is sparse, its inverse is not so we cannot solve the elastic problem by “inverting” the matrix. Particular methods to represent and more importantly to solve linear systems involving these kind of matrices have been developed, which are the methods used by finite-element (and the other finite-cousins) programs. In general there are two approaches
+
+ * [direct](https://en.wikipedia.org/wiki/Frontal_solver), where a rather complex algorithm for building partial decompositions ([LU](https://en.wikipedia.org/wiki/LU_decomposition), [QR](https://en.wikipedia.org/wiki/QR_decomposition), [Cholesky](https://en.wikipedia.org/wiki/Cholesky_decomposition), etc.) are used in order to avoid needing aforementioned the 700\ Gb of RAM, or
+ * [iterative](https://en.wikipedia.org/wiki/Iterative_method), meaning that at each step of the iteration the numerical answer improves, i.e. $\left|K \cdot \mathbf{u} - \mathbf{b}\right| \rightarrow 0$. Even though in particular for [Krylov-subspace](https://en.wikipedia.org/wiki/Krylov_subspace) methods, $K \cdot \mathbf{u} - \mathbf{b} = \mathbf{0}$ after\ $NG$ steps, a few dozen of iterations are usually enough to assume that the problem is effectively solved (i.e. the residual is less than a certain threshold).
 
-So the question is... how hard is to solve a sparse linear problem? Well, the number of iterations needed to attain convergence depends on the [spectral radius](https://en.wikipedia.org/wiki/Spectral_radius) of the stiffness matrix\ $K$, which in turns depend on the elements themselves, which depend mainly on the parameters of the elastic problem. But it also depends on how sparse\ $K$ is, which changes with the element topology and order. [@Fig:test] shows the structure of two stiffness matrices for the same linear elastic problem discretised using exactly the same number of nodes but using linear and quadratic elements respectively. The matrices have the same size (because the number of nodes is the same) but the former is more sparse than the latter. Hence, it would be harder (in computational terms of CPU and RAM) to solve the second-order case.
+So the question is... how hard is to solve a sparse linear problem? Well, the number of iterations and the complexity of the partial decompositions needed to attain convergence depends on the [spectral radius](https://en.wikipedia.org/wiki/Spectral_radius) of the stiffness matrix\ $K$, which in turns depend on the elements themselves, which depend mainly on the parameters of the elastic problem. But it also depends on how sparse\ $K$ is, which changes with the element topology and order. [@Fig:test] shows the structure of two stiffness matrices for the same linear elastic problem discretised using exactly the same number of nodes but using linear and quadratic elements respectively. The matrices have the same size (because the number of nodes is the same) but the former is more sparse than the latter. Hence, it would be harder (in computational terms of CPU and RAM) to solve the second-order case.
 
 ::::: {#fig:test}
 ![42k first-order elements](test1.png){#fig:test1 width=48%}\  
@@ -823,7 +865,7 @@ In a similar way, different types of elements will give rise to different sparsi
 
 ### Stress computation  {#sec:stress-computation}
 
-In the [displacement-based formulation](http://web.mit.edu/kjb/www/Books/FEP_2nd_Edition_4th_Printing.pdf), the solver finds the displacements\ $\mathbf{u}(\mathbf{x})$ that satisfy [@eq:kub], which are the principal unknowns. But from\ [@sec:tensor] we know that we actually have solved the problem after we have the stress tensors at every location\ $\mathbf{x}$, which are the secondary unknowns. So the FEM program has to compute the stresses from the displacements using the materials’ [strain-stress constitutive equations](https://en.wikipedia.org/wiki/Constitutive_equation#Stress_and_strain) which involve the Young Modulus\ $E$, the Poisson ratio\ $\nu$ and the spatial derivatives of the displacements\ $\mathbf{u}=[u,v,w]$. This sounds easy, as we (well, the solver) knows what the shape functions are for each element and then it is a matter of computing nine derivatives (that of the three displacements in each of the three directions) and multiplying by something involving\ $E$ and\ $\nu$. Yes, but there is a catch. As the displacements\ $u$, $v$ and\ $w$ are computed at the nodes, we would like to also have the stresses at the nodes. However,
+In the [displacement-based formulation](http://web.mit.edu/kjb/www/Books/FEP_2nd_Edition_4th_Printing.pdf), the solver finds the displacements\ $\mathbf{u}(\mathbf{x})$ that satisfy [@eq:kub], which are the principal unknowns. But from\ [@sec:tensor] we know that we actually have solved the problem after we have the stress tensors at every location\ $\mathbf{x}$, which are the secondary unknowns. So the FEM program has to compute the stresses from the displacements. It first computes the strain tensor, which is composed of the nine partial derivatives of the three displacements with respect to the three coordinates. Then it computes the stress tensor (atready introduced in\ [@sec:tensor]) using the materials’ [strain-stress constitutive equations](https://en.wikipedia.org/wiki/Constitutive_equation#Stress_and_strain) which involve the Young Modulus\ $E$, the Poisson ratio\ $\nu$ and the spatial derivatives of the displacements\ $\mathbf{u}=[u,v,w]$. This sounds easy, as we (well, the solver) knows what the shape functions are for each element and then it is a matter of computing nine derivatives and multiplying by something involving\ $E$ and\ $\nu$. Yes, but there is a catch. As the displacements\ $u$, $v$ and\ $w$ are computed at the nodes, we would like to also have the stresses at the nodes. However,
 
  i. the displacements\ $\mathbf{u}(\mathbf{x})$ are not differentiable at the nodes, and
  ii. if the node belongs to a material interface, neither\ $E$ nor\ $\nu$ are defined.