Model stops calculating when multiple materials are defined

I have a model here that, in my opinion, behaves a bit strangely. To keep things simple, my model is currently just a rectangular block, meshed with C3D10 tetrahedral elements. There is a fixed constraint on both long sides, and a small surface in the middle for applying forces.

This is my solve.inp script:

*INCLUDE, INPUT=model.inp
*INCLUDE, INPUT=material.inp

*STEP, NLGEOM
*STATIC
*BOUNDARY
fixed, 0, 3, 0
*CLOAD
load, 3, -100
*NODE FILE
U
*EL FILE
S
*END STEP

The “model.inp” file contains the mesh information, including the nodes and element sets.

The “material.inp” file contains the material properties. If I assign the same material to all elements, as shown below, it works perfectly, and the model is solved in about 1–2 seconds.

*MATERIAL, NAME=P_100
*ELASTIC
210000.0, 0.3
*SOLID SECTION, ELSET=P_100, MATERIAL=P_100

In the example above, the element set P_100 contains all elements.

However, if I assign different materials to different elements, the model doesn’t produce a solution. I stopped the solver after 40 minutes. For clarification, here is a shortened version of a modified material.inp file.

*ELSET, ELSET=P_020
1, 5, 6, 14, 16, 17, 19, 22, 23, 24, 25, 29, 30, 31, 32, 38, 
42, 44, 45, 48, 64, 65, 67, 68, 82, 83, 86, 91, 94, 96, 101, 104, 
(...)
*MATERIAL, NAME=P_020
*ELASTIC
42000.0, 0.3
*SOLID SECTION, ELSET=P_020, MATERIAL=P_020
*ELSET, ELSET=P_040
3, 4, 7, 11, 12, 13, 20, 26, 28, 33, 34, 35, 36, 37, 39, 40, 
41, 47, 49, 50, 51, 52, 53, 54, 56, 57, 59, 60, 62, 63, 66, 71, 
*MATERIAL, NAME=P_040
*ELASTIC
42000.0, 0.3
*SOLID SECTION, ELSET=P_040, MATERIAL=P_040
(...)
*ELSET, ELSET=P_100
174, 210, 971, 987, 1276, 1292, 1303, 1430, 1454, 1585, 1686, 1702, 1708, 2166, 2266, 2388, 
2450, 2473, 2489, 2686, 2703, 2765, 2779, 2827, 2992, 2998, 3007, 3053, 3205, 3488, 3937, 4125, 
4263, 4265, 4321, 4450, 4482, 4569, 4671, 4938, 5357, 5411, 5560, 6253, 6254, 6293, 6331, 6332
*MATERIAL, NAME=P_100
*ELASTIC
210000.0, 0.3
*SOLID SECTION, ELSET=P_100, MATERIAL=P_100

I hope this snippet makes the structure clear. I’m assigning different elastic moduli to various elements. The elastic moduli range from that of steel for P_100 (210 GPa) to P_020, which has only 20% of the elastic modulus (42 GPa) of steel. The e-modules are graded in 20% increments, so I have a total of 5 levels.

The question is: Why does the model with a uniform modulus of elasticity compute in a matter of seconds, while the model with different moduli of elasticity (probably) doesn’t compute at all?

P.S.: First of all, I created 100 different materials, each with 1% increments, and ran into the same problem. So I narrowed it down to just 5 different materials because I was worried there might be an upper limit on the number of materials that can be defined. However, I couldn’t find any information on that, and the result is the same. Is there an upper limit on the number of materials that can be defined?

1 Like

Can you share the whole input deck (merged or separate files) for testing ? You could upload it to some hosting website and share the link here.

Multiple materials in the same model certainly do work normally. That’s not not a known bug.

A possibility consistent with your description is that softer materials don’t converge because they’re overloaded but as Calc_em says, we’d have to see the whole deck to work it out.

Hi,

I can definitely share my files. Since I don’t know of any web hosting services, I just went with the first result on Google. Even though it’s just a test model, I still protected the archive with the password “multi_material.”

https://limewire.com/d/OQqoN#rXQRLSMVp1

The range of DOFs should start from 1.

But there’s also one element with missing material assignment - element number 6446.

Young Modulus can be temperature dependant. If you set up an initial temperature profile T=T(x,y,z) you could use just one material definition.The Young modulus will be computed internally. You can also directly define a temperature gradient in beams and shells (I have never used it).

Okay, thanks a lot, everyone. I’ve made a few changes:

  • Corrected the typo involving the degree of freedom 0 (I defenetly switch to much between scripts based on indices of 1 and 0…)
  • The missing element 6446 was the one with the highest number. I guess I didn’t count the “from…to” loop all the way to the end somewhere. But that’s been fixed by the last point.
  • I’m no longer setting the minimum modulus of elasticity quite as low (20%), but only to 50%. I’ll iterate to find the smallest modulus of elasticity at which the simulation runs stably.
  • The idea of a temperature-dependent modulus of elasticity is absolutely brilliant. I will modify the script to include this feature shortly. However, I’ll need to make a few adjustments first. In the end, it might even run more smoothly than it does now. That said, I haven’t run a temperature simulation with Calculix yet, so I’m sure I’ll have a question about that soon. :smile:

In any case, the model is working now. Thanks a million.

This approach is quite common in Abaqus to avoid having to use subroutines for time (e.g. aging) or spatial variation of material properties. Of course, it only makes sense if you don’t simulate actual thermal phenomena. And you don’t need coupled temperature-displacement analysis for that, you can still use *STATIC and *TEMPERATURE to apply ΔT (add *INITIAL CONDITIONS, TYPE=TEMPERATURE too).

there was another issue with ELSET P_060 whose list of elements startet with element label “0”. After elimination of the label “0” your analysis should work:

*ELSET, ELSET=P_060
** original:

0, 2, 15, 21, 27, 43, 46, 55, 69, 70, 73, 80, 93, 95, 97, 99,
** corrected line by elimination of “0”:
**2, 15, 21, 27, 43, 46, 55, 69, 70, 73, 80, 93, 95, 97, 99,
1 Like

Yikes, that’s why element 6446 was missing. I obviously numbered the elements from 0 to N-1 instead of from 1 to N. As I said before, I’ve been switching too much between scripting languages with different indexing schemes.

Thanks for pointing that out!

That’s another reason why CalculiX should have better warning/error handling. Here, it simply failed (or even froze) without any feedback. There is a recent post about it, but it got flagged for some reason: Can’t make sense of CalculiX errors? I built an open-source tool to help diagnose them

Abaqus throws an error about missing section definition, but also warns about that elset:

***WARNING: The following entries for set P_060 are not found and are
removed: 0

Even with the problems corrected, I don’t think that the stress output will be correct since the elements with different material properties are sharing nodes.

The temperature dependence trick is also useful in CCX if you have a huge number of different material property values, like a FGM because there’s a performance problem with too many *MATERIALs.

One could evaluate the integration point stresses, but if it’s a kind of RVE then volume averaged quantities are usually of interest.

Hi all,
from what I see here, you are trying to get something like RVEs to work with your input deck. A similar approach is used for topology optimization. Therfore I highly recommend using following repository with the ccx-executable:

You have to create a file “topOpt_xPhys.txt” with following lines:

1.0
1.0
0.9
0.2

0.1

Where each line N corresponds to the elastic modulus value in percent for the actual element number N of your “uniform” material. This results in 1.0 to be 210e3 and 0.1 in 21e3.

Additionally I would try to get rid of the NLGEOM option, becaus this migh lead to nonconvergence.
Here three examples. The left one is with all “1.0” in the topOpt_xPhys.txt file. The middle one with all “0.5”. The last one has a linear gradient from 0.1 to 1.0 for all Elements from 1 to 6446.

I also tried to use a temperature dependend material in the past, but instead I’ve changed the actual ccx source code. This works quite well - no need to define infinite amount of sets and materials.

Have fun with this!

Cheers,
Stoli

Yes, you’re right. But achieving the correct stress distribution isn’t the goal of my simulation.

I don’t want to simulate an RVE; this is actually an topological optimization attempt. And it’s actually working pretty well now with the temperature-dependent modulus of elasticity. I still need to tweak my code in a few places, but the basic framework seems to be working.

Here are the first 9 iterations of the Cube example. I need to come up with something better than just a block with a surface load. But I’m always so uncreative, and of course I can’t use real components.

animation

2 Likes

Well well, this looks great to me. Finally somebody, who is also working on topology optimization within ccx. Could you describe your framework briefly? Are you using the level-set or the SIMP method? To me, it doesn’t look like BESO.

Here is my attempt (on the very coarse mesh):

output

Currently I’m working on the publication of the results. The code will be available soon.

And to have some more fancy parts … search for “GE Bracket Challenge” or “Airbus Nacelle”. These are typical parts for benchmarking topology optimization algorithms.

Hi stoli,

your animation looks impressive, even with this coarse mesh. What post processor do you use for the animation? (my animation was made with paraview)

I recreated the SKO (Soft-Kill-Option) method from C. Mattheck’s book Design in Nature using CalculiX. Unfortunately, I’m not that well-versed in the theories of topology optimization, so I’d first have to research what “level-set,” “SIMP,” and “BESO” actually are in order to compare how they stack up against the SKO method.

The core idea behind SKO is to iteratively adjust the modulus of elasticity of individual elements based on the local stress state. There are various strategies for doing this, but this is the one I liked best at the moment:

E_n+1 = E_n + k*(sig_n-sig_ref)

with

sig_ref = X*(1-exp(-n/tau)

where X is a target stress level and tau is a constant for numeric stability (smaller values mean faster iterations, but more unstable models.

Thanks for pointing me out to the “GE Bracket Challenge”, I know this geometry but didn’t know the name. Unfortunately, I cant find anything related by google “Airbus Nacelle”.