cmp

Author: claudioperez

assets/img/examples/e0020.png

@@ -69,7 +69,6 @@ is created using steel and concrete fibers.
 
 ### Materials
 
-
 ![](img/RCsection0.svg)
 
 {{< tabs tabTotal="2" >}}

content/en/examples/Example3/index.md

@@ -1,17 +1,25 @@
 ---
-title: "Linearized Buckling"
-tags: ["Python", "Tcl", "Frame"]
+#title: "Linearized Buckling"
+title: "Column Buckling"
+meta:
+  title: "A nonlinear buckling simulation of a column with varied boundary conditions and finite element formulations"
+tags: ["Frame"]
 #render: ./model.glb
-description: Corotational frame elements are used to approximate Euler's buckling load.
+description: >-
+  This example performs a nonlinear buckling analysis of a column under axial load,
+  using various element formulations and end conditions.
 downloads:
   Python: ["buckling.py"]
   Tcl:    ["buckling.tcl"]
+weight: 20
 ---
 
 
 ## Problem
 
-Corotational frame elements are used to approximate Euler's buckling load which is given by:
+This example analyzes a column subject to axial loading until buckling occurs.
+Different element types, shear deformation settings, and boundary conditions are compared
+by computing the critical load and comparing it to the theoretical Euler load which is given by:
 
 \[
 P_{\mathrm{euler}} = \frac{\pi^2 EI}{L^2}
@@ -57,3 +65,276 @@ $$
 $$
 \tan \lambda_{\text {cr }}=\chi \lambda_{\text {cr }}=\frac{\lambda_{\text {cr }}}{1+\lambda_{\text {cr }} 2 \varphi / 12}
 $$
+
+## Implementation
+
+
+## Model
+
+We begin by importing the required modules and defining some material and section properties.
+
+{{< tabs tabTotal="2" >}}
+{{% tab name="Python" %}}
+```python
+import opensees.openseespy as ops
+import numpy as np
+from math import pi
+
+FACTORS = {
+    "pin-pin":     1,
+    "fix-roll":    1,
+    "fix-fix":     0.5,
+    "fix-pin":     0.7,
+    "fix-free":    2,
+    "pin-roll":    2,
+}
+```
+{{% /tab %}}
+{{< /tabs >}}
+
+We define a function `create_column` that constructs the finite element model of the column
+based on options like element type, shear deformation flag, and spatial dimension:
+
+{{< tabs tabTotal="2" >}}
+{{% tab name="Python" %}}
+```python
+def create_column(boundary="pin-pin", elem_data=None, ndm=2):
+    if elem_data is None:
+        elem_data = {}
+
+    elem_type  = elem_data.get("type",      "forceBeamColumnCBDI")
+    geom_type  = elem_data.get("transform", "Corotational")
+    use_shear  = elem_data.get("shear",     False)
+    kinematics = elem_data.get("order",     0)
+
+    E  = 29000.0
+    G  = 11200.0
+    A  = 112.0
+    I  = 110.0
+    Ay = 3/6*A
+    Az = 3/6*A
+    L  = 60.0
+    ne = 10
+    nIP = 5
+
+    model = ops.Model(ndm=ndm)
+
+    for i in range(1, ne+2):
+        y = (i-1)/float(ne)*L
+        location = (0.0, y, 0.0) if ndm == 3 else (0.0, y)
+        model.node(i, location)
+        model.mass(i, *[1.0]*model.getNDF())
+
+    fix_node(model, 1, boundary.split("-")[0])
+    fix_node(model, ne+1, boundary.split("-")[1])
+```
+{{% /tab %}}
+{{< /tabs >}}
+
+The section is defined as an `ElasticFrame` and assigned to each element in the column:
+
+{{< tabs tabTotal="2" >}}
+{{% tab name="Python" %}}
+```python
+    sec_tag = 1
+    properties = [E, A, I]
+    if ndm == 3:
+        properties += [2*I, "-J", 100*I, "-G", G]
+    if use_shear:
+        properties += ["-Ay", Ay, "-G", G]
+        if ndm == 3:
+            properties += ["-Az", Az]
+    model.section("ElasticFrame", sec_tag, *properties)
+
+    geo_tag = 1
+    vector = (0, 0, 1) if ndm == 3 else ()
+    model.geomTransf(geom_type, geo_tag, *vector)
+
+    for i in range(1, ne+1):
+        model.element(elem_type, i, (i, i+1),
+                      section=sec_tag,
+                      transform=geo_tag,
+                      order=kinematics,
+                      shear=int(use_shear))
+```
+{{% /tab %}}
+{{< /tabs >}}
+
+The axial load is applied at the top node, using the Euler buckling expression as a reference:
+
+{{< tabs tabTotal="2" >}}
+{{% tab name="Python" %}}
+```python
+    if use_shear:
+        phi = 12*E*I/(Ay*G*L**2)
+        lam = buckle_factor(boundary, phi)
+        kL = L/lam
+        euler_load = E*I/kL**2 / (1 + lam**2*phi/12)
+    else:
+        kL = L / buckle_factor(boundary)
+        euler_load = E*I/kL**2
+
+    model.pattern("Plain", 1, "Linear")
+    load = (0.0, -euler_load, 0.0) if ndm == 2 else (0.0, -euler_load, 0.0, 0, 0, 0)
+    model.load(ne+1, *load, pattern=1)
+
+    return model, euler_load
+```
+{{% /tab %}}
+{{< /tabs >}}
+
+## Buckling Analysis
+
+A nonlinear static analysis is performed to determine the limit load. The first
+eigenvalue of the tangent stiffness matrix is tracked and interpolated to find the
+load factor at which it reaches zero.
+
+{{< tabs tabTotal="2" >}}
+{{% tab name="Python" %}}
+```python
+def buckling_analysis(model, peak_load):
+    load_step = 0.01
+    PeakLoadRatio = 1.5
+
+    model.system('UmfPack')
+    model.test("EnergyIncr", 1e-8, 20, 9)
+    model.algorithm("Newton")
+    model.integrator("LoadControl", load_step)
+    model.analysis("Static")
+
+    lam_0 = model.getTime()
+    eig_0 = model.eigen("-fullGenLapack", 1)[0]
+
+    limit_load = None
+    for i in range(1, int(PeakLoadRatio/load_step)+1):
+        if model.analyze(1) != 0:
+            break
+        lam = model.getTime()
+        eig = model.eigen("-fullGenLapack", 1)[0]
+        if eig <= 0.0:
+            lam_i = lam_0 + (lam - lam_0)*eig_0/(eig_0-eig)
+            limit_load = lam_i * peak_load
+            break
+        lam_0 = lam
+        eig_0 = eig
+
+    return limit_load
+```
+{{% /tab %}}
+{{< /tabs >}}
+
+## Results
+
+The following script compares element formulations and boundary conditions in 3D with and without shear deformation:
+
+{{< tabs tabTotal="2" >}}
+{{% tab name="Python" %}}
+```python
+if __name__ == "__main__":
+    for ndm in (3,):
+        for shear in (False, True):
+            for elem in ("PrismFrame", "ForceFrame", "ExactFrame"):
+                print(f"\n{elem:10}      Shear    Order     Theory   Computed       Error")
+                for boundary in FACTORS:
+                    orders = (0,)
+                    for order in orders:
+                        if not shear and "Exact" in elem:
+                            continue
+                        elem_data = {
+                            "type": elem,
+                            "shear": shear,
+                            "order": order,
+                            "transform": "Corotational"
+                        }
+                        model, euler_load = create_column(boundary, elem_data, ndm=ndm)
+                        limit_load = buckling_analysis(model, euler_load)
+                        model.wipe()
+                        del model
+                        print(f"  {boundary:10} {shear:8} {order:8} ", end="")
+                        if limit_load is None:
+                            print("No singularity found.")
+                        else:
+                            err = 100*(limit_load/euler_load - 1)
+                            print(f"{euler_load:10.2f} {limit_load:10.2f} {err:10.3f} %")
+```
+{{% /tab %}}
+{{< /tabs >}}
+
+
+### Euler
+
+| Boundary    | Critical Load |
+|-------------|---------------|
+| pin-pin     |   8745.57 |
+| fix-roll    |   8745.57 |
+| fix-fix     |  34982.26 |
+| fix-pin     |  17891.23 |
+| fix-free    |   2186.39 |
+| pin-roll    |   2186.39 |
+
+PrismFrame
+
+|  Boundary   |   Computed    |   Error   |
+|-------------|---------------|-----------|
+|  pin-pin    |    8841.80    |  1.100 %  |
+|  fix-roll   |    8841.80    |  1.100 %  |
+|  fix-fix    |   36559.09    |  4.507 %  |
+|  fix-pin    |   18296.87    |  2.267 %  |
+|  fix-free   |    2192.37    |  0.273 %  |
+|  pin-roll   |    2192.37    |  0.273 %  |
+
+ForceFrame 
+
+|  Boundary   |  Computed  |   Error  |
+|-------------|------------|----------|
+|  pin-pin    |   8841.80  |  1.100 % |
+|  fix-roll   |   8841.80  |  1.100 % |
+|  fix-fix    |  36559.09  |  4.507 % |
+|  fix-pin    |  18296.87  |  2.267 % |
+|  fix-free   |   2192.37  |  0.273 % |
+|  pin-roll   |   2192.37  |  0.273 % |
+
+### Shear
+
+PrismFrame
+
+| Boundary  |  Computed  |   Error   |
+|-----------|------------|-----------|
+| pin-pin   |   8718.89  |  1.085 %  |
+| fix-roll  |   8718.88  |  1.085 %  |
+| fix-fix   |  34545.22  |  4.259 %  |
+| fix-pin   |  17727.65  |  2.192 %  |
+| fix-free  |   2184.73  |  0.272 %  |
+| pin-roll  |   2184.73  |  0.273 %  |
+
+ForceFrame
+
+| Boundary  |  Computed  |   Error   |
+|-----------|------------|-----------|
+| pin-pin   |   8718.89  |  1.085 %  |
+| fix-roll  |   8718.88  |  1.085 %  |
+| fix-fix   |  34545.22  |  4.259 %  |
+| fix-pin   |  17727.65  |  2.192 %  |
+| fix-free  |   2184.73  |  0.272 %  |
+| pin-roll  |   2184.73  |  0.273 %  |
+
+ExactFrame
+
+| Boundary  |  Computed  |   Error   |
+|-----------|------------|-----------|
+| pin-pin   |   8792.02  |  1.933 %  |
+| fix-roll  |   8792.01  |  1.933 %  |
+| fix-fix   |  35773.31  |  7.965 %  |
+| fix-pin   |  18068.91  |  4.159 %  |
+| fix-free  |   2189.24  |  0.480 %  |
+| pin-roll  |   2189.25  |  0.480 %  |
+
+
+| Boundary  | Critical Load |
+|-----------|------------|
+| pin-pin   |   8625.30
+| fix-roll  |   8625.30
+| fix-fix   |  33134.20
+| fix-pin   |  17347.40
+| fix-free  |   2178.80
+| pin-roll  |   2178.80