Thick pipe subjected to thermal stress - axial stress inaccuracy in plane strain


I’m analyzing a simple benchmark problem - thick pipe (external diameter 80 mm, internal diameter 60 mm) with a length of 400 mm is made of steel (E=210 [GPa], ν=0.3 [-], λ=16 [mW/mm*K], α=1.2e-5 [1/°C], T_ref=20°C) and subjected to temperature difference - the inner surface has a temperature of 100°C while the outer surface has a temperature of 20°C.

I have an analytical solution for this problem and want to compare it with simulation. I’m using CPE6 plane strain elements (plane strain state is also assumed in the analytical solution) and making use of the symmetry (so that the model is 1/4 of the pipe’s cross-section). The mesh is quite dense. Of course, I applied symmetry boundary conditions and temperature BCs in a coupled thermo-mechanical analysis step. For postprocessing, I’m using ParaView to calculate radial and tangential stresses (those are in good agreement with the analytical solution) and plot the stress values along the path going through the thickness of the pipe.

The problem is that the axial stress obtained in the analysis is very different from the analytical result:

  • simulation:

  • analytical:

As you can see, according to the analytical solution, the axial stress should be somewhere in the range -237 MPa to 37 MPa across the thickness of the pipe while the simulation shows that it varies from -295 MPa to -8 MPa. That’s a big difference, especially since the remaining two stress components are in really good agreement (within a few MPa).

Any ideas what might be causing this discrepancy ?

I attached the input file here:

How are you plotting these stresses along the path in Paraview? Are you doing any transformation w.r.t. the cutting angle?

Thank you for the reply. I’m using formulas to transform the stresses in cartesian coordinates into radial and tangential stress components. But for the axial stress (which is problematic here) I’m just using the SZZ (normal stress in the Z direction) output. To exclude the possibility that the Plot Over Line tool, used to create the path plot in ParaView, is the source of the problem, I checked the values at the selected points on the inner and outer edge. The results are still incorrect there.

your calc. is working with plane strain elements and with axisymmetric elements?


I haven’t tried axisymmetric elements yet. The analytical solution is based on the plane strain assumption and with this approach, I don’t have to worry about BCs like I would have to with a 3D or axisymmetric model.

I submitted the same input file in Abaqus (just had to change the element type to CPE6MT because the same type as the one used in ccx is not available in Abaqus). Interestingly, the results differ from analytical values in a similar manner:

So it’s not a fault of CalculiX.


I have tried with ccx v2.19 , CPE6 plain strain but full model and in two separate steps. (Mecway doesn’t have Thermo mechanical coupling yet). First, I have obtain the thermal profile followed by a transfer of node temperatures to a new step in which I compute the mechanical stress induced by that thermal profile.
I’m obtaining [-249,39] Mpa which doesn’t look bad to be a “direct two step” aproach. Maybe it reflects better the way the theoretical result was obtained.


Thank you for checking this. Your results are in really good agreement with the theory. I will try using uncoupled thermo-mechanical analysis although I think that it should work with coupled as well. Unless some effects occur in coupled simulation and are not accounted for in the theoretical approach…

By full model do you mean no use of symmetry (the whole ring modeled) ? What are the boundary conditions ?

For this kind of problems I like very much constraining rotations in the plane.
That is obtained by imposing displacement in tangential direction (-y,x,0) equal to zero.
With just one BC the whole model is forced to move radially. It’s like a constrain in cilindrical coordinates (Theta=0), but written in Cartesian.

EDITED: 1/4 and conventional symmetry BC should work too. Just tested.

Interesting approach. Exactly how do you realize it ? Do you apply it to all nodes in the model ? Using equation MPC ? Can you share a screenshot showing the settings for this BC in Mecway or its definition in keywords ?

these is a very interesting example!?
maybe it would be nice, to share the hand- and fem calculation!?
maybe Mr. Kraska is interested to put in his collection of examples?

1 Like


If Mr. Kraska is interested in this BC, it would be better to show it’s real potential reducing the model to a simple slice. It is done by means of a transform card. Each node has a new local coordinate system with X’ oriented tangentially in the direction determined by its own coordinates (-y,x,0) and Y’ in the orthogonal/radial (x,y,0). Then a simple constrain of the first axis BC is applied. It is a pure TRANSFORM/rotation of the polar coordinates.
Then you can easily constrain the tangential or radial direction with
1,1, , 0

If your model is not centered just add ((-y+y0,x-x0,0) and (x-x0,y-y0,0) to the transformation.

This can be extended to 3D (Useful to do nonlinear “axisymmetric” which is not available in CalculiX now). There are more complex and amazing BC out there for Parabolic, Spherical, Elliptical,…

I have posted many files with those BC before on the Mecway forum.

NOTE: Victor ¿wouldn’t be nice to implement this thermo-mechanical coupling in Mecway?.You already have all the ingredients.

I tried using the *Uncoupled temperature-displacement step in CalculiX and got the same results as before. That’s really interesting since this type of step should work like a sequence consisting of thermal and mechanical step. Moreover, I performed sequential thermo-mechanical analysis in Abaqus (exactly what you described - thermal analysis first and them mechanical analysis with temperature field imported from the previous simulation). The results are again pretty much the same as before (different than the analytical solution, with even a bit larger difference). I wonder what causes your results to agree with the theory.

This case is taken from the Polish book by M. Bijak-Żochowski although I calculated the analytical solution myself based on the equations given in the book (it agrees with the values shown there for this case). Before publishing anywhere it would be necessary to change the inputs and then reference that book.

Thank you for explaining this. It seems to be really complex. Is this approach based on some literature ?

I’m aware that very good engineers don’t know about this option. I haven’t seen that in any book, but I don’t think it is new.
X=0 and Y=0 axisymmetric condition at the boundaries of a ¼ model is nothing else than two directions where the radial and tangential directions coincide with the Global coordinate system. That’s why is mostly used (I guess). It’s the simplest way of theta=0. The closest I have seen is frictionless support. In fact, I do this because Mecway doesn’t support frictionless for ccx.
It has some minor drawbacks like X=Y=0 singularity , eliminates torsional modes, less visual appeal …

Regarding Paraview and equations to transform,… if you look into any of the global coordinate systems axis, there, SXX and SYY are radial and tangential stress.

Regarding your results,

Problems with mesh or BC: ¿Could you increase “a lot” the deformation scale of your model to see if it is radially growing or some asymmetry has emerged?. Sometimes the expansion is made through a single line, like a “fissure”.

Problems with mesh or BC:¿Could you check that X=0 side has all the nodes on it with coordinates x=0 and Y=0 the same?.

Problem with solver:
¿Could you change Poisson Ratio to 0 to see if it changes something?.

If there is no difference try structured mesh.

That’s right but the transformation gives me plots with radial/tangential stresses that can be read also away from the X=0 and Y=0 edges. Especially since the path used for plotting is inclined at an angle of 45 degrees.

I checked it with a very large deformation scale factor and there are no such spurious effects.

When checked in postprocessor, one side has X coordinates like 10^-15 so essentially zero. The other one is reported as Y=0 everywhere.

With PR=0 I get the axial stress ranging from -248 MPa to -48 MPa.

Hi, Calc_em
Don’t know where the difference could come from.
¿Are reaction forces perfectly perpendicular to the surface?
¿Does the Stresses show a symmetric pattern?
¿Element orientations?
¿Maybe some duplicated or disconected element?

Find attached the liml file. It can be run with the free version. (183 nodes)
Please let me know if you find the source of discrepancy.

I tried with structured meshes (different types of elements - CPE3, CPE4, CPE6 and CPE8) but still no luck.

I haven’t noticed any issues like those that you’ve mentioned.

Thank you for sharing the file. When I saw this shape in one of your screenshots, I thought that it’s an axisymmetric model but now I realized that it’s a sort of unfolded geometry with those tricky BCs that you’ve mentioned. Unfortunately, I can’t utilize this approach for my purposes (it’s supposed to be a simple test example) but thanks to your model I know that it’s possible to achieve the desired range of axial stresses. The only question is what I should do to get such results.

Could you show the liml or inp files with this example involving 1/4 of the pipe’s cross-section and standard symmetry BCs ? It’s the closest one to my approach and still gave you good results so I’d like to compare it with my model.

P.S. I also tried using an axisymmetric model and got the same incorrect axial stresses as always.

That’s the quarter model liml and inp. It is the final mechanical step.
Free version too. 533 nodes.


PD: :sweat_smile:This problem could be solved with 12 elements.When you have time give it a try to those BC’s. It’s worth.


Thank you very much for sharing the files and for your help with this topic.

Problem finally solved. In short, I just had to add:

*Initial conditions, type=temperature
nall, 20

to make my model work correctly (provide accurate axial stresses). Now it works in CalculiX and in Abaqus, with both coupled and sequential approach.

It’s interesting that the influence of this omission on the results is so hidden here. It only highlights the importance of checking all the results and comparing them with analytical solutions whenever possible.

As a side note, it’s worth mentioning that for the sequential approach in CalculiX it’s necessary to add OUTPUT=3D to *NODE FILE. With the default 2D setting CalculiX won’t be able to map temperatures correctly and mechanical analysis will give empty results.

1 Like

I’m glad you were able to resolve it.
Will you share the file or are you planning to do a tutorial?
Maybe you could submit it to that nice youtube channel where they do those nice tutorials.


That’s my YouTube channel actually :wink: And this example is meant for the next video tutorial so you will see it soon :slight_smile: I will share the link once it’s available.

1 Like