Tracking plastic forces at member ends

Author: JWock82

PyNite/FEModel3D.py

@@ -2193,7 +2193,7 @@ def _not_ready_yet_analyze_pushover(self, log=False, check_stability=True, push_
 
         # Step through each load combination
         for combo in combo_list:
-
+            
             # Skip the pushover combo
             if combo.name == push_combo:
                 continue
@@ -2202,6 +2202,11 @@ def _not_ready_yet_analyze_pushover(self, log=False, check_stability=True, push_
                 print('')
                 print('- Analyzing load combination ' + combo.name)
 
+            # Reset nonlinear material member end forces to zero
+            for member in self.Members.values():
+                member._fxi, member._myi, member._mzi = 0, 0, 0
+                member._fxj, member._myj, member._mzj = 0, 0, 0
+
             # Get the pushover load step and initialize the load factor
             load_step = list(self.LoadCombos[push_combo].factors.values())[0]
             load_factor = load_step
@@ -2220,13 +2225,14 @@ def _not_ready_yet_analyze_pushover(self, log=False, check_stability=True, push_
             P1_push, P2_push = Analysis._partition(self, self.P(push_combo), D1_indices, D2_indices)
 
             # Solve the current load combination without the pushover load applied
-            Analysis._PDelta_step(self, combo.name, P1, FER1, D1_indices, D2_indices, D2, log, sparse, check_stability, max_iter, tol, first_step=True)
+            Analysis._PDelta_step(self, combo.name, P1, FER1, D1_indices, D2_indices, D2, log, sparse, check_stability, max_iter, first_step=True)
 
-            # Delete this next line used for debugging
-            fx, my, mz = self.Members['M1'].f(combo.name)[0, 0], self.Members['M1'].f(combo.name)[4, 0], self.Members['M1'].f(combo.name)[5, 0]
-            dummy = 1
+            # Since a P-Delta analysis was just run, we'll need to correct the solution to flag it
+            # as 'pushover' instead of 'PDelta'
+            self.solution = 'pushover'
 
-            # Apply the pushover load in steps, summing deformations as we go, until the full pushover load has been analyzed
+            # Apply the pushover load in steps, summing deformations as we go, until the full
+            # pushover load has been analyzed
             while load_factor <= 1:
 
                 # Inform the user which pushover load step we're on
@@ -2236,8 +2242,14 @@ def _not_ready_yet_analyze_pushover(self, log=False, check_stability=True, push_
                 # Run the next pushover load step
                 Analysis._pushover_step(self, combo.name, push_combo, step_num, P1_push, FER1_push, D1_indices, D2_indices, D2, log, sparse, check_stability)
 
-                # Delete this next line used for debugging
-                fx, my, mz = self.Members['M1'].f(combo.name)[0, 0], self.Members['M1'].f(combo.name)[4, 0], self.Members['M1'].f(combo.name)[5, 0]
+                # Update nonlinear material member end forces for each member
+                for member in self.Members.values():
+                    member._fxi = member.f(combo.name, push_combo, step_num)[0, 0]
+                    member._myi = member.f(combo.name, push_combo, step_num)[4, 0]
+                    member._mzi = member.f(combo.name, push_combo, step_num)[5, 0]
+                    member._fxj = member.f(combo.name, push_combo, step_num)[6, 0]
+                    member._myj = member.f(combo.name, push_combo, step_num)[10, 0]
+                    member._mzj = member.f(combo.name, push_combo, step_num)[11, 0]
 
                 # Move on to the next load step
                 load_factor += load_step

PyNite/Member3D.py

@@ -49,6 +49,14 @@ def __init__(self, name, i_node, j_node, material, model, Iy, Iz, J, A, auxNode=
             self.Iz = model.Sections[section_name].Iz
             self.J = model.Sections[section_name].J
 
+        # Variables used to track nonlinear material member end forces
+        self._fxi = 0
+        self._myi = 0
+        self._mzi= 0
+        self._fxj = 0
+        self._myj = 0
+        self._mzj = 0
+
         self.auxNode = auxNode # Optional auxiliary node used to define the member's local z-axis
         self.PtLoads = []   # A list of point loads & moments applied to the element (Direction, P, x, case='Case 1') or (Direction, M, x, case='Case 1')
         self.DistLoads = [] # A list of linear distributed loads applied to the element (Direction, w1, w2, x1, x2, case='Case 1')
@@ -59,14 +67,10 @@ def __init__(self, name, i_node, j_node, material, model, Iy, Iz, J, A, auxNode=
         self.tension_only = tension_only # Indicates whether the member is tension-only
         self.comp_only = comp_only # Indicates whether the member is compression-only
 
-        # Members need to track whether they are active or not for any given load combination.
-        # They may become inactive for a load combination during a tension/compression-only
-        # analysis. This dictionary will be used when the model is solved.
+        # Members need to track whether they are active or not for any given load combination. They may become inactive for a load combination during a tension/compression-only analysis. This dictionary will be used when the model is solved.
         self.active = {} # Key = load combo name, Value = True or False
 
-        # The 'Member3D' object will store results for one load combination at a time. To reduce repetative calculations
-        # the '__solved_combo' variable will be used to track whether the member needs to be resegmented before running
-        # calculations for any given load combination.
+        # The 'Member3D' object will store results for one load combination at a time. To reduce repetative calculations the '__solved_combo' variable will be used to track whether the member needs to be resegmented before running calculations for any given load combination.
         self.__solved_combo = None # The current solved load combination
 
         # Members need a link to the model they belong to
@@ -230,9 +234,6 @@ def km(self, combo_name='Combo 1', push_combo='Push', step_num=1):
         :return: The plastic reduction matrix for the element
         :rtype: array
         """
-
-        # Get the total elastic end forces applied to the element
-        f = self.f(combo_name) - self.fer(combo_name) - self.fer(push_combo)*step_num
 
         # Get the elastic local stiffness matrix
         ke = self.k()
@@ -251,9 +252,8 @@ def km(self, combo_name='Combo 1', push_combo='Push', step_num=1):
         if self.section is None:
             raise Exception('Nonlinear material analysis requires member sections to be defined. A section definition is missing for element ' + self.name + '.')
         else:
-            # TODO: Note that we have assumed that `f` is based on the total elastic load acting on the member at this load step, rather than just the change in load for this load step. This seems appropriate, since G is the gradient to the yield surface, and where we are heading along that surface depends on the total load applied, rather than a part of the total load. Need to verify this works correctly with test cases.
-            Gi = self.section.G(f[0, 0], f[4, 0], f[5, 0])
-            Gj = self.section.G(f[6, 0], f[10, 0], f[11, 0])
+            Gi = self.section.G(self._fxi, self._myi, self._mzi)
+            Gj = self.section.G(self._fxj, self._myj, self._mzj)
 
         # Expand the gradients for a 12 degree of freedom element
         zeros_array = zeros((6, 1))
@@ -427,8 +427,8 @@ def _partition(self, unp_matrix):
             return  m11, m12, m21, m22
 
 #%%   
-    def f(self, combo_name='Combo 1'):
-        """Returns the member's local end force vector for the given load combination.
+    def f(self, combo_name='Combo 1', push_combo='Push', step_num=1):
+        """Returns the member's elastic local end force vector for the given load combination.
 
         :param combo_name: The load combination to get the local end for vector for. Defaults to 'Combo 1'.
         :type combo_name: str, optional
@@ -441,6 +441,9 @@ def f(self, combo_name='Combo 1'):
             # Back-calculate the axial force on the member from the axial strain
             P = (self.d(combo_name)[6, 0] - self.d(combo_name)[0, 0])*self.A*self.E/self.L()
             return add(matmul(add(self.k(), self.kg(P)), self.d(combo_name)), self.fer(combo_name))
+        elif self.model.solution == 'Pushover':
+            P = self._fxj  -self._fxi
+            return add(matmul(add(self.k(), self.kg(P), self.km(combo_name, push_combo, step_num)), self.d(combo_name)), self.fer(combo_name))
         else:
             return add(matmul(self.k(), self.d(combo_name)), self.fer(combo_name))