Refactoring and Revising State Determination Process and Force Increment Calculation

[z-ichinohe]

Hello, Mr. Ahmed.

In order to deal with the problem I encountered, I refactored the code and revised state determination process and force increment calculation. Almost all of the strength deterioration process is preserved. The previous code was a bit redundant for me. I think lighter and clearer code will be helpful for fix in future, and state determination and force calculation process also should be simple. Unfortunately I couldn't specify all source of the bugs but the new code works well for me in both static and dynamic analysis. It lacks tests for a spring with residual strength or ultimate capacity at this time. I also lacked information about your coding principles. Could you help me refine the codes and contribute to the community? Of course, your question, opinion and modifications are welcome. Ask me whatever you notice since I'm not experienced well to create a pull request.

[z-ichinohe]

This is a result from the revised code in contrast with the figure attached to the issue. After Before

[amaelkady]

@z-ichinohe Thank you for the modifications. I'm currnetly on leave and will be back in the office on August 11th. Once I'm back I will review your code. Best regards,A.

[z-ichinohe]

I've finished more thorough refactoring and reconsideration of variables and working on pinching behaviors. Then, this is what I got for the strain history and the parameters of the issue raised by Ms. Mazzoni.

[amaelkady]

@z-ichinohe hi Kazuki, this looks very good. I have been modifying the code of the original models to address some of the raised issues. I'll continue to run more tests. Meanwhile, it would be good that you check your model with various deformation histories to make sure it is robust. Also, try to briefly document the conecpt of your modifications; particularly the "branch" flags.

[z-ichinohe]

My concepts are below. I shall add some illustrations to demonstrate how branches works and additional modifications to implement these concepts thoroughly. It would be appreciated if you have a displacement set to confirm the robustness and share them with me.

Concept of Modifications

2022 Sep 13 Kazuki ICHINOHE

Refactoring

Omitted variables

ui_1
du
Epj
sPCsp
Excursion_Flag

Variables which don't preserve their values to the next step are removed from header files and declared in the source file.

_j_1 _0

Up_pos              -> posUp_0
Uy_neg_j_1          -> negUy
Kul_j_1             -> Kunload

"_j_1" are removed and initial values are distinguished with "_0".

pos neg eng

Up_pos              -> posUp_0
Uy_neg_j_1          -> negUy
Energy_Acc          -> engAcml

"pos", "neg" and "eng" are moved to the front because it just looks prettier.

max -> cap

FmaxFy_pos          -> posFcapFy_0

"max" are substituted with "cap" because it may be confused with "peak".

local global

Upeak_pos_j_1       -> posUglobal
ULastPeak_pos_j_1   -> posUlocal

You have mentioned "peak" as "global peak" and "local" may be matched with "global" though "the last peak" is an exact expression for the implementation.

State Determination

"Branch" is an advanced "Flag"s and straight-forward way to express the "CASE"s. Reversal_Flag is raised only when unloading starts, while Unloading_Flag is raised until the unloading ends. Reversal_Flag means an occurrence of the event which needs calculating the stiffness deterioration, but Unloading_Flag means the branch where it is. Moreover Unloading_Flag and Reloading_Flag would never be raised at the same time. Then, we can use an integer to express where it is, which is "Branch" and the flags for whether updating deterioration is necessary or not.

The last branch where it was is very helpful for the state determination. When it's on the elastic branch, the next possible branch will be one of the two branches: the positive backbone branch towards the capping point and the negative one. We should only consider whether it's beyond the yield point or not. If it is on the same branch as the last, the force increment will just be calculated by multipling the displacement increment and the tangent stiffness. It sometimes need more lines than other methods but it's more straight-forward to express our perception of the hysteresis model.

Branch

0:  Elastic
1:  Unloading Branch
2:  Towards Local Peak      +
3:  Towards Pinching Point  +
4:  Towards Global Peak     +
5:  Towards Capping Point   +
6:  Towards Residual Point  +
7:  Residual Branch         +

0 and 1 are (pseudo-)elastic, 2-4 are reloading phases and 5-7 are on the backbone curve. 2-7 have their counterparts 12-17 in negative.

Branch Shift Summary

0   ->         5,              15
1   ->   2,3,4,       12,13,14
2   -> 1,  3,4                
3   -> 1,    4
4   -> 1,      5
5   -> 1,        6
6   -> 1,          7
7   -> 1

Considering concecutive shift consecutive "if" clauses instead of "else if" clauses is desirable, which has not been thoroughly covered yet.

1 to 12-14 are accompanied by Excursion_Flag and 2-7 to 1 are accompanied by Reversal_Flag.

If branch is unchanged, the force increment will just be calculated as a product of the displacement increment and the tangent stiffness. If not,the tangent stiffness is updated and the force increment calculated respectively.

Others

• Failure or zero displacement increment should be caught at the beginning.
• Ambiguous evaluation with 1.05 or 1.01 should not be used.
• Tangent stiffnesses should be and can be calculated accurately.
• If pinching, the local peak should be compared with the pinching point.

[amaelkady]

@z-ichinohe Thanks for the clarification. These updates make sense and for sure seem better than the original code.

** For the displacement protocols, you should consider the assymetric protocol; see example below set disp [list 0 2 -2 2 -2 5.5 1.5 4.5 1.5 5.5 1.5 3.5 5.5 1.5 5.5 1.5 9 5 8 5 9 5 7];

** For any protocol that you use, try different levels of resolution (i.e., disp increment size "du"). Many issues typically appear with too large increments at backbone pivot points; see example below:

** Check the response with the residual strength branch because several issues can be traced to this part.

** Finally, there are several special cases that need to be addressed by the code (see attached document). You can investigate these cases with custom protocols. These cases are perhaps what made the current code a bit busy but may be this will not be an issue with your new branch approach. Cases_Pinching.pdf

Let me know if you have questions. A.

[z-ichinohe]

Thank you for your precise suggestion! This loading protocol is designed to cover almost all of the branch shifts but one-step unloading completion is not included.

Peak-Oriented

First of all, though Λ is 1000 contrary to the yield displacement 1 as cyclic deteriorations are ignored your functions have some constant deteriorations. After unloading at the point 15, as the local peak has been updated and it goes towards the global peak. It may look strange but preserving all local peaks is unrealistic.

Pinching

Pinching behaviors towards negative are different. I think the pinching displacement will be zero at the first excursion. At the loop 20-21-22-23-24, I go towards the global peak, and you go towards the local one. Yours is right. I should add a possible branch shift from 3(Pinching) to 2(Local Peak). However, at the loop 24-25-26-27-28, you go towards the local peak again but at the last, I can't guess why it happened. Moreover, it looks like lacking the continuity across the unloading 29. At the last, yours successfully caught the capacity loss in pinching one, but failed in peak-oriented one.

Reloading towards the peak is one of the complicated behaviors. The pdf file you attached will be helpful for me to understand that, but I think only p1 is correct. p2 and p3 should go towards the "Global peak" and p4 will go to the pinch point in my opinion. Could you add some explanations to the file?

[amaelkady]

@z-ichinohe Thank you for checking:

For the loading cases in the PDF files:

**Case in P2, the dashed lines represent the different options to move to 9local vs global peak). I was using these illutrations back then to discuss with other researchers the proper branch to be used. In that case, i beleive loading should be towards the global peak (i think this is also what the original code does).

** Case in P3, is a bit trickey, I believe it makes sense -physically- that minor unloading should not affect the original pathtowards the pinch point before moving to the global peak. This is important to maintain the pinching behvaiour.

** With repect to cyclic deterioration, it is strange that there is a difference between the two models. You are using the same energy rules and formulas as the original code, correct?

** Also, can you send me the model definition you are using in the figures above as well as the disp history.

[z-ichinohe]

These are a test script written in Python and a JSON file including parameters and displacement history. Python file works alone and JSON file is a mere supplementary. Please download and modify their extensions since github doesn't support .py or .json to attach. To be honest, I made some modifications to the deterioration section, but it doesn't change the logic. I've just introduced some += or *=, and replace "max" or "min" with if clauses. imktest.txt imktestp.txt

[z-ichinohe]

This is a result of 0.01 displacement increment while the figure above was 0.25. Mine is also altered to meet the criteria below. I noticed that the "constant deterioration behavior" is caused by very sparse displacement increment because you don't make the point on the skeleton curve thoroughly when the branch changes. They only have slight differences made by that.

[z-ichinohe]

Pinching Behavior Summary

Though I have not understood the PDF file you attached yet, I summarized the pinching behavior varieties affected by the local peak location(L) and the excursion starting point(E). They are drawn as elastic-perfectly plastic and the reloading towards negative might not be correct but for the sake of brevity. The global peak(G) deterioration is ignored except for case 1. In case 1-4, the unloading starts after passing the pinching point(P) and otherwise in case 5-8. The excursion starting points(E) move towards negative from case 2 to case 4 and from case 5 to case 8 except for case 7.

Case 1: The deteriorated global peak is smaller than the local peak.

If the pinching point or the local peak is not in the square between the E and G, ignore them and move to the global peak directly.

Case 2: The pinching point is behind the excursion starting point.

If P or L is in the square, go through that.

Case 3: The local peak locates between the pinching point and the global peak.

The orange line is correct.

Case 4: The local peak locates between the pinching point and the global peak but beneath the line from P to G.

The blue one is correct.

Case 5: The new pinching point is lower than the local peak.

Case 6: The new pinching point is higher than the local peak.

I think you chose the blue line in both of case 5 and 6 while I chose the orange line in case 6. However, I agree with you now. The orange line in case 6 looks so strange.

Case 7: Minor unloading

This is why I add the orange line in case 5 and the green line in case 6. If it crosses the x-axis, it will move on the blue but move on the orange if not. It's annoying to me that a slight difference affects much. However, it might not be avoidable.

Case 8: The new pinching point is higher than the local peak, but the local peak is beneath the line from E to P.

The blue line is correct.

Conclusion

• If P or L is in the square between E and G, pass by them.
• If both P and L is in the square, P has a priority. First, move towards P. Then, consider whether it move towards L or not after arriving at P.
• If the tangent of PL is larger than that of PG, move towards L.
• When L is a candidate to pass by, the tangent of EL is always larger than that of EG.
• The tangent of EP is always smaller than that of PG.

[amaelkady]

@z-ichinohe very good, I agree with the path chosen in all cases.

I understand what you mean about Case 7 but it is unavoidable if you want to stay consistent, I can check this later on with some colleagues; ideally, we want the rule to reflect the realistic reponse of structural members in these cases.

In your conculsions, you mention that the tangent of EL is always larger than that of EG. I'm not sure if that always true; ex. Case 4

Finally, let me know once you finalize the code and feel confident about it, so that I can push it to OpenSees master branch and inform Frank Mckenna and Michael Scott of the updates.

[z-ichinohe]

I had to discard all commits, catch up your latest commit and upload the files again due to my poor branch management. Then, GitHub automatically close the pull request when all commits are discarded.

Thank you for your confirmation. It'll take a while to check the behavior thoroughly.

I forgot to mention that whether the tangent of EL is larger than that of EG or not have importance if only L is a candidate to be the reloading target. Thus, it means that P is not in the square. In such a circumstance, the tangent of EL is always larger than that of EG and we don't need to check that.

[amaelkady]

@z-ichinohe Ok good. I merged your commit.

There are couple of things:

** Can you change the variable name: posResF_0 to posFresFy_0 (and same for negative one) to be consistent with the rest of the name convention (ex, posFmaxFy_0)

** Send me your e-mail, position and affiliation, so that I can communicate these changes to M. Scott and others.

[z-ichinohe]

I changed the names (...)ResF_0 to (...)FresFy_0. Would it be better to reopen a pull request? I'm a graduate student of the department of architecture, The University of Tokyo. e-mail: z-ichinohe@g.ecc.u-tokyo.ac.jp

[amaelkady]

@z-ichinohe Thanks. Yes, commit these changes.