瀏覽代碼

fixes

tags/CANDIDATE
gtheler 6 年之前
父節點
當前提交
8ced971114
共有 2 個文件被更改,包括 35 次插入30 次删除
  1. +1
    -1
      make.sh
  2. +34
    -29
      nafems4.md

+ 1
- 1
make.sh 查看文件

@@ -145,7 +145,7 @@ for format in ${formats}; do
# post-process
if [[ "${format}" == "html" ]]; then
sed -i 's/<table>/<table class="table table-responsive table-striped">/' ${name}.${format}
sed -i 's/<table>/<table class="table table-responsive-md table-striped table-hover">/' ${name}.${format}
sed -i 's/<blockquote>/<blockquote class="blockquote">/' ${name}.${format}
fi
done

+ 34
- 29
nafems4.md 查看文件

@@ -516,9 +516,9 @@ The ASME code says that these accelerations (depicted in [@fig:acceleration]) ou

## 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 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.
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 a couple of 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 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 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 both 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?

\medskip

@@ -534,7 +534,7 @@ 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” ($x>0$), “back” ($y>0$ and “top,” ($z>0$). 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}
@@ -649,7 +649,7 @@ Now, if we indeed were facing a fully linear problem then the results of the sum

\medskip

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.
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 face 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 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 (try to avoid black boxes as explained in\ [@sec:two-materials]).

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

The computation of these membrane and bending stresses are called [“stress linearisation”](https://www.ramsay-maunder.co.uk/knowledge-base/technical-notes/stress-linearisation-for-practising-engineers/) because (I am guessing) it is like computing the [Taylor expansion](https://en.wikipedia.org/wiki/Taylor_series) (or for the case, in [Legendre polynomials](https://en.wikipedia.org/wiki/Legendre_polynomials)) of a real stress distribution along a line, and retaining the first two terms. That is to say, to obtain a linear approximation. More (optional) mathematical details below.
The computation of these membrane and bending stresses are called [“stress linearisation”](https://www.ramsay-maunder.co.uk/knowledge-base/technical-notes/stress-linearisation-for-practising-engineers/) because (I am guessing) it is like computing the [Taylor expansion](https://en.wikipedia.org/wiki/Taylor_series) (or for the case, expansion in [Legendre polynomials](https://en.wikipedia.org/wiki/Legendre_polynomials)) of an arbitrary stress distribution along a line, and retaining the first two terms. That is to say, to obtain a linear approximation. More (optional) mathematical details below.

dnl **figures**

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

\medskip

Now the optional (but recommended) mathematical details. Accoding to ASME\ VIII Div.\ 2 Annex\ 5-A, the expression for computing the $i$-$j$-th element of the membrane tensor is

$$
\text{M}_{ij} = \frac{1}{t} \cdot \int_0^t \sigma_{ij}(t^\prime) \, dt^\prime
$$
where $t$ is the length and $t^\prime$ is the parametrisation variable of the SCL. The other linearised stress, namely the _bending_ stress tensor $\text{B}$ is
where $t$ is the length and $0<t^\prime<t$ is the parametrisation variable of the SCL. The other linearised stress, namely the _bending_ stress tensor $\text{B}$ is

$$
\text{B}_{ij} = \frac{6}{t^2} \cdot \int_0^t \sigma_{ij}(t^\prime) \cdot \left( \frac{t}{2}-t^\prime\right) \, dt^\prime
$$

For the fatigue assessment, it is the sum\ $\text{MB}$ which is taken as a measure of the stress
For the fatigue assessment, it is the sum\ $\text{MB}$ that measures the stress

$$
\text{MB}_{ij} = \text{M}_{ij} \pm \text{B}_{ij}
@@ -740,7 +742,7 @@ In any case, $\text{MB}_{ij}$ should be taken as a measure of the primary stress
$$
\sigma(t^\prime) = a\cdot t^\prime + b
$$
because any other higher-order polynomial term would self-balance.
because any other higher-order polynomial term would self-balance (exercise: prove it).
The membrane and bending stresses, according to ASME would then be

$$
@@ -750,7 +752,7 @@ $$
\text{B} = \frac{6}{t^2} \int_{0}^{t} \big[ a\cdot t^\prime + b \big] \cdot \left( \frac{t}{2} - t^\prime \right) \, dt^\prime = -\frac{1}{2} \cdot a \cdot t
$$

We can see that
In effect, we can see that
$$
\sigma(0) = b = \text{M} + \text{B}
$$
@@ -789,21 +791,21 @@ Two of the hundreds of different ways the infinite pressurised pipe can be solve
:::::


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

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

::::: {#fig:error-vs-cpu}
@@ -811,16 +813,18 @@ You can get both the exponential nature of each added bullet and how easily we c

![Membrane plus bending $\text{MB}$](error-MB-vs-cpu.svg){#fig:error-MB-vs-cpu width=90%}

Error in the computation of the linearised stresses vs. CPU time needed to solve the problem using FEM
Error in the computation of the linearised stresses vs. CPU time needed to solve the infinite pipe problem using the finite element method
:::::

An additional note should be added. The FEM solution, which not only gives the nodal displacements but also a method to interpolate these values inside the elements, does not fully satisfy the original equilibrium equations at every point (i.e. the strong formulation). It is an approximation to the solution of the [weak formulation](https://en.wikipedia.org/wiki/Weak_formulation) that is close (measured in the vector space spanned by the [shape functions](https://www.quora.com/What-is-a-shape-function-in-FEM) to the real solution. Mechanically, this means that the FEM solution leads only to nodal equilibrium but not pointwise equilibrium.
An additional note should be added. The FEM solution, which not only gives the nodal displacements but also a method to interpolate these values inside the elements, does not fully satisfy the original equilibrium equations at every point (i.e. the strong formulation). It is an approximation to the solution of the [weak formulation](https://en.wikipedia.org/wiki/Weak_formulation) that is close (measured in the vector space spanned by the [shape functions](https://www.quora.com/What-is-a-shape-function-in-FEM)) to the real solution. Mechanically, this means that the FEM solution leads only to nodal equilibrium but not pointwise equilibrium.

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

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_.
The last two bullets above lead to an issue that has come many times when discussing the issue of convergence with respect to the mesh size with other colleagues. There apparently exists a common misunderstanding that the number of elements is the main parameter that defines how complex a FEM model is. This is strange, because even in college we are taught that the most important parameter is the _size_ of the stiffness matrix, which is three times (for 3D problems with the [displacement-based formulation](http://web.mit.edu/kjb/www/Books/FEP_2nd_Edition_4th_Printing.pdf) formulation) the number of _nodes_.

Let us pretend we are given the task of comparing two different FEM programs. So we solve the same problem in each one and see what the results are. I have seen many times the following situation: the user loads the same geometry in both programs, run the meshing step in both of them so that the number of elements is more or less them same (because she wants to be “fair”) and then solves the problem. Voilà! It turns out that the first program defaults to first-order elements and the second one to second-order elements. So if the first one takes one minute to obtain a solution, the second one should take nearly four minutes. How come that is a fair comparison? Or it might be the case that one program uses tetrahedra while the other one defaults to hexahedra. Or any other combination. In general, there is no single problem parameter that can be fixed to have a “fair” comparison, but if there was, it would definitely be the number of _nodes_ rather than the number of _elements_. Let us see why.

Let us pretend we are given the task of comparing two different FEM programs. So we solve the same problem in each one and see what the results are. I have seen many times the following situation: the user loads the same geometry in both programs, run the meshing step in both of them so that the number of elements is more or less them same (because she wants to be “fair”) and then solves the problem. Voilà! It turns out that the first program defaults to first-order elements and the second one to second-order elements. So if the first one takes one minute to obtain a solution, the second one should take nearly four minutes. How come that is a fair comparison? Or it might be the case that one program uses tetrahedra while the other one defaults to hexahedra. Or any other combination. In general, there is no single problem parameter that can be fixed to have a “fair” comparison, but if there was, it would definitely be the number of nodes rather than the number of elements. Let us see why.
\medskip

Fire up your imagination again and make a [thought experiment](https://en.wikipedia.org/wiki/Thought_experiment) in which you have to compare say a traditional FEM approach with a new radical formulation that a crazy mathematician from central Asia came up with claiming it is a far superior theory than our beloved finite elements (or for the case, any other formulation from\ [@sec:formulations]). How can we tell if the guy is really a genius or purely nuts? Well, we could solve a problem which we can compute the analytical solution (for example the infinite pipe from\ [@sec:infinite-pipe]) first with the traditional method ([@sec:infinite-pipe-fem]) and then with the program which uses the new formulation. Say the traditional FEM gives an error between 1% and 5% running in a few seconds depending on the mesh size. The new program from the crazy guy takes no input parameters and gives an error of 0.1%, but it takes one week of computation to produce a result. Would you say that the new radical formulation is really far superior?

@@ -835,13 +839,13 @@ To fix ideas, let us stick to a linear elastic FEM problem. The CPU time needed

### 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.
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 (3D). 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 {#sec:building}

The [stiffness matrix]([stiffness matrix](https://en.wikipedia.org/wiki/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.
The [stiffness matrix]([stiffness matrix](https://en.wikipedia.org/wiki/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. This is called the assembly of the 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!). 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.
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!). So the sweep is larger for linear grids. But the effort needed to integrate quadratic shape functions is greater than for the linear case, so these two effects almost cancel out.

### Solving {#sec:solving}

@@ -863,11 +867,11 @@ So the question is... how hard is to solve a sparse linear problem? Well, the nu
Structure of the stiffness matrices for the same FEM problem with 10k nodes. Red (blue) are positive (negative) elements
:::::

In a similar way, different types of elements will give rise to different sparsity patterns which change the effort needed to solve the problem. In any case, the base parameter that controls the problem size and thus provides a basic indicator of the level of difficulty the problem poses is the number of nodes. And again, not the number of elements.
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. Again, not the number of elements, as the solver does not even know if the matrix comes from FEM, FVM or FDM.

### Stress computation {#sec:stress-computation}

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

i. the displacements\ $\mathbf{u}(\mathbf{x})$ are not differentiable at the nodes, and
ii. if the node belongs to a material interface, neither\ $E$ nor\ $\nu$ are defined.
@@ -948,15 +952,15 @@ What we as responsible engineers who have to sign a report stating that a nuclea

Without resorting into philosophical digressions about the difference between [free and open-source software](https://en.wikipedia.org/wiki/Free_and_open-source_software) (not because it is not worth it, but because it would take a whole book), the programs that make their source code available for their users are called [open-source software](https://en.wikipedia.org/wiki/Open-source_software). If the users can also modify and re-distribute the modified versions, they are called [free software](https://www.fsf.org/about/what-is-free-software). Note that the important concept here is freedom, not price. In Spanish (my native language) it would have been easier because there are two separate words for free as in freedom (“libre”) and for free as in price (“gratis”).

In effect, a couple of years ago Angus Ramsay noted [a weird behaviour](https://www.ramsay-maunder.co.uk/knowledge-base/technical-notes/asmeansys-potential-issue-with-thermal-expansion-calculations/) in the results given by a certain commercial [non-free](https://en.wikipedia.org/wiki/Proprietary_software) FEA software regarding the handling of expansion coefficients from ASME data. To understand what was going on, we had to guess what the program was doing to [reproduce the allegedly weird results](https://www.seamplex.com/docs/SP-WA-17-TN-F38B-A.pdf). Finally, it was a [matter how the data was rounded up to be presented in a paper table](https://www.ramsay-maunder.co.uk/knowledge-base/technical-notes/accuracy-of-thermal-expansion-properties-in-asme-bpv-code/) rather than a programming flaw. Nevertheless, we were lucky our guesses lead us to a reasonable answer. If we had access to the program’s source code, we could have thoroughly analysed the issue in a more efficient way. Sure, we might not have the same programming skills the original authors of the software have, but if it had been [free software](https://en.wikipedia.org/wiki/Free_software) we would have had the _freedom_ to hire a programmer to help us out. That is what _free_ means. In [Eric Raymond](https://en.wikipedia.org/wiki/Eric_S._Raymond)’s words, [“given enough eyeballs, all bugs are shallow.”](http://www.catb.org/~esr/writings/cathedral-bazaar/) This is rather important in engineering software where [verification and validation](https://en.wikipedia.org/wiki/Software_verification_and_validation) is a must, especially in regulated fields like the nuclear industry. First, think how can a piece of software be verified if the source code is not available for independent analysis. And then, ask yourself another question:
In effect, a couple of years ago Angus Ramsay noted [a weird behaviour](https://www.ramsay-maunder.co.uk/knowledge-base/technical-notes/asmeansys-potential-issue-with-thermal-expansion-calculations/) in the results given by a certain commercial [non-free](https://en.wikipedia.org/wiki/Proprietary_software) FEA software regarding the handling of expansion coefficients from ASME data. To understand what was going on, Angus and I had to guess what the program was doing to [reproduce the allegedly weird results](https://www.seamplex.com/docs/SP-WA-17-TN-F38B-A.pdf). Finally, it was a [matter how the data was rounded up to be presented in a paper table](https://www.ramsay-maunder.co.uk/knowledge-base/technical-notes/accuracy-of-thermal-expansion-properties-in-asme-bpv-code/) rather than a programming flaw. Nevertheless, we were lucky our guesses lead us to a reasonable answer. If we had access to the program’s source code, we could have thoroughly analysed the issue in a more efficient way. Sure, we might not have the same programming skills the original authors of the software have, but if it had been [free software](https://en.wikipedia.org/wiki/Free_software) we would have had the _freedom_ to hire a programmer to help us out. That is what _free_ means. In [Eric Raymond](https://en.wikipedia.org/wiki/Eric_S._Raymond)’s words, [“given enough eyeballs, all bugs are shallow.”](http://www.catb.org/~esr/writings/cathedral-bazaar/) This is rather important in engineering software where [verification and validation](https://en.wikipedia.org/wiki/Software_verification_and_validation) is a must, especially in regulated fields like the nuclear industry. First, think how can a piece of software be verified if the source code is not available for independent analysis. And then, ask yourself another question:

> Do you trust your favourite FEM program?

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

## A parametric tee {#sec:parametric}

Time for another experiment. We know (more or less) what to expect from an infinite pressurised pipe from\ [@sec:infinite-pipe;sec:infinite-pipe-fem]. What if we added a branch to such pipe? Even more, what if successively varied the diameter of the branch to see what happens? This is called parametric analysis, and sooner or later (if not now) you will find yourself performing this kind of computations more often than any other one.
Time for another experiment. We know (more or less) what to expect from an infinite pressurised pipe from\ [@sec:infinite-pipe] and [@sec:infinite-pipe-fem]. What if we added a branch to such pipe? Even more, what if we successively varied the diameter of the branch to see what happens? This is called parametric analysis, and sooner or later (if not now) you will find yourself performing this kind of computations more often than any other one.

So here come the five Feynmann-Ohno questions:

@@ -1019,7 +1023,7 @@ Parametric membrane and bending stresses as a function of the nominal diameter\
Von\ Mises stress and 400x warped displacements for three values of\ $d_b$.
:::::

[@Fig:tee-post] illustrates how the pipes deform when subject to the internal pressure. When the branch is small, the problem resembles the infinite-pipe problem where the main pipe expands radially outward and there is only traction. For large values of\ $d_b$, the pressure in the branch bends down the main pipe, generating a complex mixture of traction and compression. The tipping point seems to be around\ $d_b\approx 5$\ in.
[@Fig:tee-post] illustrates how the pipes deform when subject to the internal pressure. When the branch is small, the problem resembles the infinite-pipe problem where the main pipe expands radially outward and there is only traction. For large values of\ $d_b$, the pressure in the branch bends down the main pipe, generating a complex mixture of traction and compression. The tipping point seems to be around a branch diameter\ $d_b\approx 5$\ in.

Do you now see the added value of training throw-ins with watermelons? We might go on...

@@ -1029,9 +1033,9 @@ Do you now see the added value of training throw-ins with watermelons? We might
* using distributed forces from earthquakes,
* etc.

Most of the time at college we would barely do what is needed to be approved in one course and move on to the next one. If you have the time and consider a career related to finite-element analysis, please do not. Clone the repository from <https://bitbucket.org/seamplex/tee> and start playing. If you are stuck, do not hesitate asking for help in [wasora’s mailing list](https://www.seamplex.com/lists.html). It will pay off later on.
Most of the time at college we would barely do what is needed to be approved in one course and move on to the next one. If you have the time and consider a career related to finite-element analysis, please do not. Clone the repository^[<https://bitbucket.org/seamplex/tee>] with the input files for [Fino](https://www.seamplex.com/fino) and start playing. If you are stuck, do not hesitate asking for help in [wasora’s mailing list](https://www.seamplex.com/lists.html).

One further detail: it is always a sane check to try to explain the numerical results based on physical reasoning as we did two paragraphs above. Most of the time you will be solving problems whilst already knowing what the result would (or ought to) be.
One further detail: it is always a sane check to try to explain the numerical results based on physical reasoning (i.e. “with your fingers”) as we did two paragraphs above. Most of the time you will be solving problems whilst already knowing what the result would (or ought to) be.

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

@@ -1300,9 +1304,10 @@ A list of the tips that arose throughout the text:
* follow the “five whys rule” before compute anything, probably you do not need to
* use engineering judgment and make sure understand the [“wronger than wrong”](https://en.wikipedia.org/wiki/Wronger_than_wrong) concept
* play with your favourite FEM solver (mine is [CAEplex](https://caeplex.com)) solving simple cases, trying to predict the results and picturing the stress tensor and its eigenvectors in your imagination
* try to show that even though the principal stresses are not linear with respect to summation, they are linear with respect to multiplication (with pencil and paper)
* grab any stress distribution from any of your FEM projects, compute the iso-stress curves and the draw normal lines to them to get acquainted with SCLs
* first search online for “stress linearisation” (or “linearization” if you want) and then get a copy of ASME\ VIII Div\ 2 Annex 5-A
* try to show that even though the principal stresses are not linear with respect to summation, they are linear with respect to multiplication (with pencil and paper)
* remember that FEM solutions lead only to nodal equilibrium but not pointwise equilibrium
* measure the time needed to generate grids of different sizes and kinds with your favourite mesher
* learn this by heart: the complexity of a FEM problem is given mainly by the number of _nodes_, not by the number of elements
* remember that welded materials with different thermal expansion coefficients may lead to fatigue under cyclic temperature changes

Loading…
取消
儲存