gtheler
/
nafems4
Seguir 1
Destacar 0
Fork 0
Código Incidencias 0 Pull Requests 0 Lanzamientos 0 Wiki Actividad
Explorar el Código

cosas

tags/CANDIDATE
gtheler hace 7 años
padre
ba15326782
commit
2216e66e1b
Se han modificado 1 ficheros con 68 adiciones y 45 borrados
Dividir vista Mostrar estadísticas de diff
  1. +68
    -45
      nafems4.md

+ 68
- 45
nafems4.md Ver fichero

@@ -2,8 +2,7 @@

First of all, please take this text as a written chat between you an me, i.e. an average engineer that have already taken the journey from college to performing actual engineering using finite element analysis and has something to say about it. Picture yourself in a coffee bar, talking and discussing concepts and ideas with me. Maybe needing to go to a blackboard (or notepad?). Even using a tablet to illustrate some three-dimensional results. But always as a chat between colleagues.

Please also note that I am not a mechanical engineer, although I shared many undergraduate courses with some of them. I am a nuclear engineer with a strong background on mathematics and computer programming. I went to college between 2002 and 2008. Probably a lot of things have changed since then---at least that is what these “millenials” guys and girls seem to be boasting about---but chances are we all studied solid mechanics and heat transfer with a teacher using a piece of chalk on a blackboard and students writing down notes with pencils on paper sheets. And there is really not much that one can do with pencil and paper regarding mechanical analysis. Any actual case worth the time of an engineer need to be more complex than an ideal canonical case with analytical solution.

Please also note that I am not a mechanical engineer, although I shared many undergraduate courses with some of them. I am a nuclear engineer with a strong background on mathematics and computer programming. I went to college between 2002 and 2008. Probably a lot of things have changed since then---at least that is what these “millennial” guys and girls seem to be boasting about---but chances are we all studied solid mechanics and heat transfer with a teacher using a piece of chalk on a blackboard and students writing down notes with pencils on paper sheets. And there is really not much that one can do with pencil and paper regarding mechanical analysis. Any actual case worth the time of an engineer need to be more complex than an ideal canonical case with analytical solution.

::::: {#fig:pendulum}
![Simple pendulum](simple.svg){#fig:simple width=35%}\
@@ -20,7 +19,7 @@ An infinitely-long pressurised thick pipe as taught in college and a section of
:::::


Whether you are a student or a seasoned engineer with many years of experience, you might recall from first year physics courses the introduction of the [simple pendulum](https://en.wikipedia.org/wiki/Pendulum) as case study\ ([@fig:simple]). You learned that the period does not depend on the hanging mass because the weight and the inertia exactly canceled each other. Also, that Galileo said (and [Newton proved](https://www.seamplex.com/wasora/realbook/real-012-mechanics.html)) that for small oscillations the period does not even depend on the amplitude. Someone showed you why it worked this way: because if\ $\sin \theta \approx \theta$ then the motion equations converge to an [harmonic oscillator](https://en.wikipedia.org/wiki/Harmonic_oscillator). It might have been a difficult subject for you back in those days when you were learning physics and calculus at the same time. You might later study the [Lagrangian](https://en.wikipedia.org/wiki/Lagrangian_mechanics) and even the [Hamiltonian](https://en.wikipedia.org/wiki/Hamiltonian_mechanics) formulations, added a [parametric excitation](https://en.wikipedia.org/wiki/Parametric_oscillator) and studied the [chaotic double pendulum](https://www.seamplex.com/wasora/realbook/real-017-double-pendulum.html). But it was probably after college, say when you took your first son to a swing on a windy day\ ([@fig:hamaca]), that you were faced with a real pendulum worth your full attention. Ok, this is my personal story but could easily be yours as well. My point is that the very same distance between what I imagined as a student what studying a pendulum was and what I saw that day at the swing (namely that the period does depend on the hanging mass) is the same distance between the mechanical problems studied in college and the actual cases encountered during a professional engineer’s lifetime. In this regard, I am referring only to technical issues. The part of dealing with clients, colleagues, bosses, etc. which is definitely not taught in engineering schools (you can get a heads up in business schools, but again it would be a theoretical pendulum) is way beyond the scope of both this article and my own capacities.
Whether you are a student or a seasoned engineer with many years of experience, you might recall from first year physics courses the introduction of the [simple pendulum](https://en.wikipedia.org/wiki/Pendulum) as case study\ ([@fig:simple]). You learned that the period does not depend on the hanging mass because the weight and the inertia exactly cancelled each other. Also, that Galileo said (and [Newton proved](https://www.seamplex.com/wasora/realbook/real-012-mechanics.html)) that for small oscillations the period does not even depend on the amplitude. Someone showed you why it worked this way: because if\ $\sin \theta \approx \theta$ then the motion equations converge to an [harmonic oscillator](https://en.wikipedia.org/wiki/Harmonic_oscillator). It might have been a difficult subject for you back in those days when you were learning physics and calculus at the same time. You might later study the [Lagrangian](https://en.wikipedia.org/wiki/Lagrangian_mechanics) and even the [Hamiltonian](https://en.wikipedia.org/wiki/Hamiltonian_mechanics) formulations, added a [parametric excitation](https://en.wikipedia.org/wiki/Parametric_oscillator) and studied the [chaotic double pendulum](https://www.seamplex.com/wasora/realbook/real-017-double-pendulum.html). But it was probably after college, say when you took your first son to a swing on a windy day\ ([@fig:hamaca]), that you were faced with a real pendulum worth your full attention. Ok, this is my personal story but could easily be yours as well. My point is that the very same distance between what I imagined as a student what studying a pendulum was and what I saw that day at the swing (namely that the period does depend on the hanging mass) is the same distance between the mechanical problems studied in college and the actual cases encountered during a professional engineer’s lifetime. In this regard, I am referring only to technical issues. The part of dealing with clients, colleagues, bosses, etc. which is definitely not taught in engineering schools (you can get a heads up in business schools, but again it would be a theoretical pendulum) is way beyond the scope of both this article and my own capacities.

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.

@@ -38,11 +37,11 @@ There are some useful tricks that come handy when trying to solve a mechanical p
One of the most important ones is using your _imagination_. You will need a lot of imagination to “see” what it is actually going on when analysing an engineering problem. This skill comes from my background in nuclear engineering where I had not choice but to imagine a [positron-electron annihilation](https://en.wikipedia.org/wiki/Electron%E2%80%93positron_annihilation) or an [Spontaneous fission](https://en.wikipedia.org/wiki/Spontaneous_fission). But in mechanical engineering, it is likewise important to be able to imagine how the loads “press” one element with the other, how the material reacts depending on its properties, how the nodal displacements generate stresses (both normal and shear), how results converge, etc. And what these results actually mean besides the pretty-coloured figures (a former boss once told me “I need the CFD” when I handed in some results. I replied that I did not do computational fluid-dynamics but computed the neutron flux kinetics within a nuclear reactor core. He joked “I know, what I need are the _Colors For Directors_, those pretty coloured figures along with your actual results.”).
This journey will definitely need your imagination. We will see equations, numbers, plots, schematics, 3D geometries, interactive 3D views, etc. Still, when the theory says “thermal expansion produces linear stresses” you have to picture in your head three little arrows pulling away from the same point in three directions, or whatever mental picture you have about what you understand are thermally-induced stresses. What comes to your mind when someone says that out of the nine elements of the stress tensors there are only six that are independent? Whatever it is, try to practice that kind of graphical thoughts with every concept.

Another heads up is that we will dig into some math. Probably it would be be simple and you would deal with it very easily. But probably you do not like equations. No problem! Just ignore them for now. Read the text skipping them, it should work.
Łukasz Skotny says [you do not need to know math to perform finite-element analysis](https://enterfea.com/math-behind-fea/). And he is right, in the sense that you do not need to know thermodynamics to drive a car. It is fine to ignore math for now.
Another heads up is that we will dig into some mathematics. Probably it would be be simple and you would deal with it very easily. But probably you do not like equations. No problem! Just ignore them for now. Read the text skipping them, it should work.
Łukasz Skotny says [you do not need to know maths to perform finite-element analysis](https://enterfea.com/math-behind-fea/). And he is right, in the sense that you do not need to know thermodynamics to drive a car. It is fine to ignore maths for now.
But, eventually, a time will come in which it cannot (or should not) be avoided. If you want to go to space, you will definitely have to learn thermodynamics.

So here comes another experience tip: do not fear mathematics. Even more, keep exercising. You have used differences of squares in high school. You know (or at least knew) how to integrate by parts. Remember what Laplace transforms are used for? Once in a while, perform a division of polynomials using [Ruffini’s rule](https://en.wikipedia.org/wiki/Ruffini's_rule). Or compute the second derivative of the quotient of two functions. Whatever. It should be like doing crosswords on the newspaper. Grab those old physics college books and read the exercises at the end of each chapter. It will pay off later on.
So here comes another experience tip: do not fear equations. Even more, keep exercising. You have used differences of squares in high school. You know (or at least knew) how to integrate by parts. Remember what Laplace transforms are used for? Once in a while, perform a division of polynomials using [Ruffini’s rule](https://en.wikipedia.org/wiki/Ruffini's_rule). Or compute the second derivative of the quotient of two functions. Whatever. It should be like doing crosswords on the newspaper. Grab those old physics college books and read the exercises at the end of each chapter. It will pay off later on.

Throught the text I will be referring to “your favourite FEM program.” I bet you do have one. We will be using it to perform some tests and play a little bit. And we will use it to think about what it means to use a FEM program to generate results that will eventually end up in a written project with your signature. Keep that in mind.

@@ -66,7 +65,7 @@ In each of the countries that have at least one nuclear power plant there exists

How come that pipes are subject to fatigue? Well, on the one hand and without getting into many technical details, the most common nuclear reactor design uses liquid water as coolant and moderator. On the other hand, nuclear power plants cannot by-pass the thermodynamics of the Carnot cycle, and in order to maximise the efficiency of the conversion between the energy stored in the uranium nuclei into electricity they need to reach temperatures as high as possible. So, if we want to have liquid water in the core as hot as possible, we need to increase the pressure. The limiting temperature and pressure are given by the [critical point of water](https://en.wikipedia.org/wiki/Critical_point_(thermodynamics)), which is around 374ºC and 22\ MPa. It is therefore expected to have temperature and pressures near those values in many systems of the plant, especially in the primary circuit those that directly interact with it, such as pressure and inventory control system, decay power removal system, feedwater supply system, emergency core-cooling system, etc.

Nuclear power plants are not always working at 100% power. They need to be maintained and refueled, they may undergo operational transients, they might operates at a lower power due to load following conditions, etc. These transient cases involved changes both in temperatures and in pressures that the pipes are subject to, which in turn give rise to changes in the stresses within the pipes. As the transients are postulated to occur conservatively cyclically during a number of times during the life-time of the plant (plus its extension period), mechanical fatigue in these piping systems arise especially at the interfaces between materials with different thermal expansion coefficients.
Nuclear power plants are not always working at 100% power. They need to be maintained and refuelled, they may undergo operational transients, they might operates at a lower power due to load following conditions, etc. These transient cases involved changes both in temperatures and in pressures that the pipes are subject to, which in turn give rise to changes in the stresses within the pipes. As the transients are postulated to occur conservatively cyclically during a number of times during the life-time of the plant (plus its extension period), mechanical fatigue in these piping systems arise especially at the interfaces between materials with different thermal expansion coefficients.

An important part of the analysis that almost always applies to nuclear power plants but usually also to other installations is the consideration of a possible seismic event. Given a postulated design earthquake, the civil structures that hold the pipes have a response spectra for each floor level. One has to combine these spectra with the natural oscillation modes and frequencies of the piping system using one of several available methodologies such as the “Square Root of Sum of Squares” or SRSS method. This way, an static-equivalent distributed internal load can be computed and then applied as extra loads conservatively (see [@sec:kinds]) at the moment of highest mechanical demand.

@@ -88,11 +87,11 @@ So, let us start our journey. Our starting place: undergraduate solid mechanics
2. in order for a solid not to move, the sum of all the forces ought to be equal to zero, and
3. for every external load there exists an internal reaction with the same magnitude but opposite direction.

We have to accept that there is certain intellectual beauty when complex stuff can be expressed in simple term. Yet, from now on, everything can be complicated at will. We can take the mathematical path like D’Alembert and his virtual displacements ideas (in his mechanical treatise, D’Alembert brags that he does not need to use a single figure throughout the book). Or we can go graphical following Cullman. Or whatever other logic reasoning to end up with a set of actual equations which we need to solve in order to obtain engineering results.
We have to accept that there is certain intellectual beauty when complex stuff can be expressed in simple term. Yet, from now on, everything can be complicated at will. We can take the mathematical path like [D’Alembert](https://en.wikipedia.org/wiki/Jean_le_Rond_d'Alembert) and his virtual displacements ideas (in his mechanical treatise, he brags that he does not need to use a single figure throughout the book). Or we can go graphical following [Culmann](https://en.wikipedia.org/wiki/Carl_Culmann). Or whatever other logic reasoning to end up with a set of actual equations which we need to solve in order to obtain engineering results.

## The stress tensor {#sec:tensor}

In any case, what we should understand (and imagine) is that external forces lead to internal stresses. And in any three-dimensional body subject to such external loads, the best way to represent internal stresses is through a $3 \times 3$ _stress tensor_. This is the first point in which we should not fear math. Trust me, it will pay back later on.
In any case, what we should understand (and imagine) is that external forces lead to internal stresses. And in any three-dimensional body subject to such external loads, the best way to represent internal stresses is through a $3 \times 3$ _stress tensor_. This is the first point in which we should not fear mathematics. Trust me, it will pay back later on.

Does the term _tensor_ scare you? It should not. A tensor is just a slightly more complex vector, and I assume you are not afraid of vectors. If you recall, a vector somehow generalises the idea of a scalar in the following sense: a given vector $\vec{v}$ can be projected into any direction $\vec{n}$ to obtain a scalar $p$. We call this scalar $p$ the “projection” of the vector $\vec{v}$ in the direction $\vec{n}$. Well, a tensor can be also projected into any direction $\vec{n}$. The difference is that instead of a scalar, a vector is now obtained.

@@ -129,7 +128,7 @@ Let us proceed to a our second step, and consider an infinite pipe subject to un

### Displacements

Remember that when any solid body is subject to external forces, it has to reach in such a way to satisfy the equilibrium conditions. The way solids do this is by deforming a little bit in such a way that the whole body acts as a compressed (or elongated) spring balancing the load. So it is worth to ask how a pressurized pipe deforms to counteract the internal pressure\ $p$.
Remember that when any solid body is subject to external forces, it has to reach in such a way to satisfy the equilibrium conditions. The way solids do this is by deforming a little bit in such a way that the whole body acts as a compressed (or elongated) spring balancing the load. So it is worth to ask how a pressurised pipe deforms to counteract the internal pressure\ $p$.

* There are no longitudinal displacements\ $u_l$ because the pipe is infinite in the axial direction.
* There are no azimuthal displacements\ $u_\theta$ because the pipe is fully symmetric around the axis.
@@ -197,12 +196,12 @@ There are literally dozens of ways to numerically solve the equilibrium equation
2. Finite volumes
3. Finite elements

Each of these methods (also called schemes) have of course their own features, pros and cons. They all exploit the fact that the equations are easy to solve in simple geometries (say a cube). Then the actual geometry is divided into a yuxtaposition of these cubes, the equations are solved in each one and then a global solution is obtained by sewing the little simple solutions one to another. The process of dividing the original domain into simple geometries is called _discretization_, and the resulting collection of these simple geometries is called a mesh or grid. They are composed of volumes, called cells (or elements) and vertices called nodes. Now, grids can be either
Each of these methods (also called schemes) have of course their own features, pros and cons. They all exploit the fact that the equations are easy to solve in simple geometries (say a cube). Then the actual geometry is divided into a juxtaposition of these cubes, the equations are solved in each one and then a global solution is obtained by sewing the little simple solutions one to another. The process of dividing the original domain into simple geometries is called _discretization_, and the resulting collection of these simple geometries is called a mesh or grid. They are composed of volumes, called cells (or elements) and vertices called nodes. Now, grids can be either

a. structured, or
b. unstructured

Figure [-@fig:grids] illustrates how the same domain can be discretized using these two kind of grids. In the first case, we could identify any single cell by using just two indexes. We could even tell which nodes define each cell just from these indexes. In the second case, we need an explicit list first to know how many cells there are. Even more, there is no way to link the nodes with the cells (back and forth) other than having a list of nodes and cells. Again, there are pros and cons for each of the grid types such as simplicity, flexibility, etc. In general, unstructured grids and better represent a certain geometry with the same number of cells. Structured grids suffer the so-called “staircase effect” that makes the unusable for discretizing mechanical parts.
Figure [-@fig:grids] illustrates how the same domain can be discretised using these two kind of grids. In the first case, we could identify any single cell by using just two indexes. We could even tell which nodes define each cell just from these indexes. In the second case, we need an explicit list first to know how many cells there are. Even more, there is no way to link the nodes with the cells (back and forth) other than having a list of nodes and cells. Again, there are pros and cons for each of the grid types such as simplicity, flexibility, etc. In general, unstructured grids and better represent a certain geometry with the same number of cells. Structured grids suffer the so-called “staircase effect” that makes the unusable for discretising mechanical parts.

::::: {#fig:grids}
![Continuous domain](dominio-continuo.svg){#fig:continuous width=30%}\
@@ -246,7 +245,7 @@ This section is not (just) about different kinds of elements like tetrahedra, he
- sub-integrated elements
- incomplete elements

And then there exist different pre-processors, meshers, solvers, pre-conditioners, post-processing steps, etc. A similar list can be made for the heat conduction problem, electromagnetism, the Schröedinger equation, neutron transport, etc. But there is also another level of “kind of problem,” which is related to how much accuracy and precision we are to willing sacrifice in order to have a (probably very much) simpler problem to solve. Again, there are different combinations here but a certain problem can be solved using any of the following three approaches, listed in increasing amount of difficulty and complexity:
And then there exist different pre-processors, meshers, solvers, pre-conditioners, post-processing steps, etc. A similar list can be made for the [heat conduction problem](https://en.wikipedia.org/wiki/Thermal_conduction), [electromagnetism](https://en.wikipedia.org/wiki/Electromagnetism), the [Schröedinger equation](https://en.wikipedia.org/wiki/Schr%C3%B6dinger_equation), [neutron transport](https://en.wikipedia.org/wiki/Neutron_transport), etc. But there is also another level of “kind of problem,” which is related to how much accuracy and precision we are to willing sacrifice in order to have a (probably very much) simpler problem to solve. Again, there are different combinations here but a certain problem can be solved using any of the following three approaches, listed in increasing amount of difficulty and complexity:

i. conservative
ii. best-estimate
@@ -264,7 +263,7 @@ Finally, when then uncertainties associated to the parameters, methods and model
2. performing a large number of runs for different combination of parameters, and
3. combining all the results into to obtain a best-estimate value plus uncertainty.

This kind of computation is usually required by the nuclear regulatory authorities when power plant designers need to address the safety of the reactors. What is the heat capacity of uranium above 1000ºC? What is the heat transfer coefficient when approaching the [critical heat flux](https://en.wikipedia.org/wiki/Critical_heat_flux) before the [Leidenfrost effect](https://en.wikipedia.org/wiki/Leidenfrost_effect) occurs? A certain statistical analysis has to be done prior to actually parametrically swifting the input parameters so as to obtain a distribution of possible outcomes.
This kind of computation is usually required by the nuclear regulatory authorities when power plant designers need to address the safety of the reactors. What is the heat capacity of uranium above 1000ºC? What is the heat transfer coefficient when approaching the [critical heat flux](https://en.wikipedia.org/wiki/Critical_heat_flux) before the [Leidenfrost effect](https://en.wikipedia.org/wiki/Leidenfrost_effect) occurs? A certain statistical analysis has to be done prior to actually parametrically sweeping (see\ [@sec:parametric]) the input parameters so as to obtain a distribution of possible outcomes.


## Five whys {#sec:five}
@@ -291,14 +290,14 @@ Here is an [original example](https://www.toyota-global.com/company/toyota_tradi
5. Why is the intake clogged with metal shavings?
Because there is no filter on the pump.

You get the point, even though we know thanks to Richard Feynmann that to answer a “why” question at some point we need to rely on the questioner’s previous experience. We usually assume we have to do what we usually do (i.e. perform finite element analysis). But do we? Do we add a filter or just replace the fuse?
You get the point, even though we know thanks to [Richard Feynmann](https://en.wikipedia.org/wiki/Richard_Feynman) that to answer a “why” question at some point we need to rely on the questioner’s previous experience. We usually assume we have to do what we usually do (i.e. perform finite element analysis). But do we? Do we add a filter or just replace the fuse?

Getting back to the case study: do we need to do FEM analysis? Well, it does not look like we can obtain the stresses the transient cases with just pencil and paper. But how much complexity should we add? We might do as little as axysimmetric linear steady-state conservative studies or as much as full three-dimensional non-linear transient best-estimate plus uncertainties computations. And here is where good engineers should appear: in putting their engineering judgment (call it experience or hunches) into defining what to solve. And it is not (just) because the first option is faster to solve than the latter. Involving many complex methods need more engineering time
Getting back to the case study: do we need to do FEM analysis? Well, it does not look like we can obtain the stresses the transient cases with just pencil and paper. But how much complexity should we add? We might do as little as axisymmetric linear steady-state conservative studies or as much as full three-dimensional non-linear transient best-estimate plus uncertainties computations. And here is where good engineers should appear: in putting their engineering judgement (call it experience or hunches) into defining what to solve. And it is not (just) because the first option is faster to solve than the latter. Involving many complex methods need more engineering time

1. to prepare the input data and set up the algorithms, and
2. to analyse the output data and write engineering reports.

In the first years of the history of computer, when programs where written in decks and output results were printed in continuous paper sheets, it made sense for computer programs to calculate and write as much data as possible even if it was not needed. One would never know if it would not be needed in the future, and CPU time was so expensive that re-running engineering computations because a results was not included in the output was forbidden. But that is not remotely true in the XXI century anymore. Computing time is far cheaper than engineering time (result known as the [UNIX Rule of Economy](http://www.catb.org/~esr/writings/taoup/html/ch01s06.html)) that it should be neglected with respect to the time spent by a cognizant engineer searching and sorting thousands of hard-to-read floating-point numbers.
In the first years of the history of computer, when programs where written in decks and output results were printed in continuous paper sheets, it made sense for computer programs to calculate and write as much data as possible even if it was not needed. One would never know if it would not be needed in the future, and CPU time was so expensive that re-running engineering computations because a results was not included in the output was forbidden. But that is not remotely true in the XXI century anymore. Computing time is far cheaper than engineering time (result known as the [UNIX Rule of Economy](http://www.catb.org/~esr/writings/taoup/html/ch01s06.html)) that it should be neglected with respect to the time spent by a cognisant engineer searching and sorting thousands of hard-to-read floating-point numbers.

divert(-1)

@@ -351,7 +350,7 @@ So we need to address the issue of fatigue in nuclear reactor pipes that
b. heat transients, and
c. seismic loads.

As I wanted to illustrate in [@sec:five], it is very important to decide what kind of problem (actually problems) we should be dealing with. As a nuclear engineer, I learned (theoretically in college but practically after college) that there are some models that let you see some effects and some that let you see other effects (please [say “modeling” not “simulation.”](https://www.seamplex.com/blog/say-modeling-not-simulation.html)]). And even if, in principle, it is true that more complex models should let you see more stuff, they definitely might show you nothing at all if the model is so big and complex that it does not fit into a computer (say because it needs hundreds of gigabytes of RAM to run) or because it takes more time to compute than you may have before the final report is expected.
As I wanted to illustrate in [@sec:five], it is very important to decide what kind of problem (actually problems) we should be dealing with. As a nuclear engineer, I learned (theoretically in college but practically after college) that there are some models that let you see some effects and some that let you see other effects (please [say “modelling” not “simulation.”](https://www.seamplex.com/blog/say-modeling-not-simulation.html)]). And even if, in principle, it is true that more complex models should let you see more stuff, they definitely might show you nothing at all if the model is so big and complex that it does not fit into a computer (say because it needs hundreds of gigabytes of RAM to run) or because it takes more time to compute than you may have before the final report is expected.

First of all, we should note that we need to solve

@@ -416,27 +415,27 @@ An example case where the SCLs are located around the junction between stainless



As an example, let us consider the system depicted in\ [@fig:valve-cad1;@fig:valve-scls1] where there is a stainless-carbon steel interface at the discharge of the valves. Instead of solving the transient heat-conduction problem with the internal temperature of the pipes equal to the temperature of the water in the reference transient condition of the power plant and an external condition of natural convection to the ambient temperature in the whole mesh of\ [@fig:valve-mesh1], a reduced model consisting of half of one of the two valves and a small length of the pipes at both the valve inlet and outlet is used. Once the temperature distribution\ $\hat{T}(\vec{x},t)$ for each time is obtained in the reduced mesh ([@fig:valve-temp], which has the origin at the center of the valve), the actual temperature distribution\ $T(\vec{x},t)$ is computed by an algebraic genearalisation of $\hat{T}(\vec{x},t)$ in the full coordinate system (where the origin is shown in\ [@fig:valve-cad1]). As stated above, those locations which are not covered by the reduced model are generalised with a time-dependent uniform temperature which is the average of the inner and outer temperatures at the inlet and outlet of the reduced mesh. The result is illustrated in figure\ [@fig:valve-gen].
As an example, let us consider the system depicted in\ [@fig:valve-cad1;@fig:valve-scls1] where there is a stainless-carbon steel interface at the discharge of the valves. Instead of solving the transient heat-conduction problem with the internal temperature of the pipes equal to the temperature of the water in the reference transient condition of the power plant and an external condition of natural convection to the ambient temperature in the whole mesh of\ [@fig:valve-mesh1], a reduced model consisting of half of one of the two valves and a small length of the pipes at both the valve inlet and outlet is used. Once the temperature distribution\ $\hat{T}(\vec{x},t)$ for each time is obtained in the reduced mesh ([@fig:valve-temp], which has the origin at the centre of the valve), the actual temperature distribution\ $T(\vec{x},t)$ is computed by an algebraic generalisation of $\hat{T}(\vec{x},t)$ in the full coordinate system (where the origin is shown in\ [@fig:valve-cad1]). As stated above, those locations which are not covered by the reduced model are generalised with a time-dependent uniform temperature which is the average of the inner and outer temperatures at the inlet and outlet of the reduced mesh. The result is illustrated in figure\ [@fig:valve-gen].

::::: {#fig:valve-temp-gen}
![Reduced mesh around the valve refined around the interface where the transient heat conduction problem is solved.](valve-temp.png){#fig:valve-temp width=80%}

![Generalization of the resulting temperature to the original mesh from figure [@fig:valve-mesh1]](valve-gen.png){#fig:valve-gen}
![Generalisation of the resulting temperature to the original mesh from figure [@fig:valve-mesh1]](valve-gen.png){#fig:valve-gen}

Computation of the thermal problem in a reduced mesh and generalisation of the result to the full original 3D mesh of figure\ [@fig:valve-mesh1].
:::::


Please note that there is no need to have a one-to-one correspondence between the elements from the reduced mesh with the elements from the original one. Actually, the reduced mesh contains first-order elements whilst the former has quadratic tetrahedra. Also the grid density is different. Nevertheless, the finite-element solver [Fino](https://www.seamplex.com/fino) used both to solve the heat and the mechanical problems, allows to read functions of space and time defined over one mesh and continuously evaluate and use them into another one even if the two grids have different elements, orders or even dimensions. In effect, in the system from figure [@fig:real-life] the material interface is between a orifice plate made in stainless steel that is welded to a carbon-steel pipe ([@fig:real-mesh2]). The thermal problem can be modeled using a two-dimensional axisymmetric grid [@fig:real-temp] and then generalized to the full three-dimensional mesh using the algebraic manipulation capabilities provided by [Fino](https://www.seamplex.com/fino) (actually by [wasora](https://www.seamplex.com/wasora)) as shown in\ [@fig:real-gen].
Please note that there is no need to have a one-to-one correspondence between the elements from the reduced mesh with the elements from the original one. Actually, the reduced mesh contains first-order elements whilst the former has quadratic tetrahedra. Also the grid density is different. Nevertheless, the finite-element solver [Fino](https://www.seamplex.com/fino) used both to solve the heat and the mechanical problems, allows to read functions of space and time defined over one mesh and continuously evaluate and use them into another one even if the two grids have different elements, orders or even dimensions. In effect, in the system from figure [@fig:real-life] the material interface is between a orifice plate made in stainless steel that is welded to a carbon-steel pipe ([@fig:real-mesh2]). The thermal problem can be modelled using a two-dimensional axi-symmetric grid [@fig:real-temp] and then generalised to the full three-dimensional mesh using the algebraic manipulation capabilities provided by [Fino](https://www.seamplex.com/fino) (actually by [wasora](https://www.seamplex.com/wasora)) as shown in\ [@fig:real-gen].

dnl 33410 10-12D-24

::::: {#fig:real}
![Original full 3D mesh, carbon steel is magenta, stainless steel is green.](real-mesh2.png){#fig:real-mesh2 width=75%}

![Temperature distribution for a certain time within the transient computed on a reduced two-dimensional axisymmetric mesh modeling half the orifice plate and a length of the carbon pipe.](real-temp.png){#fig:real-temp width=100%}
![Temperature distribution for a certain time within the transient computed on a reduced two-dimensional axi-symmetric mesh modelling half the orifice plate and a length of the carbon pipe.](real-temp.png){#fig:real-temp width=100%}

![Generalization of the temperature from the reduced mesh to the full three-dimensional grid.](real-gen.png){#fig:real-gen width=80%}
![Generalisation of the temperature from the reduced mesh to the full three-dimensional grid.](real-gen.png){#fig:real-gen width=80%}

The material interface in the system from [@fig:real-life] is configured by an orifice plate made of stainless steel welded to a carbon-steel pipe.
:::::
@@ -448,7 +447,7 @@ Before considering the actual mechanical problem that will give us the stress te

### Earthquake spectra

In case you are wondering, the answer is yes: all nuclear power plants are designed to withstand earthquakes. Of course, not all plants need the same level of reinforcements. Those built in large quiet plains will be, seismically speaking, cheaper than those located in more geologically active zones. Keep in mind that all the 54 Japanese nuclear power plants did structurally resist the 2011 earthquake, and all of the reactores were safely shut down. What actually happened in [Fukushima](http://www.world-nuclear.org/information-library/safety-and-security/safety-of-plants/fukushima-accident.aspx) is that one hour after the main shake, a 14-metre tsunami splashed on the coast, jumping over the 9-metre defenses and flooding the emergency Diesel generators that provided power to the pumps in charge of removing the remaining [decay power](https://en.wikipedia.org/wiki/Decay_heat) from the already-stopped reactor core.
In case you are wondering, the answer is yes: all nuclear power plants are designed to withstand earthquakes. Of course, not all plants need the same level of reinforcements. Those built in large quiet plains will be, seismically speaking, cheaper than those located in more geologically active zones. Keep in mind that all the 54 Japanese nuclear power plants did structurally resist the 2011 earthquake, and all of the reactors were safely shut down. What actually happened in [Fukushima](http://www.world-nuclear.org/information-library/safety-and-security/safety-of-plants/fukushima-accident.aspx) is that one hour after the main shake, a 14-metre tsunami splashed on the coast, jumping over the 9-metre defences and flooding the emergency Diesel generators that provided power to the pumps in charge of removing the remaining [decay power](https://en.wikipedia.org/wiki/Decay_heat) from the already-stopped reactor core.

Back to our case study, the point is that each site where nuclear power plants are built must have a geological study where a postulated design-basis earthquake is to be defined. In other words, a theoretical earthquake which the plant ought to withstand needs to be specified. How? By giving a set of three spectra (one for each coordinate direction) giving acceleration as a function of the frequency for each level of the building. That is to say, once the earthquake hits the power plant, depending on soil-structure interactions the energy will shake the building foundations in a way that depends on the characteristics of the earthquake, the soil and the concrete structure. Afterwards, they way the oscillations travel upward and shake each of the mechanical components erected in each floor level depends on the design of the civil structure in a way which is fully determined by the floor response spectra like the ones depicted in\ [@fig:spectrum].

@@ -468,7 +467,7 @@ Practically, these problems are solved using the same mechanical finite-element
1. The computation of the natural frequencies is “load free”, i.e. there can be no surface nor volumetric loads, and
2. The displacement boundary conditions ought to be homogeneous, i.e. only displacements equal to zero can be given. One may fix only one of the three degrees of freedom in certain surfaces and leave the others free though, as long as all the rigid body motions are removed as usual.

A real continuous solid has infinite modes of oscillation. A discretized one (using the most common and efficient FEM formulation, the [displacement-based formulation](http://web.mit.edu/kjb/www/Books/FEP_2nd_Edition_4th_Printing.pdf)) has three times the number of nodes modes of oscillation. In any case, one is usually interested only in a few of them, namely those with the lower frequencies because they take most of the energy with them. Each mode has two associated parameters called modal mass and excitation parameter that reflect how “important” the mode is regarding the absorption of energy from an external oscillatory source. Usually a couple of dozens of modes are enough to take up more than 90% of the earthquake energy. Figure\ [-@fig:modes] shows the first six natural modes of a sample piping section.
A real continuous solid has infinite modes of oscillation. A discretised one (using the most common and efficient FEM formulation, the [displacement-based formulation](http://web.mit.edu/kjb/www/Books/FEP_2nd_Edition_4th_Printing.pdf)) has three times the number of nodes modes of oscillation. In any case, one is usually interested only in a few of them, namely those with the lower frequencies because they take most of the energy with them. Each mode has two associated parameters called modal mass and excitation parameter that reflect how “important” the mode is regarding the absorption of energy from an external oscillatory source. Usually a couple of dozens of modes are enough to take up more than 90% of the earthquake energy. Figure\ [-@fig:modes] shows the first six natural modes of a sample piping section.

::::: {#fig:modes}
![$i=1$](mode1.png){width=48%}\
@@ -483,7 +482,7 @@ A real continuous solid has infinite modes of oscillation. A discretized one (us
First six natural oscillation modes for a piping section.
:::::

These first modes that take up most of the energy are then used to take into account the earthquake load. There are several ways of performing this computation, but the ASME\ III code states that the method known as SRSS (for Square Root of Sum of Squares) can be safely used. This method mixes the eigenvectors with the floor response spectra through the eigenvalues and gives an spatial (actually nodal) distribution of three accelerations (one for each direction) that, when multiplied by the density give a vector of a distributed force (in units of Newton per cubic millimeter for example) which is statically equivalent to the load coming from the postulated earthquake.
These first modes that take up most of the energy are then used to take into account the earthquake load. There are several ways of performing this computation, but the ASME\ III code states that the method known as SRSS (for Square Root of Sum of Squares) can be safely used. This method mixes the eigenvectors with the floor response spectra through the eigenvalues and gives an spatial (actually nodal) distribution of three accelerations (one for each direction) that, when multiplied by the density give a vector of a distributed force (in units of Newton per cubic millimetre for example) which is statically equivalent to the load coming from the postulated earthquake.

::::: {#fig:acceleration}
![$a_x$](ax.png){width=80%}
@@ -498,7 +497,7 @@ The equivalent accelerations for the piping section of [@fig:modes] for the spec
The ASME code says that these accelerations (depicted in [@fig:acceleration]) are to be applied twice. Once with the original sign and once with all the elements with the opposite sign during two seconds of the transient each time.


## Linearity (not yet linearization)
## Linearity (not yet linearisation)

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.

@@ -615,7 +614,7 @@ Spatial distribution of principal stress\ 3 for cases\ B and\ C. If linearity ap
:::::

In the second case, the principal stresses are not uniform and have a non-trivial distribution. Indeed, the distribution of\ $\sigma_3$ obtained by CAeplex is shown in\ [@fig:cube-shear].
Now, if we indeed were facing a fully linear problem then the results of the sum of two inputs would be equal to the sum of the individual inputs. And\ [@fig:cube-full], which shows the principal stress\ 3 of case\ C is not the result from case\ B plus any of the three constants from case\ A. Had it been, the color distribution would have been _exactly_ the same as the scale goes automatically from the most negative value in blue to the most positive value in red. And 7+30\ $\neq$ 33. Alas, it seems that there exists some kind of unexpected non-linearity (the feared master-pet interaction?) that prevents us from from fully splitting the problem into simpler chunks.
Now, if we indeed were facing a fully linear problem then the results of the sum of two inputs would be equal to the sum of the individual inputs. And\ [@fig:cube-full], which shows the principal stress\ 3 of case\ C is not the result from case\ B plus any of the three constants from case\ A. Had it been, the colour distribution would have been _exactly_ the same as the scale goes automatically from the most negative value in blue to the most positive value in red. And 7+30\ $\neq$ 33. Alas, it seems that there exists some kind of unexpected non-linearity (the feared master-pet interaction?) that prevents us from from fully splitting the problem into simpler chunks.

So what is the source of this unexpected non-linear effect in an otherwise nice and friendly linear formulation? Well, probably you already know it because after all it is almost high-school mathematics. But I learned it way after college when facing a real engineering problem and not just back-of-the-envelope pencil-and-paper trivial exercises.

@@ -661,14 +660,14 @@ ans =
5.767255
```

Did I convince you? More or less, right? The third eigenvalue seems to fit. Let us not throw all of our beloved linearity away and dig in further into the subject. There are still two important issues to discuss which can be easily addressed using fresh-year linear algebra (remember, do not fear math!). First of all, even though principal stresses are not linear with respect to the sum they are linear with respect to pure multiplication. Once more, think what happens to the the eigenvalues and eigenvectors of a single stress tensor as all its elements are scaled up or down by a real scalar. They are the same! So, for example, the [Von Mises stress](https://en.wikipedia.org/wiki/Von_Mises_yield_criterion) (which is a combination of the principal stresses) of a beam loaded with a force\ $\alpha \cdot \vec{F}$ is\ $\alpha$ times the stress of the beam loaded with a force\ $\vec{F}$. Please test this hypothesis by playing with your favorite FEM solver an play. Or even better, take a look at the stress invariants $I_1$, $I_2$ and $I_3$ (you can search online or peek into the source code of [Fino](https://www.seamplex.com/fino) and grep for the routine called `fino_compute_principal_stress()`) and see (using paper and pencil!) how they scale up if the individual elements of the stress tensor are scaled by a real factor\ $\alpha$.
Did I convince you? More or less, right? The third eigenvalue seems to fit. Let us not throw all of our beloved linearity away and dig in further into the subject. There are still two important issues to discuss which can be easily addressed using fresh-year linear algebra (remember, do not fear maths!). First of all, even though principal stresses are not linear with respect to the sum they are linear with respect to pure multiplication. Once more, think what happens to the the eigenvalues and eigenvectors of a single stress tensor as all its elements are scaled up or down by a real scalar. They are the same! So, for example, the [Von Mises stress](https://en.wikipedia.org/wiki/Von_Mises_yield_criterion) (which is a combination of the principal stresses) of a beam loaded with a force\ $\alpha \cdot \vec{F}$ is\ $\alpha$ times the stress of the beam loaded with a force\ $\vec{F}$. Please test this hypothesis by playing with your favorite FEM solver an play. Or even better, take a look at the stress invariants $I_1$, $I_2$ and $I_3$ (you can search online or peek into the source code of [Fino](https://www.seamplex.com/fino) and grep for the routine called `fino_compute_principal_stress()`) and see (using paper and pencil!) how they scale up if the individual elements of the stress tensor are scaled by a real factor\ $\alpha$.

The other issue is that even though in general the eigenvalues of the sum of two matrices are not the same as the eigenvalues of the matrix sum, there are some cases when they are. In effect, if two matrices\ $A$ and\ $B$ commute, i.e. their product is commutative

$$
A \cdot B = B \cdot A
$$
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.
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 diagonalisable using the same basis. That is to say, the diagonalising matrix\ $P$ such that $P^{-1} A P$ is diagonal is the same that renders\ $P^{-1} B P$ also diagonal. What does this mechanically 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.


## ASME stress linearisation (not linearity!)
@@ -699,11 +698,11 @@ No need to know or even understand these two integrals which for sure are not in

# 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.
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 pressurised pipe once again. It is time to solve the problem with a computer using finite elements and obtain some funny coloured pictures instead of just equations.

The first thing that has to be said is that, as with any interesting problem, there are literally hundreds or different ways of solving it, each of them throwing particular conclusions. For example, one can:

1. solve a real 3D problem or a 2D axisymmetric case (or even a 1D case using the collocation method to solve the radial dependence),
1. solve a real 3D problem or a 2D axi-symmetric case (or even a 1D case using the collocation method to solve the radial dependence),
2. have a full cylindrical geometry or just a half (180º), or a quarter (90º), or a thin slice (a small amount of degrees),
3. use a structured or an unstructured grid,
4. uniformly or locally-refine the mesh (with several choices of refinement),
@@ -718,22 +717,22 @@ The first thing that has to be said is that, as with any interesting problem, th
![Structured second-order incomplete hexahedra](quarter-struct.png){#fig:cube-struct width=48%}\
![Unstructured second-order tetrahedra](quarter-caeplex.png){#fig:quarter-caeplex width=48%}

Two of the hundreds of different ways the infinite pressurized pipe can be solved using FEM
Two of the hundreds of different ways the infinite pressurised pipe can be solved using FEM
:::::


You can get both the exponential nature of each added bullet and how easily we can add new further choices to solve a FEM problem. And each of these choices will reveal you something about the nature of either the mechanical problem or the numerical solution. It is not possible to teach any possible lesson from every outcome in college, so you will have to learn them by yourself getting hands at them. I have already tried to address the particular case of the infinite pipe in a [recent report](https://www.seamplex.com/fino/doc/pipe-linearized/) that is worth reading before carrying on with this article. The conclusions of the report are:

* Engineering problems ought not to be solved using black-boxes (i.e. privative software whose source code is not freely available)---more on the subject below in\ [@sec:two-materials].
* The pressurized infinite pipe has only one independent variable (the radius $r$) and one primary dependent variable (the radial displacement $u_r$).
* The pressurised infinite pipe has only one independent variable (the radius $r$) and one primary dependent variable (the radial displacement $u_r$).
* The problem has analytical solution for the radial displacement and the radial, tangential and axial stresses.
* There are no shear stresses, so these three stresses are also the principal stresses.
* Analytical expressions for the membrane and membrane plus bending stresses along any radial SCL can be obtained.
* The spatial domain can be discretized using linear or higher-order elements. In particular first and second-order elements have been used in this report.
* The spatial domain can be discretised using linear or higher-order elements. In particular first and second-order elements have been used in this report.
* For the same number of elements, second-order grids need more nodes than linear ones, although they can better represent curved geometries.
* The discretized problem size depends on the number of nodes and not on the number of elements.
* The discretised problem size depends on the number of nodes and not on the number of elements.
* The finite-element results for the displacements are similar to the analytical solution, with second-order grids giving better results for the same number of elements (this is expected as they involved far more nodes).
* The three stress distributions computed with the finite-element give far more reasonable results for the second-order case than for the first-order grid. This is qualitatively explained by the fact that first-order tetrahedra have uniform derivatives and such the elements located in both the external and external faces represent the stresses not at the actual faces but at the barycenter of the elements.
* The three stress distributions computed with the finite-element give far more reasonable results for the second-order case than for the first-order grid. This is qualitatively explained by the fact that first-order tetrahedra have uniform derivatives and such the elements located in both the external and external faces represent the stresses not at the actual faces but at the barycentre of the elements.
* Membrane stresses converge well for both the first and second-order cases because they represent a zeroth-order moment of the stress distribution and such the excess and defect errors committed by the internal and external elements approximately cancel out.
* Membrane plus bending stresses converge very poorly in first-order grids because the excess and defect errors do not cancel out because it is a first-order moment of the stress distribution.
* The computational effort to solve a given problem, namely the CPU time and the needed RAM depend non-linearly on various factors, but the most important one is the problem size which is three times the number of nodes in the grid.
@@ -750,11 +749,11 @@ Error in the computation of the linearised stresses vs. CPU time needed to solve

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 with the [displacement-based formulation](http://web.mit.edu/kjb/www/Books/FEP_2nd_Edition_4th_Printing.pdf) formulation) the number of _nodes_.

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.
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 definitely be the number of nodes rather than the number of elements. Let us see why.

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?

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.
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 abscissas rather than the number of elements to compare first and second-order formulations.

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:

@@ -775,13 +774,13 @@ For a fixed number of nodes, first-order grids have far more elements than secon

### 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):
The linear FEM problem leads of course of a system of\ $NG$ linear equations, cast 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).
The objective of the solver is to find the vector\ $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 \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.
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 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%}\
@@ -803,7 +802,7 @@ Let us focus on the first item and leave the second one for a separate discussio

**FIGURA 1D**

Math shows that the location where the derivatives of the interpolated displacements are closer to the real (i.e. the analytical in problem that have it) solution are the elements’ [Gauss points](https://en.wikipedia.org/wiki/Gaussian_quadrature). Even better, the material properties at these points are continuous (they are usually uniform but they can depend on temperature for example) because, unless we are using weird elements, there are no material interfaces inside elements. But how to manage a set of stresses given at the Gauss points instead of at the nodes? Should we use one mesh for the input and another one for the output? What happens when we need to know the stresses on a surface and not just in the bulk of the solid? There are still no one-size-fits-all answers. There is a very interesting [blog post](http://tor-eng.com/2017/11/practical-tips-dealing-stress-singularities/) by Nick Stevens that addresses the issue of stresses computed at sharp corners. What does your favourite FEM program do with such a case?
Detailed mathematics show that the location where the derivatives of the interpolated displacements are closer to the real (i.e. the analytical in problem that have it) solution are the elements’ [Gauss points](https://en.wikipedia.org/wiki/Gaussian_quadrature). Even better, the material properties at these points are continuous (they are usually uniform but they can depend on temperature for example) because, unless we are using weird elements, there are no material interfaces inside elements. But how to manage a set of stresses given at the Gauss points instead of at the nodes? Should we use one mesh for the input and another one for the output? What happens when we need to know the stresses on a surface and not just in the bulk of the solid? There are still no one-size-fits-all answers. There is a very interesting [blog post](http://tor-eng.com/2017/11/practical-tips-dealing-stress-singularities/) by Nick Stevens that addresses the issue of stresses computed at sharp corners. What does your favourite FEM program do with such a case?

In any case, this step takes a non-negligible amount of time. The most-common approach, i.e. the node-averaging method is driven mainly by the number of nodes of course. So all-in-all, these are the reasons to use the number of nodes instead of the numbers of elements as a basic parameter to measure the complexity of a FEM problem.

@@ -865,9 +864,33 @@ In effect, a couple of years ago Angus Ramsay noted [a weird behaviour](https://

> Do you trust your favourite FEM program?

Back to the two-material interface, all the discussion above about non-continuous derivatives applies to a sharp abrupt interface. In the study case the junctions are welded so there is a heat-affected zone with changes in the material microstructure. Therefore, there exist a smooth transition from the mechanical properties of one material to the other one in a way that is very hard to predict and to model. In principle, the assumption of a sharp interface is conservative (in the sense of\ [@sec:kinds]). There cannot be an SCL exactly on a material interface so there should be at least two SCLs, one at each side of the junction as\ [@fig:weldolet-scls;@fig:valve-scls1] illustrate. The actual distance would have to be determined first as an educated guess, then with trial and error and finally in accordance with the regulator.
Back to the two-material interface, all the discussion above about non-continuous derivatives applies to a sharp abrupt interface. In the study case the junctions are welded so there is a heat-affected zone with changes in the material micro structure. Therefore, there exist a smooth transition from the mechanical properties of one material to the other one in a way that is very hard to predict and to model. In principle, the assumption of a sharp interface is conservative (in the sense of\ [@sec:kinds]). There cannot be an SCL exactly on a material interface so there should be at least two SCLs, one at each side of the junction as\ [@fig:weldolet-scls;@fig:valve-scls1] illustrate. The actual distance would have to be determined first as an educated guess, then with trial and error and finally in accordance with the regulator.

## A parametric tee {#sec:parametric}

Time for another experiment. We know (more or less) what to expect from an infinite pressurised pipe. What if we added a branch to such pipe? Even more, what if we added successively (one at a time) branches with different diameters? This is called parametric analysis, and sooner or later (if not now) you will find yourself performing this kind of computations more often than any other one.

So here comes the five Feynmann-Ohno questions:

1. Why do you want to perform a parametric computation?
So we can study how the linearised stresses change with the inclusion of a branch of a certain diameter.

2. Why do you want to know how the stresses change with the inclusion of a branch?
Because it is like training throw-ins with watermelons during the week so you can reach the area with a regular ball on the weekends (football joke).

3. Why do you want to do something harder?
Because [we can](https://www.youtube.com/watch?v=BVd-rYIqSy8), and getting out of our comfort zone once in a while is a healthy habit.

4. Why is getting out of our comfort zone a healthy habit?
Because it fires up a certain part of our brains that keep us active and prevents deseases so we can live longer.

5. Why do you want to live longer?
To perform more parametric finite-element analysis.


Let us solve the following mock-up case: a long main 12-inch schedule\ 80 pipe has an orthogonal branch of a certain nominal diameter of\ $d_b$\ inches (it seems that the [SI](https://en.wikipedia.org/wiki/International_System_of_Units) did not do well amongst the piping engineering community). Both the main line and the branch are pressurised with\ $p=10$\ MPa. The main pipe is aligned with the\ $x$ axis and the branch coincides with the\ $z$ axis. Thus, the $x$-$z$ plane acts as a symmetry plane and we only need to model one octant of the full geometry.


## A parametric tee

## Bake, shake and break {#sec:break}


Escribir Vista previa
Cargando…
Cancelar
Guardar