cosas

nafems4.md

@@ -9,7 +9,7 @@ We will be swinging back and forth between a case study about fatigue analysis i
 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 preconditioner for the inversion of the stiffness matrix passing through Sobolev spaces and the grid 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 assembly, the numerical solution of the system of equations and the computation of the gradient of the solution. Yet, the fact that all these a-priori unconnected steps once gets a pretty picture that resembles 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.
+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, we should aim at writing.”]
 
 
 ## Tips and tricks
@@ -98,10 +98,10 @@ It looks (and works) like a regular $3 \times 3$ matrix. Some brief comments abo
 
 What does this all have to do with mechanical engineering? Well, once we know what the stress tensor is for every point of a solid, in order to obtain the internal forces per unit area acting in a plane passing through that point and with a normal given by the direction $\vec{n}$, all we have to do is “project” the stress tensor through $\vec{n}$. In plain simple words:
 
- * If you can compute the stress tensor at each point of our geometry, then... Congratulations! You have solved the solid mechanics problem.
+> If you can compute the stress tensor at each point of our geometry, then... Congratulations! You have solved the solid mechanics problem.
 
 
-## An infinitely-long pressurised pipe
+## An infinitely-long pressurised pipe {#sec:infinite-pipe}
 
 Let us proceed to a our second step, and consider an infinite pipe subject to uniform internal pressure. Actually, we are going to solve the mechanical problem on an infinite hollow cylinder, which looks like pipe. This case is usually tackled in college courses, and chances are you already solved it. Actually, the first (and simpler) problem is the “thin cylinder problem.” Then, the “thick cylinder problem” is introduced, which is slightly more complex. Nevertheless, it has an analytical solution which is derived [here](https://www.seamplex.com/fino/doc/pipe-linearized/). For the present case, Let us consider an infinite pipe (i.e. a hollow 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$.
 
@@ -165,7 +165,7 @@ That is all what we can say about an infinite pipe with uniform material propert
 
 Besides infinite pipes (both thin and thick), spheres and a couple of other geometries, there are not other cases for which we can obtain analytical expressions for the elements of the stress tensor. To get results for a solid with real engineering interest, we need to use numerical methods to solve the equilibrium equations. It is not that the equations are hard _per se_. It is that the mechanical parts we engineers like to design (which are of course better than cylinders and spheres) are so intricate that render simple equations into monsters which are unsolvable with pencil and paper. Hence, finite elements enter into the scene.
 
-## The name of the game
+## The name of the game {#sec:formulations}
 
 But before turning our attention into finite elements (and leaving college, at least undergraduate) it is worth some time to think about other alternatives. Are we sure we are tackling your problems in the best possible way? I mean, not just engineering problems. Do we take a break, step back for a while and see the whole picture looking at all the alternatives so we can choose the best cost-effective one?
 
@@ -479,7 +479,7 @@ The ASME code says that these accelerations (depicted in [@fig:acceleration]) ar
 
 Even though we did not yet discuss it in detail, we want to solve an elastic problem subject to an internal pressure condition, with a non-uniform temperature distribution that leads to both thermal stresses and variations in the mechanical properties of the materials. And as if this was not enough, we want to add at some instants a statically-equivalent distributed load that comes from a design earthquake. This last point means that at the transient instant where the stresses (from the fatigue’s point of view) are maximum we have to add the distributed loads that we computed from the seismic spectra to the other thermal and pressure loads. But we have a linear elastic problem (well, we still do not have it but we will in\ [@sec:break]), so we might be tempted exploit the problem’s linearity and compute all the effects separately and them sum them up to obtain the whole combination. We may thus compute just the stresses due to the seismic loads and then add them up to the stresses of any instant of the transient, in particular to the one with the highest ones. After all, in linear problems the result of the sum of two cases is the results of the sum of the cases, right? Wrong.
 
-Let us jump out of our nuclear piping problem and step back into the general finite-element theory ground for a moment (remember we were going to jump back and forth). Assume you want to know how much your dog weights. One thing you can do is weight yourself (let us say you weight 81.2\ kg), then grab your dog and weight yourself and your dog (let us say you and your dog weight 87.3\ kg). Do you swear your dog weights 6.1\ kg plus/minus the scale’s uncertainty? I can tell you that the weight of two individual protons and two individual neutrons in not the same as the weight of an\ [$alpha$ particle](Alpha_particle). Will not there be a master-pet interaction that renders the weighting problem non-linear?
+Let us jump out of our nuclear piping problem and step back into the general finite-element theory ground for a moment (remember we were going to jump back and forth). Assume you want to know how much your dog weights. One thing you can do is weight yourself (let us say you weight 81.2\ kg), then grab your dog and weight yourself and your dog (let us say you and your dog weight 87.3\ kg). Do you swear your dog weights 6.1\ kg plus/minus the scale’s uncertainty? I can tell you that the weight of two individual protons and two individual neutrons in not the same as the weight of an\ [$\alpha$ particle](Alpha_particle). Will not there be a master-pet interaction that renders the weighting problem non-linear?
 
 Let us both (i.e. you and me) make an experiment. Grab a FEM program of your choice (mine is [CAEplex](https://caeplex.com)) and load a 1mm $\times$ 1mm $\times$ 1mm cube. Set any values for the Young Modulus and Poisson ratio as you want. I chose\ $E=200$MPa and\ $\nu=0.28$. Restrict the three faces pointing to the negative axes to their planes, i.e.
 
@@ -648,7 +648,33 @@ $$
 then the sums of their eigenvalues are equal to the eigenvalues of the sums. In order for this to happen, both\ $A$ and\ $B$ need to be diagonalizable using the same basis. That is to say, the diagonalizing matrix\ $P$ such that $P^{-1} A P$ is diagonal is the same that renders\ $P^{-1} B P$ also diagonal. What does this mecanically mean? Well, if case\ A has loads that are somehow “independent” from the ones in case\ B, then the principal stresses of the combination might be equal to the sum of the individual principal stresses. A notable case is for example a beam that is loaded vertically in case\ A and horizontally in case\ B. I dare you to grab your FEM program one more time, run a test and picture the eigenvalues and eigenvectors of the three cases in your head.
 
 
-# The infinite pipe revisited after college
+## ASME stress linearisation (not linearity!)
+
+After discussing linearity, let us now dig into linearisation. The name is similar but these two guys 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 provide 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 (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.
+
+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 first two terms of 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.
+
+**figures**
+
+So what about the SCLs? Well, the ASME standard says that they are lines that go through a wall of the ipe (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 color 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
+
+$$
+\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 a parametrization 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
+
+$$
+\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
+$$
+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”).
+
+
+
+# The infinite pipe revisited after college {#sec:infinite-pipe-fem}
 
 Let us now make a (tiny) step from the general and almost philosophical subject from the last section down to the particular case, and reconsider the infinite pressurized pipe once again. It is time to solve the problem with a computer using finite elements and obtain some funny colored pictures instead of just equations.
 
@@ -697,35 +723,66 @@ 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
 :::::
 
-The last two bullets lead to an issue that has come many times when discussing the issue of convergence with respect to the mesh size with other colleagues. There apparently exists a common misunderstanding that the number of _elements_ is the main parameter that defines how complex a FEM model is.
+## Elements, nodes and CPU
 
-meshing, building, solving, stress
+The last two bullets above lead to an issue that has come many times when discussing the issue of convergence with respect to the mesh size with other colleagues. There apparently exists a common misunderstanding that the number of elements is the main parameter that defines how complex a FEM model is. This is strange, because even in college we are taught that the most important parameter is the size of the stiffness matrix, which is three times (for 3D problems) the number of _nodes_.
 
-matrices test
+Let us pretend we have to compare two different FEM programs so we solve the same problem in each and see what the results are. I have seen many times the following situation: the user loads the same geometry in both programs, meshes them so that the number of elements is more or less them same (because she wants to be “fair”) and solve the problem. Voilà! It turns out that the first program defaults to first-order elements and the second one to second-order elements. So if the first one takes one minute to obtain a solution, the second one should take nearly four minutes. How come that is a fair comparison? Or it might be the case that one program uses tetrahedra while the other one hexahedra. Or any other combination. In general, there is no single problem parameter that can be fixed to have a “fair” comparison, but if there were it would definetely be the number of nodes rather than the number of elements. Let us see why.
 
-## ASME stress linearisation (not linearity!)
+Fire up your imagination again and make a [thought experiment](https://en.wikipedia.org/wiki/Thought_experiment) in which you have to compare say a traditional FEM approach with a new radical formulation that a crazy mathematician from central Asia came up with claiming it is a far superior theory than our beloved finite elements (or for the case, any other formulation from\ [@sec:formulations]). How can we tell if the guy is really a genius or purely nuts? Well, we could solve a problem which we can compute the analytical solution (for example the infinite pipe from\ @[sec:infinite-pipe]) first with the traditional method (@[sec:infinite-pipe-fem]) and then with the program which uses the new formulation. Say the traditional FEM gives an error between 1% and 5% running in a few seconds depending on the mesh size. The new program from the crazy guy takes no input parameters and gives an error of 0.001%, but it takes one week of computation to produce a result. Would you say that the new radical formulation is really far superior?
 
-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 provide 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 (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.
+What I would do is to run a FEM program that takes also one week to compute, and only then compare the errors. So that is why\ [@fig:error-vs-cpu] uses the CPU time in the abscissae rather than the number of elements to compare first and second-order formulations.
 
-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 first two terms of 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.
+To fix ideas, let us stick to a linear elastic FEM problem. The CPU time needed to completely solve such a problem can be divided into four steps:
 
-**figures**
+ 1. meshing the continuous geometry,
+ 2. building the stiffness matrix,
+ 3. solving the equations to obtain the displacements, and
+ 4. computing the stress from the displacements.
 
-So what about the SCLs? Well, the ASME standard says that they are lines that go through a wall of the ipe (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 color 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.
+### Meshing
 
-If you cannot wait to know, the expression for computing the $i$-$j$-th element of the membrane tensor is
+The effort needed to compute a discretization of a continuous domain depends on the meshing algorithm. But nearly all meshers first put nodes on the edges (1D), then on the surfaces (2D) and finally on the volumes. Afterwards, they join the nodes to create the elements. Depending on the topology (i.e. tetrahedra, hexahedra, pyramids, etc) and the order (i.e. linear, quadratic, etc.) this last step will vary, but the main driver here is the number of nodes. Try measuring the time needed to obtain grids of different sizes and kinds with your mesher.
 
-$$
-\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 a parametrization 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
+### Building
 
-$$
-\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
-$$
-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.
+The stiffness matrix is a square matrix that has\ $NG$ rows and\ $NG$ columns where $N$ is the number of nodes and $G$ is the number of degrees of freedom per node, which for three-dimensional problems is $G=3$. But even though FEM problems have to build a $NG\times NG$ matrix, they usually sweep through elements rather than through nodes, and then scatter the elements of the elemental matrices to the global stiffness matrix. So the effort needed here depends again on how the solver is programmed, but it is a combination of the number of elements and the number of nodes.
 
-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 a fixed number of nodes, first-order grids have far more elements than second-order grids because in the first case each node has to be a vertex while in the latter half will be vertexes and half will be located at the edges (think! FIGURE?). So the sweep is larger for linear grids. But the effort needed to integrate properties using quadratic shape functions is greater than for the linear case, so these two effects almost cancel out.
+
+### Solving
+
+The linear FEM problem leads of course of a system of\ $NG$ linear equations, casted in matrix form by the stiffness matrix\ $K$ and a right-hand size vector\ $\vec{b}$ containing the loads (both volumetric and the ones at the surfaces from the boundary conditions):
+
+$$K \cdot \vec{u} = \vec{b}$$
+
+The objective of the solver is to find the vector\ $u$ of nodal displacements that satisfies the momentum equilibriums.
+
+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 \vec{u} - \vec{b}\right| \rightarrow 0$. Even though in particular for [Krylov-subspace](https://en.wikipedia.org/wiki/Krylov_subspace) methods, $K \cdot \vec{u} - \vec{b} = \vec{0}$ after\ $NG$ steps, a few dozen of iterations are usually enough to assume that the problem is effectively solve (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 discretized using exactly the same number of nodes with first and second-order elements. 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%}
+![15k second-order elements](test2.png){#fig:test2 width=48%}
+
+Structure of the stiffness matrices for the same FEM problem with 10k nodes. Blue (red) are positive (negative) elements.
+:::::
+
+In a similar way, different types of elements will give rise to different sparsity patterns which change the effort needed to solve the problem. In any case, the base parameter that controls the problem size and thus provides a basic indicator of the level of difficulty the problem poses is the number of nodes. And again, not the number of elements.
+
+### Stress computation
+
+In the most common FEM formulation, the solver just “solves” for the displacements\ $\vec{u}(\vec{x})$ which are the principal unknowns. But from\ [@sec:tensor] we know that we actually “solve” the problem when we have the stress tensors at every location\ $\vec{x}$, which are the secondary unknowns. So the FEM program has to compute the stresses from the displacements using the materials’ constitutive equations which involve the Young Modulus\ $E$, the Poisson ratio\ $\nu$ and the spatial derivatives of the displacements\ $\vec{u}$. 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$ are computed at the nodes, we would like to also have the stresses at the nodes. But
+
+ i. the displacements\ $\vec{u}(\vec{x})$ are not differentiable at the nodes, and
+ ii. if the node belongs to a material interface, the\ $E$ and\ $\nu$ stuff is not defined.
+
+Let us focus on the first item and leave the second one for a separate discussion in\ [@sec:two-materials]. The finite-element method computes the principal unknowns at the nodes and then says “interpolate the nodal values inside each element using its shape functions.” It sounds (and it is) great, but a node belongs to more than one element. So the slope of the interpolation when we move from the node into one element might (and it never is) the same as when we move from the same node into another element. Neither in linear, nor quadratic nor higher-order elements. So what is the derivative of the displacement\ $v$ in the\ $y$ direction with respect to the $z$ coordinate? One answer would be the average of the derivatives computed in each of the elements that share a common node. But what if one of the elements is very small? Or it has a very bad quality (i.e. it is deformed in one direction) its derivatives cannot be trusted? Should the “average” be weighted? How?
+
+**FIGURA 1D**
+
+inside elements at gauss points
 
 
 
@@ -761,6 +818,9 @@ Back in College, we all learned how to solve engineering problems. But there is
  * play with your favourite FEM solver (mine is <https://caeplex.com>) solving simple cases, trying to predict the results and picturing the stress tensor and its eigenvectors in your imagination
  * grab any stress distribution from any of your FEM projects, compute the iso-stress curves and the draw normal lines to them to get acquainted with SCLs
  * first search online for “stress linearisation” (or “linearization” if you want) and then get a copy of ASME\ VIII Div\ 2 Annex 5-A
+ * try to show that even though the principal stresses are not linear with respect to summation, they are linear with respect to multiplication (with pencil and paper)
+ * learn by heart: the complexity of a FEM problem is given mainly but the number of nodes, not by the number of elements
+ * measure the time needed to generate grids of different sizes and kinds with your favourite mesher