Odd moment curvature results....

[Agent6-6-6]

I don't know if this is a bug or a misunderstanding in what is being reported. But if you review MomentCurvatureResults.n I'd expect that if undertaking an analysis at an axial load of 0 (as this is all that is currently implemented), then for equilibrium these axial loads would be zero as you're looping through all the concrete and steel geometries and adding all the force and moment contributions into the list of moment curvature results.

However this is not the case, you get a whole lot of random numbers like this:-

[0.0, -0.477399625587168, -0.7160994383370962, 0.36906589552563673, -1.042456716982997, -1.418755254032476, -1.676857578247109, 1.7723948539223784, 0.24788216741399083, 1.3830480518536206, 2.277053881669417, 2.392449132697948, 2.1055746989477484, 2.859097843946074, 3.2615732805716107, 2.495338162385451, -4.403265472748899, -5.889407494905754, 6.210181255810312, 7.746934336988488, -7.578399773250567, 10.018343768431805, 5725.446966758347, -11.60338317399146, -11.855866280500777, -13.434843217459274, -12.631551284663146, -15.377393476694124, 13.965917276596883]

If you add a few print statements to the example results in moment_curvature.ipynb examples you should see what I mean.

This seems to indicate convergence is not being achieved, as for the correct moment being reported the axial equilibrium should equal zero with no external axial load being added. Unless I'm misunderstanding here what is being reported?

I cannot seem to track down the issue (as going into the moment curvature code is like going down a rabbithole as it relies on so many other functions for evaluating stresses, strains, etc that I get a bit lost in what is being worked out/returned having not written it myself!).

But it might be related to the assumed centroid not being the geometric centroid as I noted in https://github.com/robbievanleeuwen/concrete-properties/discussions/9#discussioncomment-3779634_, and the fact this can throw out either the moment or axial force equilibrium depending on which one you're evaluating. There probably should probably be a check that the axial equilibrium is equal to the external axial load, otherwise raise an exception?

Also..... When I define my own non-linear stress/strain relationship I tend to get a moment capacity several times the maximum capacity of the section, I'm sure this is related to this issue, but cannot seem to find the culprit and guess solving the above issue my help identify the issue once code base is working as intended. EDIT - sorted this, for future reference - needed to return stress strain relationship to zero in the tension half of the relationship.

Also with the predefined ConcreteLinearNoTension & ConcreteLinear materials, there should be a horizontal part of the curve to the ultimate strain? If you address #47, then you end up with quite different moment curvature response as it is reliant on the defined curve for these classes as it tops out at 0.001 strain.

[Agent6-6-6]

Upon some further investigation, the xtol=1e-3 & rtol=1e-6 variables in the moment_curvature_analysis brentq analysis are potentially the culprit here?

Consider tightening these numbers up a bit or accepting the defaults to improve the convergence?

[robbievanleeuwen]

Hi @Agent6-6-6,

Yes these values are not zero, however noting that this is probably in Newtons assuming your units, I'm happy with numbers like 5 N or 15 N given we are trying to find equilibrium by shifting the neutral axis using brentq. If we were 15 N out when looking for equilibrium given an axial force of 1000 kN I suspect we would accept that given all the inherent variablity in concrete material properties. Although the 5725 N number in the list above seems a bit strange and out of place, is this correct? Tightening the tolerances increases runtime for what I believe to be an unrequired improvement in precision.

EDIT - sorted this, for future reference - needed to return stress strain relationship to zero in the tension half of the relationship.

Correct

Also with the predefined ConcreteLinearNoTension & ConcreteLinear materials, there should be a horizontal part of the curve to the ultimate strain? If you address https://github.com/robbievanleeuwen/concrete-properties/issues/47, then you end up with quite different moment curvature response as it is reliant on the defined curve for these classes as it tops out at 0.001 strain.

I think my comment in #47 addresses this. Currently handled by the parameter ultimate_strain for both and the compressive_strength parameter for ConcreteLinearNoTension. Basically ConcreteLinear is not really intended to be used for a moment-curvature analysis as it's too simple to capture any proper behaviour and should only be used for definining the elastic modulus for area related properties. Perhaps a warning is warranted in the docs for this one?

[Agent6-6-6]

@robbievanleeuwen, yeah the 5kN was a correct output from memory, if I recall correctly I just added a print(res.n) to one of your moment curvature examples in the moment_curvature jupyter notebook to demonstrate the results. i.e.

Hopefully you can reproduce from this example and see if it is something else? Note even if I comment out the rtol/xtol lines to get the default values in the brentq analysis, I still get the following in the same location with the yellow highlighted value, so maybe it is something specific about that particular section/curvature, etc.... every other result is like min 1e-8 convergence accuracy, but that one value is still out by a significant amount by comparison.

Understand 15N not too bad all things considered, but the 5kN thing may pop up and people are not necessarily going to check for it as they expect brentq will converge as no errors are returned.

Perhaps a warning is warranted in the docs for this one?

yeah that would make sense, saying something like stress-strain relationship not intended for section capacity analyses, and only intended for section properties analyses? I guess those particular relationships are not a lot of interest from a design standpoint if you know what you're doing, but you could imagine someone incorrectly using it if they didn't understand what they were doing!

[robbievanleeuwen]

There is something a bit funky going on with that particular example, I will investigate!

[robbievanleeuwen]

Narrowed this one down and fixed in #55.

In summary, using one gauss point to evaluate stresses for non-linear stress-strain profiles is not mesh independent so this has been changed to three gauss points per finite element. Bit of a back story to the above example, but for this particular case the neutral axis was sitting very close to a bar and slight changes in neutral axis changed the way the mesh was being generated, which had a very small impact on the results and the brentq algorithm could not get that close to 0 N.

If you run the above again this should be fixed.