**Moderator:** bernd

why do you have these hard color borders? Is it adjusted by you in Paraview?

Yes I put a threshold value for PEEQ and above that all values are colored red.

As you search to increase the speed of your code, I think it could be a good idea to look at this to solutions :

https://pythran.readthedocs.io/en/latest/ allows you to use a subset of numpy allows vectorization and parallelization very easily.

http://numba.pydata.org/numba-doc/lates ... guide.html allows you to use decorators to optimize your code, so loops are not an issues but must be explicity written, while pythran make automatically the optimization for dot for example.

I think it can produce a code much more readable that cython.

https://pythran.readthedocs.io/en/latest/ allows you to use a subset of numpy allows vectorization and parallelization very easily.

http://numba.pydata.org/numba-doc/lates ... guide.html allows you to use decorators to optimize your code, so loops are not an issues but must be explicity written, while pythran make automatically the optimization for dot for example.

I think it can produce a code much more readable that cython.

Thanks for the lead Jeff. I will have a lookJeffWitz wrote: ↑Mon Apr 29, 2019 1:45 pmAs you search to increase the speed of your code, I think it could be a good idea to look at this to solutions :

https://pythran.readthedocs.io/en/latest/ allows you to use a subset of numpy allows vectorization and parallelization very easily.

http://numba.pydata.org/numba-doc/lates ... guide.html allows you to use decorators to optimize your code, so loops are not an issues but must be explicity written, while pythran make automatically the optimization for dot for example.

I think it can produce a code much more readable that cython.

I am working on plastic collapse of reinforced concrete. The idea is to dimension the reinforcement with the routine now available in version 0.19 and use fcFEM to load the structure to collapse. This will then give an idea of reserve strength and ductility of the structure at a conceptual level, similar to what I did for steel. It will also allow visualisation of the collapse mechanism.

The challenge is how to model concrete, reinforcement and their interaction, while keeping the model simple, robust and easy to use. As I am aiming to capture failure loads, I am less worried about the accuracy of the predicted deformations. The focus is therefore more on the ultimate failure criterion and equilibrium (the ingredients for accurate collapse loads) rather than accurate description of elasticity and softening.

I am currently implementing the models for the concrete matrix and have completed the first option, i.e. a so-called non-associated Drucker-Prager plasticity model with tension cut off. Next on the program is a non-associated Mohr-Coulomb model with tension cut-off.

To show how the Drucker-Prager model behaves, I analysed two simple soil problems. The first is for a uni-axial compression test on cohesive soil with friction (cohesion = 10kPa, phi = 30 degrees) and two values of dilation angle (psi = 30 degrees and psi = 0 degrees):

In this case the dilation angle has no influene on collapse load and the predicted value is within 0.1% of theory:

The deformation at collapse, however, does very much depend on the dilation angle, where for psi = 0 degrees the volumetric strain is zero, whereas for psi = 30 degrees, the material expands. In the ultimate limit state, psi = 0 is more representative of real behaviour, because as a rule of thumb psi = phi - 30.

And finally, the Drucker Prager model applied to the trench of earlier posts:

The collapse loads are significantly higher than predicted without friction. The dilation angle also has a big influence and should therefore be selected conservatively, e.g. psi = phi - 30.

Finally, it should be noted that the collapse load continues to rise because friction causes strength to increase with weight. The determination of collapse loads of frictional materials under weight loading requires a different approach ... something I may come back on at a later stage.

The challenge is how to model concrete, reinforcement and their interaction, while keeping the model simple, robust and easy to use. As I am aiming to capture failure loads, I am less worried about the accuracy of the predicted deformations. The focus is therefore more on the ultimate failure criterion and equilibrium (the ingredients for accurate collapse loads) rather than accurate description of elasticity and softening.

I am currently implementing the models for the concrete matrix and have completed the first option, i.e. a so-called non-associated Drucker-Prager plasticity model with tension cut off. Next on the program is a non-associated Mohr-Coulomb model with tension cut-off.

To show how the Drucker-Prager model behaves, I analysed two simple soil problems. The first is for a uni-axial compression test on cohesive soil with friction (cohesion = 10kPa, phi = 30 degrees) and two values of dilation angle (psi = 30 degrees and psi = 0 degrees):

In this case the dilation angle has no influene on collapse load and the predicted value is within 0.1% of theory:

The deformation at collapse, however, does very much depend on the dilation angle, where for psi = 0 degrees the volumetric strain is zero, whereas for psi = 30 degrees, the material expands. In the ultimate limit state, psi = 0 is more representative of real behaviour, because as a rule of thumb psi = phi - 30.

And finally, the Drucker Prager model applied to the trench of earlier posts:

The collapse loads are significantly higher than predicted without friction. The dilation angle also has a big influence and should therefore be selected conservatively, e.g. psi = phi - 30.

Finally, it should be noted that the collapse load continues to rise because friction causes strength to increase with weight. The determination of collapse loads of frictional materials under weight loading requires a different approach ... something I may come back on at a later stage.

... and here is another good reason why associated plasticity (psi=phi, friction angle equals dilation angle) does not make any sense for analysing frictional collapse:

The green picture in the background was obtained with associated plasticity (psi=phi) and the blue picture in the front with the rule of thumb (psi=phi-30). Clearly the associated plasticity collapse mechanism is unrealistic, with the soil jetting out at mid-height of the trench, while the mechanism obtained with non-associated plasticity (psi=phi-30) shows a credible slip surface. The reason for the difference lies in the volumetric strain at collapse. Shear failure in frictional material with associated plasticity causes the soil to dilate (volume increase), which under plane strain conditions causes the soil to squeeze out.

The green picture in the background was obtained with associated plasticity (psi=phi) and the blue picture in the front with the rule of thumb (psi=phi-30). Clearly the associated plasticity collapse mechanism is unrealistic, with the soil jetting out at mid-height of the trench, while the mechanism obtained with non-associated plasticity (psi=phi-30) shows a credible slip surface. The reason for the difference lies in the volumetric strain at collapse. Shear failure in frictional material with associated plasticity causes the soil to dilate (volume increase), which under plane strain conditions causes the soil to squeeze out.

I implemented a tension cut-off criterion for the Drucker Prager (DP) model. It limits tensile stresses to a user-defined tensile strength. This is a next step towards modelling reinforced concrete (getting closer ). I browsed the literature and find that DP tension cut-offs are usually (if not always) applied to pressure; not principal stress. The simple reason is that the DP model itself is formulated in terms of pressure. A cut-off on principal stress is therefore much more difficult to formulate.

Anyway, IMHO putting a cut-off on pressure makes no sense for concrete. Cracks occur as a result of tensile principal stresses and not pressure. I therefore decided to go the extra mile and formulate a tension cut-off on principal stress as well as on pressure.

Tension cut-off is also relevant for the study of soil bodies. Clay, for example, can accommodate tensile stress under short term loading. The reason is that clay is rather impermeable and water cannot easily drain from it. This means that in clay under tension the water prevents the clay matrix from expanding and therefore carries the tension as suction. In the long term, however, water will flow towards areas of suction and the soil loses its tensile strength. This is the reason cracks appear at the crest of earth dams and embankments.

So lets see how the trench responds with the two versions of tension cut-off (tensile strength = 0.0):

In the foreground the model with tension cut-off on pressure and in the background the model with tension cut off on principal stress. The failure modes are completely different. In the first a normal slip surface develops, but in the second the soil wedge slides and topples into the trench due to lack of tension tying it back to the "hinterland".

The limit loads for the two mechanisms are also completely different. The model with tension cut-off on pressure predicts a limit load which is only 2.9% below the failure load without tension cut-off:

However, the model with tension cut off on principal stress predicts a limit load which is 43% lower !

So a trench which is stable in the short term can easily collapse when drainage occurs.

Next step: reinforcement steel.

Anyway, IMHO putting a cut-off on pressure makes no sense for concrete. Cracks occur as a result of tensile principal stresses and not pressure. I therefore decided to go the extra mile and formulate a tension cut-off on principal stress as well as on pressure.

Tension cut-off is also relevant for the study of soil bodies. Clay, for example, can accommodate tensile stress under short term loading. The reason is that clay is rather impermeable and water cannot easily drain from it. This means that in clay under tension the water prevents the clay matrix from expanding and therefore carries the tension as suction. In the long term, however, water will flow towards areas of suction and the soil loses its tensile strength. This is the reason cracks appear at the crest of earth dams and embankments.

So lets see how the trench responds with the two versions of tension cut-off (tensile strength = 0.0):

In the foreground the model with tension cut-off on pressure and in the background the model with tension cut off on principal stress. The failure modes are completely different. In the first a normal slip surface develops, but in the second the soil wedge slides and topples into the trench due to lack of tension tying it back to the "hinterland".

The limit loads for the two mechanisms are also completely different. The model with tension cut-off on pressure predicts a limit load which is only 2.9% below the failure load without tension cut-off:

However, the model with tension cut off on principal stress predicts a limit load which is 43% lower !

So a trench which is stable in the short term can easily collapse when drainage occurs.

Next step: reinforcement steel.

This mechanism is recognized in practice, as evidenced in the following illustration from a New Zealand government guideline on excavation safety (https://worksafe.govt.nz/dmsdocument/17 ... ion-safety):

I updated the Drucker-Prager model with reinforcement steel. The principle of the model is straightforward: the total strains in the concrete and reinforcement bars are the same, but the elastic and plastic strains (although both adding up to the same value) are different. The increments of stress for each material can then be calculated independently and added to give the total stress increment of the reinforced concrete.

The parameters for the reinforced concrete model are:

Young's Modulus for concrete (32 GPa)

Young's Modulus for steel (210 GPa)

Uniaxial compression strength for concrete (15.75 MPa - EDIT: I made a mistake in the fcFEM input parameters; see next post)

Yield strength reinforcement steel (315 MPa - this includes code material factors)

Friction angle concrete (30 degrees)

Dilation angle concrete (0 degrees)

Sustained tensile strength of concrete (0 MPa)

Reinforcement ratios rhox, rhoy, rhoz (see below, but typically 1-10%)

I applied the model to the analysis of a simply supported beam of 10m length and 0.2x0.6m cross section, loaded by a uniformly distributed load of 10kN/m and self weight.

This shows the elastic deformation obtained with Calculix and the horizontal reinforcement ratio obtained with the model in FreeCAD version 0.19. It is the intent that in the future the calculated reinforcement ratios can be automatically transferred to fcFEM as input for a collapse analysis. For now I have assumed rhox = 0.043, rhoy = 0.008, rhoz = 0.004, which are the maximum values obtained with FreeCAD for this problem.

The analysis with the DP model without tension cut-off, but with reinforcement shows a very high collapse load:

The load multiplier is applied to both weight and external load and indicates a high reserve strength of factor 4.3. Clearly the assumption that the concrete will not crack is invalid and the reason for this very high collapse load.

Subsequently I repeated the analysis with a tension cut-off on pressure, as is customary in the Drucker-Prager model:

The tension cut-off on positive pressure (tension) causes a 23% drop in collapse load. However, as discussed before, the assumption of failure under positive pressure is not realistic for concrete. For example, it would conclude that a stress state of sxx=20MPa (tension), syy=szz=-10MPa (compression) would be acceptable, because it gives to a zero pressure stress state (p=(sxx+syy+szz)/3). It also leads to high volumetric expansion at uniaxial tensile failure, which mobilises the reinforcement in all three directions.

To remedy the shortcoming of the tension cut-off on pressure I reran the analysis with cut-off on principal stress:

The collapse load is now 44% lower than without tension cut-off and 27% lower than the model with tension cut-off on positive pressure. However, there is still a healthy margin between the design load and the ultimate strength of the beam.

Finally, the deformed mesh shows the development of a plastic hinge at the center of the beam:

The parameters for the reinforced concrete model are:

Young's Modulus for concrete (32 GPa)

Young's Modulus for steel (210 GPa)

Uniaxial compression strength for concrete (15.75 MPa - EDIT: I made a mistake in the fcFEM input parameters; see next post)

Yield strength reinforcement steel (315 MPa - this includes code material factors)

Friction angle concrete (30 degrees)

Dilation angle concrete (0 degrees)

Sustained tensile strength of concrete (0 MPa)

Reinforcement ratios rhox, rhoy, rhoz (see below, but typically 1-10%)

I applied the model to the analysis of a simply supported beam of 10m length and 0.2x0.6m cross section, loaded by a uniformly distributed load of 10kN/m and self weight.

This shows the elastic deformation obtained with Calculix and the horizontal reinforcement ratio obtained with the model in FreeCAD version 0.19. It is the intent that in the future the calculated reinforcement ratios can be automatically transferred to fcFEM as input for a collapse analysis. For now I have assumed rhox = 0.043, rhoy = 0.008, rhoz = 0.004, which are the maximum values obtained with FreeCAD for this problem.

The analysis with the DP model without tension cut-off, but with reinforcement shows a very high collapse load:

The load multiplier is applied to both weight and external load and indicates a high reserve strength of factor 4.3. Clearly the assumption that the concrete will not crack is invalid and the reason for this very high collapse load.

Subsequently I repeated the analysis with a tension cut-off on pressure, as is customary in the Drucker-Prager model:

The tension cut-off on positive pressure (tension) causes a 23% drop in collapse load. However, as discussed before, the assumption of failure under positive pressure is not realistic for concrete. For example, it would conclude that a stress state of sxx=20MPa (tension), syy=szz=-10MPa (compression) would be acceptable, because it gives to a zero pressure stress state (p=(sxx+syy+szz)/3). It also leads to high volumetric expansion at uniaxial tensile failure, which mobilises the reinforcement in all three directions.

To remedy the shortcoming of the tension cut-off on pressure I reran the analysis with cut-off on principal stress:

The collapse load is now 44% lower than without tension cut-off and 27% lower than the model with tension cut-off on positive pressure. However, there is still a healthy margin between the design load and the ultimate strength of the beam.

Finally, the deformed mesh shows the development of a plastic hinge at the center of the beam:

I made a mistake in one input value for the concrete beam analyses, which results in a higher uniaxial compressive strength (37.5Mpa instead of the intended 15.75MPa).

Rather than rerunning the analysis, I will here validate the calculated collapse load against an analytical solution based on the higher compressive strength used in the analysis.

Equilibrium shows that the ultimate distributed collapse load of the beam equals qu = 8 Mp / L^2, where qu = collapse load per unit of length, Mp = plastic bending moment of the beam, L = length of the beam.

The plastic bending moment of the cross section at collapse can be derived from equilibrium of forces in the cross section:

where h = height of beam, a = crack height, omega = reinforcement ratio (sorry, I should have used rho), fy = yield strength steel, fck = uniaxial compressive strength of the concrete.

Horizontal equilibrium gives a/h = 0.79 and the resulting force couple gives Mp = 385.5 kNm and hence qu = 8 * 385.5 / 100 = 30.84 kN/m

The actual load on the beam equals q = 10 + 25 * 0.2 * 0.6 = 13 kN/m, where the first contribution is from the external load and the second from weight.

The theoretical load factor at collapse is therefore qu / q = 30.84/13 = 2.37, which is in good alignment with the value calculated with fcFEM for a tension cut off on principal stress (see load-deflection curve in my previous post).

The reason for the excessive reserve strength is the fact that the maximum reinforcement ratio required at the bottom of the beam is applied over the full height. When I complete the functionality for varying the reinforcement ratio by integration point I will rerun the analysis with the minimum required reinforcement following from the elastic analysis.

Rather than rerunning the analysis, I will here validate the calculated collapse load against an analytical solution based on the higher compressive strength used in the analysis.

Equilibrium shows that the ultimate distributed collapse load of the beam equals qu = 8 Mp / L^2, where qu = collapse load per unit of length, Mp = plastic bending moment of the beam, L = length of the beam.

The plastic bending moment of the cross section at collapse can be derived from equilibrium of forces in the cross section:

where h = height of beam, a = crack height, omega = reinforcement ratio (sorry, I should have used rho), fy = yield strength steel, fck = uniaxial compressive strength of the concrete.

Horizontal equilibrium gives a/h = 0.79 and the resulting force couple gives Mp = 385.5 kNm and hence qu = 8 * 385.5 / 100 = 30.84 kN/m

The actual load on the beam equals q = 10 + 25 * 0.2 * 0.6 = 13 kN/m, where the first contribution is from the external load and the second from weight.

The theoretical load factor at collapse is therefore qu / q = 30.84/13 = 2.37, which is in good alignment with the value calculated with fcFEM for a tension cut off on principal stress (see load-deflection curve in my previous post).

The reason for the excessive reserve strength is the fact that the maximum reinforcement ratio required at the bottom of the beam is applied over the full height. When I complete the functionality for varying the reinforcement ratio by integration point I will rerun the analysis with the minimum required reinforcement following from the elastic analysis.