ソースを参照

grammar fixes

tags/CANDIDATE
gtheler 7年前
コミット
6d43aebaed
2個のファイルの変更108行の追加89行の削除
  1. +6
    -0
      make.sh
  2. +102
    -89
      nafems4.md

+ 6
- 0
make.sh ファイルの表示

@@ -140,4 +140,10 @@ for format in ${formats}; do
if [ -z "`git check-ignore ${name}.${format}`" ]; then
echo ${name}.${format} >> .gitignore
fi
# post-process
if [[ "${format}" == "html" ]]; then
sed -i 's/<table>/<table class="table table-responsive table-striped">/' ${name}.${format}
sed -i 's/<blockquote>/<blockquote class="blockquote">/' ${name}.${format}
fi
done

+ 102
- 89
nafems4.md ファイルの表示

@@ -2,9 +2,11 @@ dnl $$\renewcommand\mathbf[1]{\mathbf{#1}}$$

# Background and introduction

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.
First of all, please take this text as a written chat between you an me, i.e. an average engineer that has 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 “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.
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 while students are 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.

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 a 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](https://en.wikipedia.org/wiki/Galileo_Galilei) 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 have later studied 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 analysed 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. Yes, this is my personal story but it could have easily been yours as well. My point is that the very same distance between what I imagined 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 ([@fig:pipes]). 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.

::::: {#fig:pendulum}
![Simple pendulum](simple.svg){#fig:simple width=35%}\
@@ -20,14 +22,12 @@ A simple pendulum from college physics courses and a real-life pendulum. Hint: t
An infinitely-long pressurised thick pipe as taught in college and an isometric drawing of a section of a real-life piping system.
:::::

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. Yes, 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.
Like the pendulums above, we will be swinging back and forth between a case study about fatigue analysis in piping systems of a nuclear power plant and more generic and even romantic topics related to finite elements and computational mechanics. These latter regressions will not remain just as abstract theoretical ideas. Not only will they be directly applicable to the development of the main case, but they will also apply to a great deal of other engineering problems tackled with the finite element method (and its cousins, [@sec:formulations]).

Finite elements are like magic to me. I mean, I can follow the whole derivation of the equations, from the strong, weak and variational formulations of the equilibrium equations for the mechanical problem (or the energy conservation for heat transfer) down to the [algebraic multigrid](https://en.wikipedia.org/wiki/Multigrid_method) [preconditioner](https://en.wikipedia.org/wiki/Preconditioner) for the inversion of the [stiffness matrix](https://en.wikipedia.org/wiki/Stiffness_matrix) passing through [Sobolev spaces](https://en.wikipedia.org/wiki/Sobolev_space) and the [grid generation](https://en.wikipedia.org/wiki/Mesh_generation). Then I can sit down and program all these steps into a computer, including the shape functions and its derivatives, the assembly of the discretised stiffness matrix 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.
Finite elements are like magic to me. I mean, I can follow the whole derivation of the equations, from the strong, weak and variational formulations of the equilibrium equations for the mechanical problem (or the energy conservation for heat transfer) down to the [algebraic multigrid](https://en.wikipedia.org/wiki/Multigrid_method) [preconditioner](https://en.wikipedia.org/wiki/Preconditioner) for the inversion of the [stiffness matrix](https://en.wikipedia.org/wiki/Stiffness_matrix) passing through [Sobolev spaces](https://en.wikipedia.org/wiki/Sobolev_space) and the [grid generation](https://en.wikipedia.org/wiki/Mesh_generation). Then I can sit down and program all these steps into a computer, including the shape functions and its derivatives, the assembly of the discretised stiffness matrix ([@sec:building]), the numerical solution of the system of equations ([@sec:solving]) and the computation of the gradient of the solution ([@sec:stress-computation). Yet, the fact that all these a-priori unconnected steps give rise to pretty pictures that resemble reality is still astonishing to me.


Again, take all this information as coming from a fellow that has already taken such a journey from college’s pencil and paper to real engineering cases involving complex numerical calculations. And developing, in the meantime, both an actual working finite-element [back-end](https://www.seamplex.com/fino) and [front-end](https://www.caeplex.com) from scratch (the dean of the engineering school I attended used to say “It is not the same to read than to write manuals, and we should aim at writing.”).
Again, take all this information as coming from a fellow that has already taken such a journey from college’s pencil and paper to real engineering cases involving complex numerical calculations. And developing, in the meantime, both an actual working finite-element [back-end](https://www.seamplex.com/fino) and [front-end](https://www.caeplex.com) from scratch.^[The dean of the engineering school I attended used to say “It is not the same to read than to write manuals, and we should aim at writing them.”]


## Tips and tricks
@@ -35,30 +35,30 @@ Again, take all this information as coming from a fellow that has already taken
There are some useful tricks that come handy when trying to solve a mechanical problem. Throughout this text, I will try to tell you some of them.


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.
One of the most important ones to use 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 thermally-induced stresses are. 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 new concept. Nevertheless, there will be particular locations of the text where imagination will be most useful. I will bring the subject up.

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.
Another heads up is that we will dig into some mathematics. Probably they would be be simple and you would deal with them 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 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.
So here comes another experience tip: do not fear equations. Even more, keep exercising. You have used [differences of squares](https://en.wikipedia.org/wiki/Difference_of_two_squares) in high school, didn’t you?. You know (or at least knew) how to [integrate by parts](https://en.wikipedia.org/wiki/Integration_by_parts). Do you remember what [Laplace transforms](https://en.wikipedia.org/wiki/Laplace_transform) 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](https://en.wikipedia.org/wiki/Quotient_rule). 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.

Throughout the text I will be referring to “your favourite FEM program.” I bet you do have one. Mine is [CAEplex](https://caeplex.com) (it works on top of [Fino](https://www.seamplex.com/fino)). We will be using it to perform some tests and play a little bit. And we will also 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.
One final comment: throughout the text I will be referring to “your favourite FEM program.” I bet you do have one. Mine is [CAEplex](https://caeplex.com) (it works on top of [Fino](https://www.seamplex.com/fino)). We will be using it to perform some tests and play a little bit. And we will also 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.


# Case study: nuclear reactors, pressurised pipes and fatigue {#sec:case}

Piping systems in sensitive industries like nuclear or oil & gas should be designed and analysed following the recommendations of an appropriate set of codes and norms, such as the ASME\ Boiler and Pressure Vessel Code.
This code of practice (book) was born during the late\ XIX century, before finite-element methods for solving partial differential equations were even developed. And much longer before they were available for the general engineering community. Therefore, much of the code assumes design and verification is not necessarily performed numerically but with paper and pencil (yes, like in college). However, it still provides genuine guidance in order to ensure pressurised systems behave safely and properly without needing to resort to computational tools. Combining finite-element analysis with the ASME code gives the cognisant engineer a unique combination of tools to tackle the problem of designing and/or verifying pressurised piping systems.
Piping systems in sensitive industries like nuclear or oil & gas should be designed and analysed following the recommendations of an appropriate set of codes and norms, such as the [ASME\ Boiler and Pressure Vessel Code](https://en.wikipedia.org/wiki/ASME_Boiler_and_Pressure_Vessel_Code).
This code of practice was born during the late\ XIX century, before finite-element methods for solving partial differential equations were even developed. And much longer before they were available for the general engineering community. Therefore, much of the code assumes design and verification is not necessarily performed numerically but with paper and pencil (yes, like in college). However, it still provides genuine guidance in order to ensure pressurised systems behave safely and properly without needing to resort to computational tools. Combining finite-element analysis with the ASME code gives the cognisant engineer a unique combination of tools to tackle the problem of designing and/or verifying pressurised piping systems.

In the years following [Enrico Fermi](https://en.wikipedia.org/wiki/Enrico_Fermi)’s demonstration that a self-sustainable fission reaction chain was possible (actually, in fact, after WWII was over), people started to build plants in order to transform the energy stored within the atoms nuclei into usable electrical power. They quickly reached the conclusion that high-pressure heat exchangers and turbines were needed. So they started to follow the ASME\ Boiler and Pressure Vessel Code. They also realised that some requirements did not fit the needs of the nuclear industry. But instead of writing a new code from scratch, they added a new chapter to the existing body of knowledge: the celebrated ASME Section\ III.
In the years following [Enrico Fermi](https://en.wikipedia.org/wiki/Enrico_Fermi)’s demonstration that a self-sustainable [fission reaction](https://en.wikipedia.org/wiki/Nuclear_fission) chain was possible (actually, in fact, after WWII was over), people started to build plants in order to transform the energy stored within the atoms nuclei into usable electrical power. They quickly reached the conclusion that high-pressure heat exchangers and turbines were needed. So they started to follow codes of practise like the aforementioned ASME\ B&PVC. They also realised that some requirements did not fit the needs of the nuclear industry. But instead of writing a new code from scratch, they added a new chapter to the existing body of knowledge: the celebrated [ASME Section\ III](https://en.wikipedia.org/wiki/ASME_Boiler_and_Pressure_Vessel_Code#ASME_BPVC_Section_III_-_Rules_for_Construction_of_Nuclear_Facility_Components).

After further years passed by, engineers (probably the same people that forked section\ III) noticed that fatigue in nuclear power plants was not exactly the same as in other piping systems. There were some environmental factors directly associated to the power plant that was not taken into account by the regular ASME code. Again, instead of writing a new code from scratch, people decided to add correction factors to the previously amended body of knowledge. This is how knowledge evolves, and it is this kind of complexities that engineers are faced with during their professional lives. We have to face it, it would be a very hard work to re-write everything from scratch every time something changes.
After further years passed by, engineers (probably the same people that forked section\ III) noticed that fatigue in nuclear power plants was not exactly the same as in other piping systems. There were some environmental factors directly associated to the power plant that were not taken into account by the regular ASME code. Again, instead of writing a new code from scratch, people decided to add correction factors to the previously-amended body of knowledge. This is how (sometimes) knowledge evolves, and it is this kind of complexities that engineers are faced with during their professional lives. We have to admit it, it would be a very difficult to re-write everything from scratch every time something changes.

Actually, this article does not focus on a single case study but on some general ideas regarding analysis of fatigue in piping systems in nuclear power plants. There is no single case study but a compendium of ideas obtained by studying many different systems which are directly related to the safety of a real nuclear reactor.

![CAD model for the piping system of fig\ [@fig:isometric]](real-piping.png){#fig:real-life}
![CAD model for the piping system of\ [@fig:isometric]](real-piping.png){#fig:real-life}

## Nuclear reactors

@@ -67,31 +67,31 @@ In each of the countries that have at least one nuclear power plant there exists
## Pressurised pipes


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.
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 of the energy stored in the uranium nuclei into electricity we 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 temperatures and pressures near those values in many systems of the plant, especially in the primary circuit and 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 refuelled, they may undergo operational (and some incidental) 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% of their maximum power capacity. They need to be maintained and refuelled, they may undergo operational (and some incidental) transients, they might operate at a lower power due to load following conditions, etc. These transient cases involve 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 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.
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, both the civil structures and the piping system itself need to be able to withstand such a load, even if it occurs at the moment of highest mechanical demand during one of the operational transients.

## Fatigue {#sec:fatigue}

Mechanical systems can fail due to a wide variety of reasons. The effect known as fatigue can create, migrate and grow microscopic cracks at the atomic level, called dislocations. Once these cracks reach a critical size, then the material fails catastrophically even under stresses lower than tensile strength limits. There are not complete mechanistic models from first principles which can be used in general situations, and those that exist are very complex and hard to use. Instead, using an experimental approach very much like the Hooke Law experiment, the stress amplitude of a periodic cycle can be related to the number of cycles where failure by fatigue is expected to occur. For each material, this dependence can be computed using normalised tests and a family of “fatigue curves” like the one depicted in\ [@fig:SN] (also called $S$-$N$ curve) for different temperatures can be obtained.
Mechanical systems can fail due to a wide variety of reasons. The effect known as fatigue can create, migrate and grow microscopic cracks at the atomic level, called dislocations. Once these cracks reach a critical size, then the material fails catastrophically even under stresses much lower than the tensile strength limits. There are not complete mechanistic models from first principles which can be used in general situations, and those that exist are very complex and hard to use. Instead, using an experimental approach very much like the [Hooke’s Law](https://en.wikipedia.org/wiki/Hooke's_law) experiment, the stress amplitude of a periodic cycle can be related to the number of cycles where failure by fatigue is expected to occur. For each material, this dependence can be computed using normalised tests and a family of “fatigue curves” like the one depicted in\ [@fig:SN] (also called $S$-$N$ curve) for different temperatures can be obtained.


![A fatigue or $S$-$N$ curve for two steels.](SN.svg){#fig:SN width=80%}
![A fatigue or $S$-$N$ curve for two steels.](SN.svg){#fig:SN width=95%}

It should be stressed that the fatigue curves are obtained in a particular load case, namely purely-periodic one-dimensional, which is not directly generalised to other three-dimensional cases. Also, any real-life case will be subject to a mixture of complex cycles given by a stress time history and not to pure periodic conditions. The application of the curve data implies a set of simplifications and assumptions that are translated into different possible “rules” for composing real-life cycles. There also exist two safety factors which increase the stress amplitude and reduce the number of cycles respectively. All these intermediate steps render the analysis of fatigue into a conservative computation scheme. Therefore, when a fatigue analysis performed using the fatigue curve method arrives at the conclusion that “fatigue is expected to occur after ten thousand cycles” what it actually means is “we are sure fatigue will not occur before ten thousand cycles, yet it may not occur before one hundred thousand or even more.”
It should be noted that the fatigue curves are obtained in a particular load case, namely purely-periodic and one-dimensional, which cannot be directly generalised to other three-dimensional cases. Also, any real-life case will be subject to a mixture of complex cycles given by a stress time history and not to pure periodic conditions. The application of the curve data implies a set of simplifications and assumptions that are translated into different possible “rules” for composing real-life cycles. There also exist two safety factors which increase the stress amplitude and reduce the number of cycles respectively. All these intermediate steps render the analysis of fatigue into a conservative computation scheme ([@sec:kinds]). Therefore, when a fatigue analysis performed using the fatigue curve method arrives at the conclusion that “fatigue is expected to occur after ten thousand cycles” what it actually means is “we are sure fatigue will not occur before ten thousand cycles, yet it may not occur before one hundred thousand or even more.”


# Solid mechanics, or what we are taught at college

So, let us start our journey. Our starting place: undergraduate solid mechanics courses. Our goal: to obtain the internal state of a solid subject to a set of movement restrictions and loads (i.e. to solve the solid mechanics problem). Our first step: Newton’s laws of motion. For each of them, all we need to recall here is that
So, let us start our journey. Our starting place: undergraduate solid mechanics courses. Our goal: to obtain the internal state of a solid subject to a set of movement restrictions and loads (i.e. to solve the solid mechanics problem). Our first step: [Newton’s laws of motion](https://en.wikipedia.org/wiki/Newton's_laws_of_motion). For each of them, all we need to recall here is that

1. a solid is in equilibrium if it is not moving in at least one inertial coordinate system,
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](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.
We have to accept that there is certain intellectual beauty when complex stuff can be expressed in such simple terms. 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}

@@ -117,26 +117,26 @@ It looks (and works) like a regular $3 \times 3$ matrix. Some brief comments abo
- $\tau_{xy} = \tau_{yx}$,
- $\tau_{yz} = \tau_{zy}$, and
- $\tau_{zx} = \tau_{xz}$.
* The elements of the tensor depends on the orientation of the coordinate system.
* There exists a particular coordinate system in which the stress tensor is diagonal, i.e. all the shear stresses are zero. In this case, the three diagonal elements are called the principal stresses, which happen to be the three [eigenvalues](https://en.wikipedia.org/wiki/Eigenvalues_and_eigenvectors) of the stress tensor. The basis of the coordinate system that renders the tensor diagonal are the eigenvectors.
* The elements of the tensor depend on the orientation of the coordinate system.
* There exists a particular coordinate system in which the stress tensor is diagonal, i.e. all the shear stresses are zero. In this case, the three diagonal elements are called the [principal stresses](https://en.wikipedia.org/wiki/Principal_stresses), which happen to be the three [eigenvalues](https://en.wikipedia.org/wiki/Eigenvalues_and_eigenvectors) of the stress tensor. The basis of the coordinate system that renders the tensor diagonal are the eigenvectors.


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 $\mathbf{n}$, all we have to do is “project” the stress tensor through $\mathbf{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 your geometry, then... Congratulations! You have solved the solid mechanics problem.


## An infinitely-long pressurised pipe {#sec:infinite-pipe}

Let us proceed to our second step, and consider the infinite pipe subject to uniform internal pressure already introduced in\ [@fig:infinite-pipe]. 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$.
Let us proceed to our second step, and consider the infinite pipe subject to uniform internal pressure already introduced in\ [@fig:infinite-pipe]. 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. In fact, 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$.

### 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 pressurised pipe deforms to counteract the internal pressure\ $p$.
Remember that when any solid body is subject to external forces, it has to react 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.
* There are only radial displacements\ $u_r$ and they depend only on the radial coordinate\ $r$ and not on the axial position\ $z$ or on the azimuthal angle\ $\theta$. This displacements are
* There are only radial displacements\ $u_r$ and they depend only on the radial coordinate\ $r$ and not on the axial position\ $z$ or on the azimuthal angle\ $\theta$. These displacements are

$$
u_r(r) = p \cdot \frac{1+\nu}{E} \cdot \frac{a^2}{b^2-a^2} \cdot \left[ 1-2\nu + \frac{b^2}{r^2} \right]\cdot r
@@ -182,7 +182,7 @@ That is all what we can say about an infinite pipe with uniform material propert
* the material properties are not uniform (say the pipe does not have an uniform temperature but a distribution), or
* the pressure was not uniform (say because there is liquid inside and its weight cannot be neglected),

\noindent then we would no longer be able to fully solve the problem with paper and pencil and draw all the conclusions above. However, at least we have a start because we know that if the pipe is finite but long enough or the temperature is not uniform but almost, we still can use the analytical equations as approximations. After all, [Enrico Fermi](https://en.wikipedia.org/wiki/Enrico_Fermi) managed to reach criticality in the [Chicago Pile-1](https://en.wikipedia.org/wiki/Chicago_Pile-1) with paper and pencil. But what happens if the pipe is short, there are branches and temperature changes like during a transient in a nuclear reactor? Well, that is why we have finite elements. And this is were what we learned at college ends.
\noindent then we would no longer be able to fully solve the problem with paper and pencil and draw all the conclusions above. However, at least we have a start because we know that if the pipe is finite but long enough or the temperature is not uniform but almost, we still can use the analytical equations as approximations. After all, [Enrico Fermi](https://en.wikipedia.org/wiki/Enrico_Fermi) managed to reach criticality in the [Chicago Pile-1](https://en.wikipedia.org/wiki/Chicago_Pile-1) with paper and pencil. But what happens if the pipe is short, there are branches and temperature changes like during a transient in a nuclear reactor? Well, that is why we have finite elements. And this is where what we learned at college pretty much ends.



@@ -192,7 +192,7 @@ Besides infinite pipes (both thin and thick), spheres and a couple of other geom

## 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?
But before turning our attention directly 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?

There are literally dozens of ways to numerically solve the equilibrium equations, but for the sake of brevity let us take a look at the three most famous ones. Coincidentally, they all contain the word “finite” in their names. We will not dig into them, but it is nice to know they exist. We might use

@@ -205,7 +205,7 @@ Each of these methods (also called schemes) have of course their own features, p
a. structured, or
b. unstructured

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.
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 just 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 can 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%}\
@@ -216,17 +216,17 @@ Discretisation of a spatial domain
:::::


The first of the three methods is based on approximating derivative (i.e. differentials) by incremental quotients (i.e. differences). The second one heavily relies on geometrical ideas rather than on pure mathematical grounds. Finally, our beloved finite elements are the most “mathematical” ones. Actually, a complete derivation of the finite element method can be written in a textbook without requiring a single figure, just like D’Alembert did more than two centuries ago. In any case, it is important to note that finite differences and elements compute results at the _nodes_ of a mesh, whilst finite volumes compute results at the _cells_ of a mesh. Finally, any method may be used in structured grids but only finite elements and volumes are especially suited for working with unstructured grids.
Back to the three numerical methods, we must say that finite differences is based on approximating derivatives (i.e. differentials) by incremental quotients (i.e. differences). The second one heavily relies on geometrical ideas rather than on pure mathematical grounds. Finally, our beloved finite elements are more strict in the mathematical sense. Actually, a complete derivation of the finite element method can be written in a textbook without requiring a single figure, just like D’Alembert did more than two centuries ago. In any case, it is important to note that finite differences and elements compute results at the _nodes_ of a mesh, whilst finite volumes compute results at the _cells_ of a mesh. Finally, any method may be used in structured grids but only finite elements and volumes are especially suited for working with unstructured grids.

There are technical reasons that justify why the finite element method is the king of mechanical analysis. But that does not mean that other methods may be employed. For instance, fluid mechanics are better solved using finite volumes. And further other combinations may be found in the literature.
There are technical reasons that justify why the finite element method is the king of mechanical analysis. But that does not mean that other methods may be employed. For instance, fluid mechanics are generally better solved using finite volumes. And further other combinations may be found in the literature.

Before proceeding, I would like to make two comments about common nomenclature. The first one is that if we exchanged the words “volumes” and “elements” in all the written books and articles, nobody would note the difference. There is nothing particular in both theories that can justify why finite volumes use volumes and finite elements use elements. Actually volumes and elements are the same geometric constructions. The names were randomly assigned.
Before proceeding, I would like to make two comments about common nomenclature. The first one is that if we exchanged the words “volumes” and “elements” in all the written books and articles, nobody would note the difference. There is nothing particular in both theories that can justify why FVM use “volumes” and FEM use “elements”. Actually volumes and elements are the same geometric constructions. As far as I know, the names were randomly assigned.

The second one is more philosophical and refers to the word “simulation” which is often used to refer to solving a problem using a numerical scheme such as the finite element method. [I am against at using this word to refer to this endeavour](https://www.seamplex.com/blog/say-modeling-not-simulation.html). The term simulation has a connotation of both “pretending” and “faking” something, that is definitely not what we are doing when solving an engineering problem with finite elements. Sure there are some cases in which we simulate, such as using the Monte Carlo method (originally used by Fermi as an attempt to understand how neutrons behave in the core of nuclear reactors). But when solving deterministic mechanical engineering problems I would rather say “modelling” than “simulation.”
The second one is more philosophical and refers to the word “simulation” which is often used to refer to solving a problem using a numerical scheme such as the finite element method. [I am against at using this word for referencing to this endeavour](https://www.seamplex.com/blog/say-modeling-not-simulation.html). The term simulation has a connotation of both “pretending” and “faking” something, that is definitely not what we are doing when we solve an engineering problem with finite elements. Sure, there are some cases in which we simulate, such as using the [Monte Carlo method](https://en.wikipedia.org/wiki/Monte_Carlo_method) (originally used by Fermi as an attempt to understand how neutrons behave in the core of nuclear reactors). But when solving deterministic mechanical engineering problems I would rather say “modelling” than “simulation.”

## Kinds of finite elements {#sec:kinds}

This section is not (just) about different kinds of elements like tetrahedra, hexahedra, pyramids and so on. It is about the different kinds of analysis there are. Indeed, there are a whole plethora of particular types of calculations we can perform, all of which can be called “finite element analysis.” For instance, for the mechanical problem, we can have different kinds of
This section is not (just) about different kinds of elements like tetrahedra, hexahedra, pyramids and so on. It is about the different kinds of analysis there are. Indeed, there is a whole plethora of particular types of calculations we can perform, all of which can be called “finite element analysis.” For instance, for the mechanical problem, we can have different kinds of

* temporal dependence
- steady-state
@@ -255,11 +255,11 @@ And then there exist different pre-processors, meshers, solvers, pre-conditioner
ii. best-estimate
iii. probabilistic

The first one is the easiest because we are allowed to choose parameters and to make engineering decisions that may simplify the computation as long as they give results towards the worse-case scenario. More often than not, an conservative _estimation_ is enough in order to consider a problem solved. Note that this is actually how fatigue results are obtained using fatigue curves, as discussed in\ [@sec:fatigue]. A word of care should be taken when considering what the “worst-case scenario” is. For instance, if we are analysing the temperature distribution in a mechanical part subject to convection boundary conditions, we might take either a very large or a very low convection coefficient as the conservative case. If we needed to design fins to dissipate heat then a low coefficient would be the choice conservative. But if the mechanical properties deteriorated with high temperatures then the conservative way to go would be to set a high convection coefficient. A common practice is to have a fictitious set of parameters, each of them being conservative leading individually to the worst-case scenario even if the overall combination is not physically feasible.
The first one is the easiest because we are allowed to choose parameters and to make engineering decisions that may simplify the computation as long as they give results towards the worst-case scenario. More often than not, a conservative _estimation_ is enough in order to consider a problem as solved. Note that this is actually how fatigue results are obtained using fatigue curves, as discussed in\ [@sec:fatigue]. A word of care should be taken when considering what the “worst-case scenario” is. For instance, if we are analysing the temperature distribution in a mechanical part subject to convection boundary conditions, we might take either a very large or a very low convection coefficient as the conservative case. If we needed to design fins to dissipate heat then a low coefficient would be the conservative choice. But if the mechanical properties deteriorated with high temperatures then the conservative way to go would be to set a high convection coefficient. A common practice is to have a fictitious set of parameters, each of them being conservative leading individually to the worst-case scenario even if the overall combination is not physically feasible.

As neat and tempting as conservative computations may be, sometimes the assumptions may be too biased toward the worst-case scenario and there might be no way of justifying certain designs with conservative computations. It is then time to sharpen our pencils and perform a best-estimate computation. This time, we should stick to the most-probable values of the parameters and even use more complex models that can better represent the physical phenomena that are going on in our problem. Sometimes best-estimate computations are just slightly more complex than conservative models. But more often than not, best-estimates get far more complicated. And these complications come not just in the finite-element model of the elastic problem but in the dependence of properties with space, time and/or temperature, in non-trivial relationships between macro and microscopic parameters, in more complicated algorithms for post-processing data, etc.

Finally, when then uncertainties associated to the parameters, methods and models used in a best-estimate calculation render the results too inaccurate for a certain regulatory body to approve a design, it might be needed to do a full set of parametric runs taking into account the probabilistic distribution of each of the input parameters. It involves
Finally, when then uncertainties associated to the parameters, methods and models used in a best-estimate calculation render the results too inaccurate, it might be needed to do a full set of parametric runs taking into account the probabilistic distribution of each of the input parameters. It involves

1. a thorough analysis of the probability densities of the parameters (and even the methods) of a problem,
2. performing a large number of runs for different combination of parameters, and
@@ -270,9 +270,9 @@ This kind of computation is usually required by the nuclear regulatory authoriti

## Five whys {#sec:five}

So we know we need a numerical scheme to solve our mechanical problem because anything slightly more complex than an infinite pipe does not have analytical solution. We need an unstructured grid because we would not use Legos to discretise pipes. We selected the finite elements method over the finite volumes method, because FEM is the king. Can we pause again and ask ourselves why is it that we want to do finite-element analysis?
So we know we need a numerical scheme to solve our mechanical problem because anything slightly more complex than an infinite pipe does not have analytical solution. We need an unstructured grid because we would not use Legos to discretise cylindrical pipes. We selected the finite elements method over the finite volumes method, because FEM is the king. Can we pause again and ask ourselves why is it that we want to do finite-element analysis?

There exists a very useful problem-solving technique coined by [Taiichi Ohno](https://en.wikipedia.org/wiki/Taiichi_Ohno), the father of the [Toyota production system](https://en.wikipedia.org/wiki/Toyota_Production_System), known as the [Five-whys rule](https://en.wikipedia.org/wiki/5_Whys). It is based on the fact people make decisions following a certain reasoning logic that most of the time is subjective and biased and not purely rational and neutral. By recursively asking (at least five times) the cause of a certain issue, it might possible to understand what the real nature of the problem (or issue being investigated) is. And it might even be possible to to take counter-measures in order to fix what seems wrong.
There exists a very useful problem-solving technique conceived by [Taiichi Ohno](https://en.wikipedia.org/wiki/Taiichi_Ohno), the father of the [Toyota production system](https://en.wikipedia.org/wiki/Toyota_Production_System), known as the [Five-whys rule](https://en.wikipedia.org/wiki/5_Whys). It is based on the fact people make decisions following a certain reasoning logic that most of the time is subjective and biased and not purely rational and neutral. By recursively asking (at least five times) the cause of a certain issue, it might possible to understand what the real nature of the problem (or issue being investigated) is. And it might even be possible to to take counter-measures in order to fix what seems wrong.

Here is an [original example](https://www.toyota-global.com/company/toyota_traditions/quality/mar_apr_2006.html):

@@ -292,14 +292,15 @@ 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](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?
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](https://fs.blog/2012/01/richard-feynman-on-why-questions/) 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 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
Getting back to the case study: do we need to do FEM analysis? Well, it does not look like we can obtain the stresses of the transient cycles 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.
1. to prepare the input data and set up the algorithms,
2. to cope with the many more errors that will inevitable appear during the computation, and
3. 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 cognisant engineer searching and sorting thousands of hard-to-read floating-point numbers.
In the first years of the [history of computers](https://en.wikipedia.org/wiki/History_of_computing_hardware), when programs were 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 particular result 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)

@@ -346,7 +347,7 @@ So we need to address the issue of fatigue in nuclear reactor pipes that

1. are not infinite and have cross-section changes, branches, valves, etc.
2. are made of different materials,
3. are fixed at different locations to the wall through piping supports,
3. are fixed at different locations through piping supports,
4. are subject to
a. pressure transients,
b. heat transients, and
@@ -357,18 +358,18 @@ As I wanted to illustrate in [@sec:five], it is very important to decide what ki
First of all, we should note that we need to solve

i. the heat transfer equation to get the temperature distribution within the pipes,
ii. a frequency analysis of the piping system to get the natural oscillation modes and use them to obtain the pseudo-accelerations created by the design earthquake, and finally
ii. the natural frequencies and oscillation modes of the piping system to obtain the pseudo-accelerations generated by the design earthquake, and finally
iii. the elastic problem to obtain the stress tensor needed to compute the alternating stress to enter into the fatigue curve.

So for each time of the operational transient, the pipes are subject to
So for each time\ $t$ of the operational transient, the pipes are subject to

a. an uniform internal pressure\ $p_i(t)$ that depends on time,
b. a uniform internal temperature $T_i(t)$ that gives rise to a non-trivial time-dependent temperature distribution\ $T(\mathbf{x},t)$ in the bulk of the pipes, and
b. a non-uniform internal temperature $T_i(t)$ that gives rise to a non-trivial time-dependent temperature distribution\ $T(\mathbf{x},t)$ in the bulk of the pipes, and
c. internal distributed forces\ $\mathbf{f}=\rho \cdot \mathbf{a}$ at those times where the design earthquake is assumed to act.

## Thermal transient {#sec:thermal}

Let us invoke our imagination once again. Assume in one part of the transients the temperature of the water inside the pipes falls from say 300ºC down to 100ºC in a couple of minutes, stays at 100ºC for another couple of minutes and then gets back to 100ºC. The temperature within the bulk of the pipes change as times evolves. The internal wall of the pipes follow the transient temperature (it might be exactly equal or close to it through the [Newton’s law of cooling](https://en.wikipedia.org/wiki/Newton%27s_law_of_cooling)). If the pipe was in a state of uniform temperature, the ramp in the internal wall will start cooling the bulk of the pipe creating a transient thermal gradient. Due to thermal inertia effects, the temperature can have a non-trivial dependence when the ramps start or end (think about it!). So we need to compute a real transient heat transfer problem with convective boundary conditions because any other usual tricks like computing a sequence of steady-state computations for different times would not be able to recover these non-trivial distributions.
Let us invoke our imagination once again. Assume in one part of the transients the temperature of the water inside the pipes falls from say 300ºC down to 100ºC in a couple of minutes, stays at 100ºC for another couple of minutes and then gets back to 100ºC. The temperature within the bulk of the pipes changes as times evolves. The internal wall of the pipes follow the transient temperature (it might be exactly equal or close to it through the [Newton’s law of cooling](https://en.wikipedia.org/wiki/Newton%27s_law_of_cooling)). If the pipe was in a state of uniform temperature, the ramp in the internal wall will start cooling the bulk of the pipe creating a transient thermal gradient. Due to thermal inertia effects, the temperature can have a non-trivial dependence when the ramps start or end (think about it!). So we need to compute a real transient heat transfer problem with convective boundary conditions because any other usual tricks like computing a sequence of steady-state computations for different times would not be able to recover these non-trivial distributions.

Remember the main issue of the fatigue analysis in these systems is to analyse what happens around the location of changes of piping classes where different materials (i.e. different expansion coefficients) are present, potentially causing high stresses due to differential thermal expansion (or contraction) under transient conditions. Therefore, even though we are dealing with pipes we cannot use beam or circular shell elements, because we need to take into account the three-dimensional effects of the temperature distribution along the pipe thickness. And even if it we could, there are some tees that connect pipes with different nominal diameters that have a non-trivial geometry, such as the weldolet-type junction shown in\ [@fig:weldolet-cad;@fig:weldolet-mesh]. In this case, there are a number of SCLs (Stress Classification Lines) that go through the pipe’s thickness at both sides of the material interface as illustrated in\ [@fig:weldolet-scls]. It is in these locations that fatigue is to be evaluated.

@@ -395,7 +396,7 @@ Three-dimensional unstructured tetrahedra-based grid for the system shown in\ [@
![Location of the six SCLs defined to analyse fatigue around the junction.](weldolet-scls.png){#fig:weldolet-scls width=75%}


On the one hand, a reasonable number of nodes (remember it is the number of nodes that defines the problem size, not the number of elements) in order to get a decent grid is around 200k for each system. On the other hand, solving a couple of dozens of transient heat transfer problems (which we cannot avoid due to the large thermal inertia of the pipes) during a couple of thousands of seconds over a couple hundred of thousands of nodes might take more time and storage space to hold the results than we might expect.
On the one hand, a reasonable number of nodes (it is the number of nodes that defines the problem size, not the number of elements as discussed in [@sec:elements-nodes]) in order to get a decent grid is around 200k for each system. On the other hand, solving a couple of dozens of transient heat transfer problems (which we cannot avoid due to the large thermal inertia of the pipes) during a few thousands of seconds over a couple hundred of thousands of nodes might take more time and storage space to hold the results than we might expect.

There is a wonderful essay by [Isaac Asimov](https://en.wikipedia.org/wiki/Isaac_Asimov) called [“The Relativity of Wrong”](https://en.wikipedia.org/wiki/The_Relativity_of_Wrong) where he introduces the idea that even if something cannot be computed exactly, there are different levels of error. For instance, believing that the Earth is a sphere is less wrong than believing that the Earth is flat, but wrong nonetheless, since it really deviates from a perfect sphere and resembles more an oblate spheroid.

@@ -419,7 +420,7 @@ 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}(\mathbf{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(\mathbf{x},t)$ is computed by an algebraic generalisation of $\hat{T}(\mathbf{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] 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}(\mathbf{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(\mathbf{x},t)$ is computed by an algebraic generalisation of $\hat{T}(\mathbf{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%}
@@ -430,7 +431,7 @@ Computation of the thermal problem in a reduced mesh and generalisation of the r
:::::


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].
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 second-order elements. Also the grid density is different. Nevertheless, the finite-element solver [Fino](https://www.seamplex.com/fino) used to solve both 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\ [@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

@@ -447,15 +448,15 @@ The material interface in the system from [@fig:real-life] is configured by an o

## Seismic loads {#sec:seismic}

Before considering the actual mechanical problem that will give us the stress tensor at the SCLs and besides needing to solve the transient thermal problem to get the temperature distributions, we need to address the loads that arise due to a postulated earthquake during a certain part of the operational transients. The full computation of a mechanical transient problem using the earthquake time-dependent displacements is off the table for two reasons. First, because the computation would take more time than we might have to deliver the report. And secondly and more importantly, because civil engineers do not compute earthquakes in the time domain but in the frequency domain using the [response spectrum method](https://en.wikipedia.org/wiki/Response_spectrum). Time to revisit our [Laplace transform](https://en.wikipedia.org/wiki/Laplace_transform) exercises from undergraduate maths courses.
Before considering the actual mechanical problem that will give us the stress tensor at the SCLs, and besides needing to solve the transient thermal problem to get the temperature distributions, we need to address the loads that arise due to a postulated earthquake during a certain part of the operational transients. The full computation of a mechanical transient problem using the earthquake time-dependent displacements is off the table for two reasons. First, because the computation would take more time than we might have to deliver the report. And secondly and more importantly, because civil engineers do not compute earthquakes in the time domain but in the frequency domain using the [response spectrum method](https://en.wikipedia.org/wiki/Response_spectrum). Time to revisit our [Fourier transform](https://en.wikipedia.org/wiki/Fourier_transform) exercises from undergraduate maths courses.

### 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 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.
Every nuclear power plant is 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 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](https://en.wikipedia.org/wiki/Diesel_generator) 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](https://en.wikipedia.org/wiki/Nuclear_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].
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, the 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].

![A sample spectrum for a certain floor level of a certain nuclear power plant.](spectrum.svg){#fig:spectrum}
![A sample spectrum for a certain floor level of a certain nuclear power plant.](spectrum.svg){#fig:spectrum width=90%}

### Natural frequencies

@@ -466,12 +467,12 @@ K \phi_i = \lambda_i \cdot M \phi_i
$$
where $K$ is the usual finite-element stiffness matrix, $M$ is the mass matrix, $\lambda_i$ is the $i$-th natural frequency of the structure and $\phi_i$ is a vector containing the nodal displacement corresponding to the $i$-th mode of oscillation.

Practically, these problems are solved using the same mechanical finite-element program one would use to solve a standard elastic problem, provided such program supports these kind of problems ([Fino](https://www.seamplex.com/fino/) does!). There are only two caveats to take into account:
Practically, these problems are solved using the same mechanical finite-element program one would use to solve a standard elastic problem, provided such program supports these kind of problems ([Fino](https://www.seamplex.com/fino/) does). There are only two caveats we need to take into account:

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.
2. The displacement boundary conditions ought to be homogeneous, i.e. only displacements equal to zero can be given. We may fix only one of the three degrees of freedom in certain surfaces and leave the others free, as long as all the rigid body motions are removed as usual.

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.
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. In any case, one is usually interested only in a few of them, namely those with the lower frequencies because they take away most of the energy with them. Each mode has two associated parameters called modal mass and excitation 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%}\
@@ -481,12 +482,12 @@ A real continuous solid has infinite modes of oscillation. A discretised one (us
![$i=4$](mode4.png){width=48%}

![$i=5$](mode5.png){width=48%}\
![$i=6$](mode6.png){width=50%}
![$i=6$](mode6.png){width=48%}

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 millimetre 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 a spatial (actually nodal) distribution of three accelerations (one for each direction) that, when multiplied by the density of the material, 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,16 +499,16 @@ These first modes that take up most of the energy are then used to take into acc
The equivalent accelerations for the piping section of [@fig:modes] for the spectra of\ [@fig:spectrum].
:::::

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.
The ASME code says that these accelerations (depicted in [@fig:acceleration]) ought to be applied twice: once with the original sign and once with all the elements with the opposite sign. Each application should last two seconds.


## Linearity (not yet linearisation) {#sec:linearity}

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.
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 during some seconds a statically-equivalent distributed load arising 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 to exploit the problem’s linearity and compute all the effects separately and then sum them up to obtain the whole combination. We may thus compute just the stresses due to the seismic loads and then add these stresses to the stresses at any time 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 result of the sum of the cases, right? Not always.

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](https://en.wikipedia.org/wiki/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 to weight yourself (let us say you weight 81.2\ kg), then to grab your dog and to 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](https://en.wikipedia.org/wiki/Alpha_particle). Would 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 your favourite FEM program for the first time (remember 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.
Let us now both (i.e. you and me) make an experiment. Grab your favourite FEM program for the first time (remember mine is [CAEplex](https://caeplex.com)) and create a 1mm $\times$ 1mm $\times$ 1mm cube. Set any values for the Young Modulus and Poisson ratio as you want. I chose\ $E=200$\ GPa and\ $\nu=0.28$. Restrict the three faces pointing to the negative axes to their planes, i.e.

* in face “left” ($x<0$), set null displacement in the $x$ direction ($u=0$),
* in face “front” ($y<0$), set null displacement in the $y$ direction ($v=0$),
@@ -519,7 +520,19 @@ Now we are going to create and compare three load cases:
b. Pure shear loads (<https://caeplex.com/p?b494>)
c. The combination of A & B (<https://caeplex.com/p?989>)

The loads in each cases are applied to the three remaining faces, namely “right”, “back” and “top,” and their magnitude in Newtons are:
The loads in each cases are applied to the three remaining faces, namely “right” ($x>0$), “back” ($y>0$ and “top,” ($z>0$). Their magnitude in Newtons are:

````{=latex}
\rowcolors{2}{black!10}{black!0}
````
| | | “right” | | | “back” | | | “top” | |
| ------------------- |:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|
| | $F_x$ | $F_y$ | $F_z$ | $F_x$ | $F_y$ | $F_z$ | $F_x$ | $F_y$ | $F_z$ |
| Case A, pure normal | +10 | 0 | 0 | 0 | +20 | 0 | 0 | 0 | +30 |
| Case B, pure shear | 0 | +15 | -15 | +25 | 0 | -5 | -15 | +25 | 0 |
| Case C, combination | +10 | +15 | -15 | +25 | +20 | -5 | -15 | +25 | +30 |

divert(-1)

````{=latex}
\begin{center}
@@ -544,14 +557,14 @@ $F_z$ \\
\hline

Case A, pure normal & +10 & 0 & 0 & 0 & +20 & 0 & 0 & 0 & +30 \\
Case B, pure shear & 0 & +15 & -15 & +25 & 0 & -5 & -15 & +25 & +30 \\
Case B, pure shear & 0 & +15 & -15 & +25 & 0 & -5 & -15 & +25 & 0 \\
Case C, combination & +10 & +15 & -15 & +25 & +20 & -5 & -15 & +25 & +30 \\
\end{tabular}
\end{center}
````

````{=html}
<table class="table table-responsive">
<table class="table table-responsive w-100">
<tr>
<th></th>
<th colspan="3">face “right” (<span class="math inline">x&gt;0</span>)</th>
@@ -584,7 +597,7 @@ Case C, combination & +10 & +15 & -15 & +25 & +20 & -5 & -15 & +25 & +30 \\
</tr>
</table>
````
divert(0)


In the first case, the principal stresses are uniform and equal to the three normal loads. As the forces are in Newton and the area of each face of the cube is 1\ mm$^2$, the usual sorting leads to
@@ -614,15 +627,15 @@ $$
![Case B, pure-shear loads](cube-shear.png){#fig:cube-shear width=48%}\
![Case C, normal plus shear loads](cube-full.png){#fig:cube-full width=48%}

Spatial distribution of principal stress\ 3 for cases\ B and\ C. If linearity applied, case\ C would be equal to case\ B plus an constant.
Spatial distribution of principal stress\ 3 for cases\ B and\ C. If linearity applied, case\ C would be equal to case\ B plus a constant.
:::::

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 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.
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](https://www.caeplex.com) 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 colour distribution would be _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.
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 long after college, when I had to fac a real engineering problem and not just back-of-the-envelope pencil-and-paper trivial exercises.

Recall that principal stresses are the eigenvalues of the stress tensor. And the fact that in a linear elastic formulation means that the stress tensor of case\ C above is the sum of the individual stress tensors from cases\ A and B does not mean that their eigenvalues can be summed (think about it!). Again, imagine the eigenvalues and eigenvectors of cases A & B. Got it? Good. Now imagine the eigenvalues and eigenvectors for case\ C. Should they sum up? No, they should not! Let us make another experiment, this time with matrices using [Octave](https://www.gnu.org/software/octave/) or whatever matrix-friendly program you want.
Recall that principal stresses are the eigenvalues of the stress tensor. And the fact that in a linear elastic formulation the stress tensor of case\ C above is the sum of the individual stress tensors from cases\ A and B does not mean that their eigenvalues can be summed (think about it!). Again, imagine the eigenvalues and eigenvectors of cases A & B. Got it? Good. Now imagine the eigenvalues and eigenvectors for case\ C. Should they sum up? No, they should not! Let us make another experiment, this time with matrices using [Octave](https://www.gnu.org/software/octave/) or whatever other matrix-friendly program you want.

First, let us create a 3 $\times$ 3 random matrix $R$ and then multiply it by its transpose\ $R^T$ to obtain a symmetric matrix\ $A$ (recall that the stress tensor from [@sec:tensor] is symmetric):

@@ -664,14 +677,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 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 \mathbf{F}$ is\ $\alpha$ times the stress of the beam loaded with a force\ $\mathbf{F}$. Please test this hypothesis by playing with your favourite 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 not to 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 \mathbf{F}$ is\ $\alpha$ times the stress of the beam loaded with a force\ $\mathbf{F}$. Please test this hypothesis by playing with your favourite FEM solver. 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), 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 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.
then the sums (in plural because there are three eigenvalues) 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 should be 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.

\medskip

@@ -680,11 +693,11 @@ The moral of this fable is that if we have a case that is the combination of two

## 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.
After discussing linearity, let us now dig into linearisation. The name is similar but these two animals are very different beasts. We said in\ [@sec:case] that the ASME Boiler and Pressure Vessel Code was born long before modern finite-elements methods were developed and of course being massively available for general engineering analysis (democratised?). Yet the code provides a comprehensive, sound and, more importantly, a widely and commonly-accepted body of knowledge as for the regulatory authorities to require its enforcement to nuclear plant owners. One of the main issues of the ASME code refers to what is known as “membrane” and “bending” stresses, which are defined in [section\ VIII](https://en.wikipedia.org/wiki/ASME_Boiler_and_Pressure_Vessel_Code#ASME_BPVC_Section_VIII_-_Rules_for_Construction_of_Pressure_Vessels) (although widely used in other sections, particularly section\ III) annex 5-A. Briefly, they give the zeroth-order (membrane) and first-order (bending) moments of the stress distribution along a so-called Stress Classification Line or SCL, which should be chosen depending on the type of problem under analysis.

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.
The computation of these membrane and bending stresses are called [“stress linearisation”](https://www.ramsay-maunder.co.uk/knowledge-base/technical-notes/stress-linearisation-for-practising-engineers/) because (I am guessing) it is like computing the Taylor expansion of a real stress distribution along a line, and retaining the first two terms. That is to say, to obtain a linear approximation.

**figures**
dnl **figures**

So what about the SCLs? Well, the ASME standard says that they are lines that go through a wall of the pipe (or vessel or pump, which is what the ASME code is for) from the inside to the outside and ought to be normal to the iso-stress curves. Stop. Picture yourself a stress field, draw the iso-stress curves (those would be the lines that have the same colour in your picture) and then imagine a set of lines that travel in a perpendicular direction to them. Finally, choose the one that seems the prettiest (which most of the time is the one that seems the easiest). There you go! You have an SCL. But there is a catch. So far, we have referred to a generic concept of “stress.” Which of the several stresses out there should you picture? One of the three normals, the three shear, Von\ Mises, Tresca? Well, actually you will have to imagine tensors instead of scalars. And there might not be such a thing as “iso-stress” curves, let alone normal directions. So pick any radial straight line through the pipe wall at a location that seems relevant and now you are done. In our case study, there will be a few different locations around the material interfaces where high stresses due to differential thermal expansion are expected to occur.

@@ -693,7 +706,7 @@ If you cannot wait to know, the expression for computing the $i$-$j$-th element
$$
\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 parametrisation of the SCL. The other linearised stress, namely the _membrane plus bending_ stress tensor \text{MB} again according to ASME VIII Annex 5-A is
where $t$ is the length and $t^\prime$ is the parametrisation variable of the SCL. The other linearised stress, namely the _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
@@ -753,7 +766,7 @@ You can get both the exponential nature of each added bullet and how easily we c
Error in the computation of the linearised stresses vs. CPU time needed to solve the problem using FEM
:::::

## Elements, nodes and CPU
## Elements, nodes and CPU {#sec:elements-nodes}

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_.

@@ -770,17 +783,17 @@ To fix ideas, let us stick to a linear elastic FEM problem. The CPU time needed
3. solving the equations to obtain the displacements, and
4. computing the stress from the displacements.

### Meshing
### Meshing {#sec:meshing}

The effort needed to compute a discretisation 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.

### Building
### Building {#sec:building}

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.

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
### Solving {#sec:solving}

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\ $\mathbf{b}$ containing the loads (both volumetric and the ones at the surfaces from the boundary conditions):

@@ -799,7 +812,7 @@ Structure of the stiffness matrices for the same FEM problem with 10k nodes. Blu

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
### Stress computation {#sec:stress-computation}

In the [displacement-based formulation](http://web.mit.edu/kjb/www/Books/FEP_2nd_Edition_4th_Printing.pdf), the solver just “solves” for the displacements\ $\mathbf{u}(\mathbf{x})$ which are the principal unknowns. But from\ [@sec:tensor] we know that we actually “solved” the problem when we have the stress tensors at every location\ $\mathbf{x}$, which are the secondary unknowns. So the FEM program has to compute the stresses from the displacements using the materials’ [strain-stress constitutive equations](https://en.wikipedia.org/wiki/Constitutive_equation#Stress_and_strain) which involve the Young Modulus\ $E$, the Poisson ratio\ $\nu$ and the spatial derivatives of the displacements\ $\mathbf{u}=[u,v,w]$. This sounds easy, as we (well, the solver) knows what the shape functions are for each element and then it is a matter of computing nine derivatives (that of the three displacements in each of the three directions) and multiplying by something involving\ $E$ and\ $\nu$. Yes, but there is a catch. As the displacements\ $u$, $v$ and\ $w$ are computed at the nodes, we would like to also have the stresses at the nodes. However,


読み込み中…
キャンセル
保存