cmp

Author: claudioperez

content/en/examples/Example2/Example2.tcl

@@ -75,3 +75,4 @@ set numIncr 100;	# Number of analysis increments
 source MomentCurvature.tcl
 MomentCurvature 1 $P [expr $Ky*$mu] $numIncr
 
+puts "[nodeDisp 2 3]"

content/en/examples/Example2/MomentCurvature.tcl

@@ -5,66 +5,65 @@
 # 
 # ===----------------------------------------------------------------------===//
 
-# A procedure for performing section analysis (only does
-# moment-curvature, but can be easily modified to do any mode
-# of section response.
+# A procedure for performing section analysis
 #
 # MHS
 # October 2000
 #
 # Arguments
-#	secTag -- tag identifying section to be analyzed
-#	axialLoad -- axial load applied to section (negative is compression)
-#	maxK -- maximum curvature reached during analysis
-#	numIncr -- number of increments used to reach maxK (default 100)
+#    secTag -- tag identifying section to be analyzed
+#    axialLoad -- axial load applied to section (negative is compression)
+#    maxK -- maximum curvature reached during analysis
+#    numIncr -- number of increments used to reach maxK (default 100)
 #
 # Sets up a recorder which writes moment-curvature results to file
 # section$secTag.out ... the moment is in column 1, and curvature in column 2
 
 proc MomentCurvature {secTag axialLoad maxK {numIncr 100} } {
-	# Define two nodes at (0,0)
-	node 1 0.0 0.0
-	node 2 0.0 0.0
+    # Define two nodes at (0,0)
+    node 1 0.0 0.0
+    node 2 0.0 0.0
 
-	# Fix all degrees of freedom except axial and bending
-	fix 1 1 1 1
-	fix 2 0 1 0
+    # Fix all degrees of freedom except axial and bending
+    fix 1 1 1 1
+    fix 2 0 1 0
 
-	# Define element
-	#                         tag ndI ndJ  secTag
-	element zeroLengthSection  1   1   2  $secTag
+    # Define element
+    #                         tag ndI ndJ  secTag
+    element zeroLengthSection  1   1   2  $secTag
 
-	# Create recorder
-	recorder Node -file section$secTag.out -time -node 2 -dof 3 disp
+    # Create recorder
+    recorder Node -file section$secTag.out -time -node 2 -dof 3 disp
 
-	# Define constant axial load
-	pattern Plain 1 "Constant" {
-		load 2 $axialLoad 0.0 0.0
-	}
+    # Define constant axial load
+    pattern Plain 1 "Constant" {
+        load 2 $axialLoad 0.0 0.0
+    }
 
-	# Define analysis parameters
-	integrator LoadControl 0.0
-	system SparseGeneral -piv;	# Overkill, but may need the pivoting!
-	test NormUnbalance 1.0e-9 10
-	numberer Plain
-	constraints Plain
-	algorithm Newton
-	analysis Static
+    # Define analysis parameters
+    integrator LoadControl 0.0
+    system SparseGeneral -piv;    # Overkill, but may need the pivoting!
+    test NormUnbalance 1.0e-9 10
+    numberer Plain
+    constraints Plain
+    algorithm Newton
+    analysis Static
 
-	# Do one analysis for constant axial load
-	analyze 1
+    # Do one analysis for constant axial load
+    analyze 1
 
-	# Define reference moment
-	pattern Plain 2 "Linear" {
-		load 2 0.0 0.0 1.0
-	}
+    # Define reference moment
+    pattern Plain 2 "Linear" {
+        load 2 0.0 0.0 1.0
+    }
 
-	# Compute curvature increment
-	set dK [expr $maxK/$numIncr]
+    # Compute curvature increment
+    set dK [expr $maxK/$numIncr]
 
-	# Use displacement control at node 2 for section analysis
-	integrator DisplacementControl 2 3 $dK 1 $dK $dK
+    # Use displacement control at node 2 for section analysis
+    integrator DisplacementControl 2 3 $dK 1 $dK $dK
 
-	# Do the section analysis
-	analyze $numIncr
+    # Do the section analysis
+    analyze $numIncr
 }
+

content/en/examples/Example2/curvature.py

@@ -0,0 +1,110 @@
+import numpy as np
+
+import xara
+import xara.units.fks as units
+from xara.units.fks import kip, ksi
+
+from xsection.library import from_aisc
+from xsection.analysis import SectionInteraction
+
+import matplotlib.pyplot as plt
+
+
+
+def section_interaction(shape, mat, axial):
+    fig, ax = plt.subplots(1, 1, figsize=(6, 4))
+
+    # test = lambda m: m.eleResponse(1, "section", "fiber", (0.0, 0.0), "stress")[0] < 0.9*50*ksi
+    for N in axial:
+        M, k = zip(*moment_curvature(shape, mat, N, 0.01, 100))
+        ax.plot(k, M, label=f"{N/kip:.0f} kip")
+    plt.show()
+
+
+def moment_curvature(shape, mat, axialLoad, maxK, numIncr):
+    """
+    A procedure for performing section analysis
+
+    Arguments
+       secTag -- tag identifying section to be analyzed
+       axialLoad -- axial load applied to section (negative is compression)
+       maxK -- maximum curvature reached during analysis
+       numIncr -- number of increments used to reach maxK (default 100)
+    """
+
+    model = xara.Model(ndm=3, ndf=6)
+    m = mat
+    model.nDMaterial(m["type"], 1, **{k: v for k, v in m.items() if k != "type"})
+
+    model.section("ShearFiber", 1)
+    
+    for fiber in shape.create_fibers():
+        model.fiber(**fiber, material=1, section=1)
+
+
+    # Define two nodes at (0,0)
+    model.node(1, (0, 0, 0))
+    model.node(2, (0, 0, 0))
+
+    # Fix all degrees of freedom except axial and bending
+    model.fix(1, (1, 1, 1,  1, 1, 1))
+    model.fix(2, (0, 1, 1,  1, 1, 0))
+
+    # Create element
+    model.element("zeroLengthSection", 1, (1, 2), 1)
+
+    # Define constant axial load
+    model.pattern("Plain", 1, "Constant", loads={2: [axialLoad, 0, 0,   0,0,0]})
+
+    # Define analysis
+    model.system("BandGeneral")
+    model.numberer("Plain")
+    model.constraints("Plain")
+    model.test("NormUnbalance", 1.0e-8, 20)
+    model.algorithm("Newton")
+    model.integrator("LoadControl", 0.0)
+    model.analysis("Static")
+
+    # Do one analysis for constant axial load
+    model.analyze(1)
+
+    # Define reference moment
+    model.pattern("Plain", 2, "Linear")
+    model.load(2, (0, 0, 0,  0,0,1), pattern=2)
+
+    # Compute curvature increment
+    dK = maxK/numIncr
+
+    # Use displacement control at node 2 for section analysis
+    model.integrator("DisplacementControl", 2, 6, dK, 1, dK, dK)
+
+    MK = []
+    for i in range(numIncr):
+        # Evaluate step
+        if model.analyze(1) != 0:
+            break
+
+        # Get moment and curvature
+        M = model.eleResponse(1, "force")[5]
+        k = model.nodeDisp(2, 6)
+
+        MK.append([-M, k])
+    return MK
+
+
+if __name__ == "__main__":
+    shape = from_aisc("W14x426", units=units)
+
+    mat = {
+            "type": "J2",
+            "nu": 0.3,
+            "E":  29000*ksi,
+            "Fy":    50*ksi,
+            "Hiso": 0.003*29e3*ksi,
+            "Fs": 0, "Fo": 0, "Hsat": 0
+    }
+
+    axial = np.linspace(5e3*kip, -5e3*kip, 10)
+    SectionInteraction((shape, mat), axial=axial).analyze(nstep=200, incr=1e-4)
+
+    section_interaction(shape, mat, axial)

content/en/examples/FrameHinges/img/curvature.png

@@ -12,3 +12,4 @@ $$
 u = -\frac{M_p L^2}{4 EI}
 $$
 
+![](img/curvature.png)