Bug? Mismatch between 2D and 3D models for steady state heat conduction

1

Hello, I uploaded a new example demonstrating a problem with handling of 2D models in heat transfer.

image
image784Ă—672 5.8 KB
image
image784Ă—672 3.24 KB

There are two versions: left: C3D8 elements, right: CPS4 elements. The red boundary has prescribed temperature, the blue boundary has prescribed heat flux density and there is a volumetric heat production.

image
image784Ă—672 8.15 KB
image
image784Ă—672 7.11 KB

The solid model gives the correct results (for the given discretization), whereas the plane model shows some skewness in the heat flux density. Also, the value of the flux density is not constant in the rightmost element in the non-averaged nodal solution.

1 Like
2

image
image786Ă—704 12.6 KB

Node ordering matters in linear elements…
just reordered your mesh:

*NODE, NSET=Nall
       1,3.000000000000e+00,0.000000000000e+00,0.000000000000e+00
       3,2.000000000000e+00,0.000000000000e+00,0.000000000000e+00
       4,2.000000000000e+00,1.000000000000e+00,0.000000000000e+00
       2,3.000000000000e+00,1.000000000000e+00,0.000000000000e+00
       5,1.000000000000e+00,0.000000000000e+00,0.000000000000e+00
       6,1.000000000000e+00,1.000000000000e+00,0.000000000000e+00
       7,0.000000000000e+00,0.000000000000e+00,0.000000000000e+00
       8,0.000000000000e+00,1.000000000000e+00,0.000000000000e+00
*ELEMENT, TYPE=CPE4, ELSET=Eall
     1,      1,      2,      4,      3
     2,      3,      4,      6,      5
     3,      5,      6,      8,      7

You’ll never see that in a fine mesh, but it happens in coarse meshes. Not a bug. This should be remarked in FEM courses as MacNeal explains in his book (1998):
Finite Elements - Richard MacNeal - Google Libros

1 Like
3

Could you share the entire input for CCX (including all.msh, nam, sur, etc)?

I don’t think node order should matter as long as the surface definitions are adjusted to refer to the correct edges.

Though I would never use 2D thermal. 3D is just as good with shells and OUTPUT=2D if you need to pretend it’s 2D.

4

Yes, that was what I would assume as well. Yet, based on JuanP74’s post I asked our library to order the book by MacNeal, that will take a while because it is out of print and obviously they ignore Google or Amazon. May be I change my mind when reading it. And you can’t ignore that the result indeed depends on the numbering scheme.

And I also agree that it is safest to use generic 3D elements because the reduced dimension versions include MPCs which might cause problems, as seen in other examples of the collection

The following link is valid until 2023-04-26 https://gigamove.rwth-aachen.de/de/download/45de641ddaed40bc7c630fb0d1febb4b

It seems that I can’t attach zip files here.

You can generate the include files by running “cgx -b run.fbd” . It doesn’t harm if separation.py is not executed properly (if you don’t have the helper scripts of the example collection installed), the averaged version is where the problem is most obvious.

5

my apologies, it’s just about reselecting elements to have same boundary conditions as before after renumbering: term1 - Google Drive

6

Sorry, I thought this was commonplace, just most of us forget after such a long time of learning the basics.
The remark in particular is in p. 87 of MacNeal’s book. But is related to the fact that linear shape functions when derived in x or y in-plane variable result in constant deformation (or corresponding variable in case of heat transfer) or linear deformation, and this is not invariant with respect to the element coordinates which are defined by node ordering.

7

@JuanP74, Was that the problem here? Considering it was fixed by re-ordering the nodes of the elements without changing the surfaces? Your link seems to only contain the one that gives symmetric results. I’m looking to recreate the asymmetric one.

That explanation still doesn’t sound right. For a simple square element with isotropic material like this, the differentiation should reduce to the same thing as simply subtracting temperatures on opposite sides and multiplying by k. That’s symmetric between the two element-aligned axes.

@mkraska I can’t get that link to begin downloading :frowning:

8

It’s because of the bilinear shape function for quadrilaterals, for example in the case of displacements in plane elasticity:
u=q1+q2 *x+q3 *y+q4 *x *y
and similar for v

if you derive u with respect to x or y to obtain the strains you’ll see the strains are not invariant with respect to the element orientation.

9

I don’t understand that equation but are you talking about arbitrary orientation or also rotating a square by 90 degrees like we have here? I don’t think there’s any asymmetry in heat flux or strain in that case.

I just tried with S4 shells with heat generation superimposed with a uniform heat flux and it gives exactly the same solution with node numbers rotated by 1 around the elements. So why is S4 OK but not CPS4? I think it’s a bug in CPS4 or a mistake in the model.

10

I apologise by my poor explanations :sweat_smile: I hope Martin can see it and he will be much better than me at explaining it.

11

The fbd file is the original one that creates the files for the “asymmetric” solution.

12

Got CGX to run run2d.fbd but the solution has a symmetric heat flux like the picture you posted, not the one in mkraska’s original.

I hope one of you can share the CCX (not CGX) input files so I can reproduce the problem. I’m not sure my CGX is working right since it shows a lot of error messages but does generate a solution. I think might have some Linux commands that don’t work on Windows??

It’s not so much your explanations as that I can’t believe your conclusion. It’s too contradictory to what I know except in buggy elements or at the edges of physical relevance where numerical error creeps in.

13

It’s funny because I was wondering why for such a simple mesh CGX was creating the ordering of nodes that created the unconverged solution in CCX.
Maybe in a different version the mesh is created in a way the peculiarity doesn’t become apparent.
The explanation is based on the shape functions and the fact that the field of any dependent variable is a function of the spatial coordinates and time. To ensure convergence of the solution the shape functions have to be “complete” but this completeness is only required for convergence when the mesh is refined enough. So, if the solution hasn’t converged then you can’t say the shape functions aren’t complete. This is the case here: depending on the orientation of the elements in the mesh the solution converges for a coarse mesh or not, depending on the boundary conditions.

14

sorry, I posted the 3D model (which has a correct symmetric solution), should have posted the 2D one, the one which exhibits the problem.

This should be the correct one (ccx input, no need to run cgx): https://gigamove.rwth-aachen.de/de/download/a80e7d42f8bf098fbf1cb22788dcb23c

I also get cgx error messages related to postscript output of pathplot. I usually ignore these, as I work with the data files and gnuplot anyways.

Yet back to the original problem. Page 87 isn’t displayed in google books, so I have to wait for my hardcopy… For now I don’t see any reason for anisotropy of the linear 2D element (that is what would be adressed by the different numbering schemes). The shape function is entirely symmetric with respect to x and y, i.e. you can swap directions and get the same shape function.

And we speak of renumbering some nodes and proper adjustment of all node and surface sets. This makes a difference for the y-dependency of the result, which is not ok for boundary conditions which are entirely symmetric.

15

there it goes: therm1_2d - Google Drive

16

Thanks @JuanP74, got it working! Now I’ll try to track this down. @mkraska, that gigamove site’s download button doesn’t work :frowning:

Just to be clear, I don’t think rotating any of the elements’ nodes in any of the 4 possible combinations (ie. 1,2,3,4, 2,3,4,1, 3,4,1,2, 4,1,2,3) should affect the heat flux in the solution, as long as the model is defined consistently with that - ie. loads on the correct faces, any rotated coordinate systems interpreted correctly, etc. Of course making the element wrong won’t work (ie. 4,3,2,1, 1,2,4,3, and friends) but CGX surely doesn’t do that.

Nor should rotating of the whole mesh by an arbitrary angle cause any difference in the solution.

17

Turns out the heat flux load is just applied to the wrong face. Maybe CGX bug since it shows the correct face?

*SURFACE, NAME=SxL
1, S4

S4 means element nodes 4-1 which is global nodes 2,1, which are both at y=0.

1 Like
18

I can’t reproduce the problems with Gigamove. When I click the link in the post, I am redirected to Gigamove and can download the file.

image
image1017Ă—391 15.7 KB

I also confirm that CGX seems to write wrong surfaces.

This is the mesh generated using run2d.fbd.

image

All.msh:

*NODE, NSET=Nall
1,3.000000000000e+00,0.000000000000e+00,0.000000000000e+00
2,2.000000000000e+00,0.000000000000e+00,0.000000000000e+00
3,2.000000000000e+00,1.000000000000e+00,0.000000000000e+00
4,3.000000000000e+00,1.000000000000e+00,0.000000000000e+00
5,1.000000000000e+00,0.000000000000e+00,0.000000000000e+00
6,1.000000000000e+00,1.000000000000e+00,0.000000000000e+00
7,0.000000000000e+00,0.000000000000e+00,0.000000000000e+00
8,0.000000000000e+00,1.000000000000e+00,0.000000000000e+00
*ELEMENT, TYPE=CPS4, ELSET=Eall
1, 1, 4, 3, 2
2, 2, 3, 6, 5
3, 5, 6, 8, 7

This results in the following elements

image

The handbook of ccx shows the face numbers versus local node numbers:
image

In our case, the right edge of the model corresponds to local node numbers 1 and 2 of element 1.
That should be face S1 and not S4 as written to xL.sur.

Contents of xL.sur is

  ** Surfaces based on xL
  *SURFACE, NAME=SxL
  1, S4

When changing that to 1, S1, then the solution is correct. So I confirm the finding by vicmw.

To me this looks like a bug in CGX, I’ll report that to Klaus Wittig.

Thank you all for engaging in the discussion.

1 Like
19

Klaus Wittig has fixed the problem. The fix will be included in the next release of CGX.

The problem was caused by the command flip, which did not update the element face numbers. I used flip to adjust the element normals.

2 Likes