Balken-Element B2

B: Beam, 2: Zwei Knoten

Lineares System

../../../_images/beam-fem1.png

Bemerkung

Herleitung dieser Gleichung: Herleitung 1 (Klassisch) und Herleitung 2 (PdvV).

Schreibweisen

Verschiedene Schreibweisen

Es gibt verschiedene Schreibweisen. Denn man kann wählen:

  • Welche Symbole man verwendet: \(\psi_1, w_1, \ldots , M_2, F_2\) und \(E, I, l\) können auch anders bezeichnet werden.

  • Welche Zählrichtungen man für dievektorwertigen Größen verwendet. D.h. in welche Richtung man die zugehörigen Pfeile zeichnet.

  • Wie man die 4 Aussagen in einer Matrix anordnet.

../../../_images/reorder.png

Oben: Hier verwendete Schreibweise. Unten: Schreibweise ähnlich wie im Lehrbuch FEM von Dr. Klein.

Erweiterte Steifigkeitsmatrix

Die erweiterte Element-Steifigkeitsmatrix bietet den Vorteil, dass die Freiheitsgrade \((\psi_1, w_1, \psi_2, w_2)\) sichtbar sind, auf die sich die Einträge der Steifigkeitsmatrix beziehen.

\begin{align*} k^| = \tfrac{EI}{l^3} \left[ \begin{array}{cccc|c} {\color{green}4} l^2 & {\color{red}-}{\color{blue}6} l & {\color{magenta}2} l^2 & {\color{magenta}6} l & \psi_1 \\ {\color{red}-}6 l & {\color{green}12} & {\color{red}-}{\color{blue}6} l & {\color{red}-}{\color{magenta}12} & w_1 \\ 2 l^2 & {\color{red}-}6 l & {\color{green}4} l^2 & {\color{blue}6} l & \psi_2 \\ 6 l & {\color{red}-}12 & 6 l & {\color{green}12} & w_2 \\ \end{array} \right] \end{align*}

Merkregel mit 5 Schritten

  • Diagonale: \(\color{green}4\,12\,4\,12\) = „Ver-Zweiflung Ver-Zweiflung“.

  • Ein Kreuz aus sechs Minuszeichen \(\color{red}-\) = „Hoffnung auf „Gott“.

  • \(\color{blue}666\) = „der Antichrist“.

  • \({\color{magenta}2} \cdot {\color{magenta}6} = {\color{magenta}12}\) aus dem „Kleinen Einmaleins“.

  • Einheitenbetrachtung liefert \(l\) und \(l^2\) an den richtigen Stellen.

Wahl der Knoten-Reihenfolge

Wie beim Stab-Element lässt sich die Knoten-Reihenfolge vertauschen, ohne dass sich die Einträge der Steifigkeitsmatrix ändern:

\begin{align*} \underbrace{ \tfrac{EI}{l^3} \begin{bmatrix} 4 l^2 & -6 l & 2 l^2 & 6 l \\ -6 l & 12 & -6 l & -12 \\ 2 l^2 & -6 l & 4 l^2 & 6 l \\ 6 l & -12 & 6 l & 12 \\ \end{bmatrix} }_k \begin{bmatrix} \psi_1 \\ w_1 \\ \psi_2 \\ w_2 \\ \end{bmatrix} &= \begin{bmatrix} M_1 \\ F_1 \\ M_2 \\ F_2 \\ \end{bmatrix} + q \begin{bmatrix} - \tfrac{1}{12}l^2 \\ \tfrac12 l\\ \tfrac{1}{12}l^2\\ \tfrac12 l \\ \end{bmatrix} \\ \underbrace{ \tfrac{EI}{l^3} \begin{bmatrix} 4 l^2 & -6 l & 2 l^2 & 6 l \\ -6 l & 12 & -6 l & -12 \\ 2 l^2 & -6 l & 4 l^2 & 6 l \\ 6 l & -12 & 6 l & 12 \\ \end{bmatrix} }_k \begin{bmatrix} \psi_2 \\ w_2 \\ \psi_1 \\ w_1 \\ \end{bmatrix} &= \begin{bmatrix} M_2 \\ F_2 \\ M_1 \\ F_1 \\ \end{bmatrix} + q \begin{bmatrix} \tfrac{1}{12}l^2\\ \tfrac12 l \\ - \tfrac{1}{12}l^2 \\ \tfrac12 l\\ \end{bmatrix} \end{align*}

Zugehörige erweiterte Steifigkeitsmatrizen:

\begin{align*} k^| &= \tfrac{EI}{l^3} \left[ \begin{array}{cccc|c} 4 l^2 & -6 l & 2 l^2 & 6 l & \psi_1 \\ -6 l & 12 & -6 l & -12 & w_1 \\ 2 l^2 & -6 l & 4 l^2 & 6 l & \psi_2 \\ 6 l & -12 & 6 l & 12 & w_2 \\ \end{array} \right] \\ k^| &= \tfrac{EI}{l^3} \left[ \begin{array}{cccc|c} 4 l^2 & -6 l & 2 l^2 & 6 l & \psi_2 \\ -6 l & 12 & -6 l & -12 & w_2 \\ 2 l^2 & -6 l & 4 l^2 & 6 l & \psi_1 \\ 6 l & -12 & 6 l & 12 & w_1 \\ \end{array} \right] \end{align*}

Postprocessing / Interpolation

\(N_1, N_2, N_3, N_4, N_5\) sind die dimensionslosen Interpolationsfunktionen:

\[\begin{split}N_1 &= - \xi^{3} + 2 \xi^{2} - \xi\\ N_2 &= 2 \xi^{3} - 3 \xi^{2} + 1 \\ N_3 &= - \xi^{3} + \xi^{2} \\ N_4 &= - 2 \xi^{3} + 3 \xi^{2} \\ N_5 &= \tfrac{\xi^{4}}{24} - \tfrac{\xi^{3}}{12} + \tfrac{\xi^{2}}{24}\\\end{split}\]

Querverschiebung w

Interpolierte Querverschiebung in Matrix-Schreibweise:

\[\begin{split}w &= w_{\psi_1} + w_{w_1} + w_{\psi_2} + w_{w_2} + w_q \\ &= \begin{bmatrix} \psi_1 & w_1 & \psi_2 & w_2 & q \end{bmatrix} \begin{bmatrix} l N_1 \\ N_2 \\ l N_3 \\ N_4 \\ \tfrac{l^4}{B} N_5 \end{bmatrix} \\ &= \begin{bmatrix} \psi_1 & w_1 & \psi_2 & w_2 & q \end{bmatrix} \begin{bmatrix} l \left( - \xi^{3} + 2 \xi^{2} - \xi \right) \\ 2 \xi^{3} - 3 \xi^{2} + 1 \\ l \left( - \xi^{3} + \xi^{2} \right) \\ - 2 \xi^{3} + 3 \xi^{2} \\ \tfrac{l^4}{B} \left( \tfrac{\xi^{4}}{24} - \tfrac{\xi^{3}}{12} + \tfrac{\xi^{2}}{24} \right) \end{bmatrix}\end{split}\]

Bemerkung

  • Eine \(1 \times 5\)-Matrix multipliziert mit einer \(5 \times 1\)-Matrix ergibt eine \(1 \times 1\)-Matrix, also eine Zahl.

  • Man kann das Produkt auch so schreiben:

    \[\begin{split}w &= \underbrace{ \begin{bmatrix} \psi_1 & w_1 & \psi_2 & w_2 & q \end{bmatrix} }_{1 \times 5} \underbrace{ \begin{bmatrix} l N_1 \\ N_2 \\ l N_3 \\ N_4 \\ \tfrac{l^4}{B} N_5 \end{bmatrix} }_{5 \times 1} \\ &= \underbrace{ \begin{bmatrix} l N_1 & N_2 & l N_3 & N_4 & \tfrac{l^4}{B} N_5 \end{bmatrix} }_{1 \times 5} \underbrace{ \begin{bmatrix} \psi_1 \\ w_1 \\ \psi_2 \\ w_2 \\ q \end{bmatrix} }_{5 \times 1}\end{split}\]
  • \(w\) ist von der Ordnung \(\xi^4\).

Verdrehung ψ, Biegemoment M, Querkraft Q

Mit \(\psi=-w'\) und \(M = - B w''\) und \(Q = - B w'''\) gilt:

\[\begin{split}\psi &= - \begin{bmatrix} \psi_1 & w_1 & \psi_2 & w_2 & q \end{bmatrix} \begin{bmatrix} l N_1' \\ N_2' \\ l N_3' \\ N_4' \\ \tfrac{l^4}{B} N_5' \end{bmatrix} \\ &= \begin{bmatrix} \psi_1 & w_1 & \psi_2 & w_2 & q \end{bmatrix} \begin{bmatrix} 3 \xi^{2} - 4 \xi + 1\\- \tfrac{1}{l} \left(6 \xi^{2} - 6 \xi\right)\\ 3 \xi^{2} - 2 \xi\\ - \tfrac{1}{l} \left(- 6 \xi^{2} + 6 \xi\right)\\ - \tfrac{l^{3}}{B} \left(\tfrac{\xi^{3}}{6} - \tfrac{\xi^{2}}{4} + \tfrac{\xi}{12}\right) \end{bmatrix}\end{split}\]

\[\begin{split}M&= - B \begin{bmatrix} \psi_1 & w_1 & \psi_2 & w_2 & q \end{bmatrix} \begin{bmatrix} l N_1'' \\ N_2'' \\ l N_3'' \\ N_4'' \\ \tfrac{l^4}{B} N_5'' \end{bmatrix} \\ &= \begin{bmatrix} \psi_1 & w_1 & \psi_2 & w_2 & q \end{bmatrix} \begin{bmatrix} - \tfrac{B}{l} \left(- 6 \xi + 4\right) \\ - \tfrac{B}{l^{2}} \left(12 \xi - 6\right) \\ - \tfrac{B}{l} \left(- 6 \xi + 2\right) \\ - \tfrac{B}{l^{2}} \left(- 12 \xi + 6\right) \\ - l^{2} \left(\tfrac{\xi^{2}}{2} - \tfrac{\xi}{2} + \tfrac{1}{12}\right) \end{bmatrix}\end{split}\]

\[\begin{split}Q&= - B \begin{bmatrix} \psi_1 & w_1 & \psi_2 & w_2 & q \end{bmatrix} \begin{bmatrix} l N_1''' \\ N_2''' \\ l N_3''' \\ N_4''' \\ \tfrac{l^4}{B} N_5''' \end{bmatrix} \\ &= \begin{bmatrix} \psi_1 & w_1 & \psi_2 & w_2 & q \end{bmatrix} \begin{bmatrix} \tfrac{6 B}{l^{2}}\\ - \tfrac{12 B}{l^{3}}\\ \tfrac{6 B}{l^{2}}\\ \tfrac{12 B}{l^{3}}\\ - l \left(\xi - \tfrac{1}{2}\right) \end{bmatrix}\end{split}\]

Graphische Darstellung

Hinweise zu den Diagrammen

Zwei Interaktive Diagramme

  • Schieberegler verschieben = Variieren des Knoten-Werts.

  • In Legende Funktion klicken = Aus- und Einblenden des Funktionsgraphen.

  • Horizontale Achse: Werte werden nach rechts größer - wie üblich.

  • Vertikale Achse: Werte werden nach unten größer.

Querverschiebung, Verdrehung, Biegemoment, Querkraft

Displacement, Rotation, Bending Mom., Shear Force (B = 1 Nm², l = 1 m)

Interpolationsfunktionen

Shape Functions

Äquivalente Knotenlasten

Video

../../../_images/arbitrary-load.png

Knotenlasten aufgrund äußerer Lasten \(q_\xi, M_\xi, F_\xi\) zwischen den Element-Knoten. Oben: Lineares System. Mitte: Generische Lasten. Unten: Äquivalente Knotenlasten,

Details

Generische Lasten -> Generische Rechte Seite

  • Momente \(M_1, M_2\) sowie Kräfte \(F_1, F_2\) an den Knoten.

  • Pro Element konstante verteilte \(q\)-Last.

  • Die FEM-Lösung gleicht der klassischen Lösung, denn das Balken-Element wurde definiert, um diese generischen Lasten abzubilden.

  • Die rechte Seite \(f\) im Gleichungssystem \(ku = f\) ist:

\[\begin{split}f= \begin{bmatrix} M_1 \\ F_1 \\ M_2 \\ F_2 \\ \end{bmatrix} + q \begin{bmatrix} - \tfrac{1}{12}l^2 \\ \tfrac12 l\\ \tfrac{1}{12}l^2\\ \tfrac12 l \\ \end{bmatrix}\end{split}\]

Sonstige Lasten -> Rechte Seite mit Äquivalenten Knotenlasten

  • Momente \(M_\xi\) an beliebiger Position \(\xi_M\) sowie Kräfte \(F_\xi\) an beliebiger Position \(\xi_F\).

  • Beliebig verteilte \(q_\xi\)-Last.

  • Die FEM-Lösung gleicht nicht der klassischen Lösung.

  • Die sonstigen Lasten werden umgerechnet in die Äquivalenten Knotenlasten. 1 Diese Äquivalenten Knotenlasten bilden die rechte Seite \(f\):

\[\begin{split}f = M_\xi \begin{bmatrix} - l N_1' \\ - N_2' \\ - l N_3' \\ - N_4' \\ \end{bmatrix}_{\xi=\xi_M} + F_\xi \begin{bmatrix} l N_1 \\ N_2 \\ l N_3 \\ N_4 \\ \end{bmatrix}_{\xi=\xi_F} + \int_{\xi=0}^1 q_\xi \begin{bmatrix} l N_1 \\ N_2 \\ l N_3 \\ N_4 \\ \end{bmatrix} l \mathsf{d}\xi\end{split}\]

mit Interpolationsfunktionen und deren Ableitungen.

Fußnote:

1

Die bei der Umrechnung verwendete Annahme ist, dass die virtuelle äußere Arbeit der sonstigen Lasten gleich ist der virtuellen äußeren Arbeit der daraus berechneten Äquivalenten Knotenlasten.

Beispiele

../../../_images/ex-1.png

Spaltenmatrizen = Äquivalente Knotenlasten: Berechnung siehe Python-Programm.

Lösung mit Python: Copy - Paste - Play

  • Copy: Source Code (siehe unten) aufklappen und kopieren.

  • Paste: Einfügen als Python-Notebook auf:

  • Play: Ausführen.

Source Code

from sympy import *

half = S(1)/2

M, F, q, l, xi = var("M, F, q, l, xi")

N1 = -xi**3 + 2*xi**2 - xi
N2 = 2*xi**3 - 3*xi**2 + 1
N3 = -xi**3 + xi**2
N4 = -2*xi**3 + 3*xi**2

N = Matrix([l*N1, N2, l*N3, N4])
# N':
Nx = diff(N, xi) / l

pprint("\nConcentrated Moment Mξ at ξ = 1/2:")
f  = - M * Nx.subs(xi, half)
pprint("\nf / M:")
tmp = f/M
pprint(tmp)

pprint("\nConcentrated Force Fξ at ξ = 1/2:")
f  = F * N.subs(xi, half)
pprint("\nf / F:")
tmp = f/F
pprint(tmp)

pprint("\nPiecewise linear distributed load qξ:")
# Piecewise defined function:
q1, q2 = 2*q*xi, 2*q*(1 - xi)
# Piecewise integration:
f  = integrate(q1*N*l, (xi, 0, half))
f += integrate(q2*N*l, (xi, half, 1))
pprint("\nf / (-ql):")
tmp = f/(-q*l)
pprint(tmp)

pprint("\nSpecial Cases:")
pprint("\nConstant distributed load qξ=q:")
f  = integrate(q*N*l, (xi, 0, 1))
pprint("\nf / (-ql):")
tmp = f/(-q*l)
pprint(tmp)

pprint("\nM1 and F1 at ξ = 0 and M2 and F2 at ξ = 1:")
M1, F1, M2, F2 = var("M1, F1, M2, F2")
f  = - M1 * Nx.subs(xi, 0)
f +=   F1 *  N.subs(xi, 0)
f += - M2 * Nx.subs(xi, 1)
f +=   F2 *  N.subs(xi, 1)
pprint("\nf:")
tmp = f
pprint(tmp)

Aufgaben-Katalog

../../../_images/ex.png

Links: Aufgaben 5.7 (Dr. Klein), B2.A , B2.B, 5.8 (Dr. Klein) = B2.C, B2.D. Rechts: Aufgaben B2.E, B2.F, B2.G, B2.H. Lösungen siehe Python-Programm.

Lösung mit Python: Copy - Paste - Play

  • Copy: Source Code (siehe unten) aufklappen und kopieren.

  • Paste: Einfügen als Python-Notebook auf:

  • Play: Ausführen.

Source Code

# -*- coding: utf-8 -*-
from sympy.physics.units import *
from sympy import *


# Rounding:
import decimal
from decimal import Decimal as DX
from copy import deepcopy
def iso_round(obj, pv, rounding=decimal.ROUND_HALF_EVEN):
    import sympy
    """
    Rounding acc. to DIN EN ISO 80000-1:2013-08
    place value = Rundestellenwert
    """
    assert pv in set([
        # place value     #  round to:
        "1",              #  closest (even) integer
        "0.1",            #  1st digit after decimal
        "0.01",           #  2nd
        "0.001",          #  3rd
        "0.0001",         #  4th
        "0.00001",        #  5th
        "0.000001",       #  6th
        "0.0000001",      #  7th
        "0.00000001",     #  8th
        "0.000000001",    #  9th
        "0.0000000001",   # 10th
        ])
    objc = deepcopy(obj)
    try:
        tmp = DX(str(float(objc)))
        objc = tmp.quantize(DX(pv), rounding=rounding)
    except:
        for i in range(len(objc)):
            tmp = DX(str(float(objc[i])))
            objc[i] = tmp.quantize(DX(pv), rounding=rounding)
    return objc

# Units:
(mm, cm)  =  ( m/1000, m/100 )
Newton    =  kg*m/s**2
kN        =  10**3*Newton
Pa        =  Newton/m**2
MPa       =  10**6*Pa
GPa       =  10**9*Pa
deg       =  pi/180
half      =  S(1)/2

a, c, l, q, F, M, EI = var("a, c, l, q, F, M, EI")

quants = False
# quants = True
# sub_list=[
#     ( M,   10 *Newton*m              ),
#     ( a,    1 *m                     ),
#     ( EI, 200*GPa * 2*mm*6*mm**3/ 12 ),
#     ]

def K_1E(EI, l):
    # l = Element length
    l2, l3 = l*l, l*l*l
    K = EI/l3
    K *= Matrix([
        [  4*l2 ,  -6*l ,  2*l2 ,   6*l  ],
        [ -6*l  ,  12   , -6*l  , -12    ],
        [  2*l2 ,  -6*l ,  4*l2 ,   6*l  ],
        [  6*l  , -12   ,  6*l    , 12   ],
        ])
    return K

def K_2E(EI, l):
    # l = Element length
    l2, l3 = l*l, l*l*l
    K = EI/l3
    K *= Matrix(
    [
    [  4*l2 ,  -6*l ,  2*l2 ,   6*l ,  0    ,   0    ],
    [ -6*l  ,  12   , -6*l  , -12   ,  0    ,   0    ],
    [  2*l2 ,  -6*l ,  8*l2 ,   0   ,  2*l2 ,   6*l  ],
    [  6*l  , -12   ,  0    ,  24   , -6*l  , -12    ],
    [  0    ,   0   ,  2*l2 ,  -6*l ,  4*l2 ,   6*l  ],
    [  0    ,   0   ,  6*l  , -12   ,  6*l  ,  12    ],
    ]
    )
    return K

def Fq(l):
    return q*l/2
def Mq(l):
    return q*l*l / 12

p0,w0,p1,w1,p2,w2 = var("ψ₀,w₀,ψ₁,w₁,ψ₂,w₂")
M0,F0,M1,F1,M2,F2 = var("M₀,F₀,M₁,F₁,M₂,F₂")

case = "5.7"
# case = "B2.A"
# case = "B2.B"
# case = "B2.C"
# case = "B2.D"
# case = "B2.E"
# case = "B2.F"
# case = "B2.G"
# case = "B2.H"
case = "AB-1"
case = "AB-2"

pprint("\nCase: "+case)

if case == "5.7":
    K = K_1E(EI, a)
    #                    1                             2
    #
    #                                                  M
    #                    |------------------------------
    #                                    a             Λ
    M = var("M")
    #
    u = Matrix([       0,0,                          p2,0      ])
    f = Matrix([      M1,F1,                          M,F2     ])
    unks =     [      M1,F1,                         p2,F2     ]

elif case == "B2.A":
    K = K_2E(EI, a)
    #                    0              1              2
    #                                   F1             F2
    #                    |------------------------------
    u = Matrix([        0,0,          p1,w1,         p2,w2       ])
    f = Matrix([       M0,F0,          0,F1,          0,F2       ])
    unks =     [       M0,F0,         p1,w1,         p2,w2       ]

elif case == "B2.B":
    K = K_1E(EI, a)
    #                    1                             2
    #                     qqqqqqqqqqqqqqqqqqqqqqqqqqqqq
    #                    |------------------------------
    #                                   a              Λ
    #
    u = Matrix([        0,0,                        p2,0         ])
    f = Matrix([ M1-Mq(a),F1+Fq(a),            0+Mq(a),F2+Fq(a)  ])
    unks =     [       M1,F1,                       p2,F2        ]

elif case == "B2.C":
    K = K_2E(EI, a/4)
    #                    0              1              2
    #                    ◆ qqqqqqqqqqqqqqqqqqqqqqqqqqqq
    #                    -------------------------------
    #                    ◆                             Λ
    #
    u = Matrix([        0,w0,        p1,w1,        p2,0             ])
    f = Matrix([M0-Mq(a/4),Fq(a/4), 0,2*Fq(a/4), Mq(a/4),F2+Fq(a/4) ])
    unks =     [       M0,w0,        p1,w1,        p2,F2            ]

elif case == "B2.D":
    K = K_1E(EI, a)
    #                    1                             2
    #                    F
    #                     qqqqqqqqqqqqqqqqqqqqqqqqqqqqq
    #                    ------------------------------|
    #                    Λ              a
    #                    |
    #                    | cw₁
    #
    u = Matrix([       p1,w1,                        0,0          ])
    f = Matrix([  0-Mq(a),F+Fq(a),           M2+Mq(a),F2+Fq(a)   ])
    # Spring instead of q:
    # f = Matrix([        0,F-c*w1,                       M2,F2         ])
    unks =     [       p1,w1,                       M2,F2         ]

elif case == "B2.E":
    K = K_2E(EI, a)
    #                    0              1              2
    #                                    qqqqqqqqqqqqqq
    #                    |-----------------------------|
    #
    u = Matrix([        0,0,          p1,w1,           0,0         ])
    f = Matrix([       M0,F0,    0-Mq(a),0+Fq(a),M2+Mq(a),F2+Fq(a) ])
    unks =     [       M0,F0,         p1,w1,          M2,F2        ]

elif case == "B2.F":
    K = K_2E(EI, a)
    #                    0              1              2
    #                     qqqqqqqqqqqqqq
    #                    ------------------------------|
    #                    Λ
    #
    u = Matrix([       p0,0,          p1,w1,          0,0        ])
    f = Matrix([  0-Mq(a),F0+Fq(a), 0+Mq(a),0+Fq(a), M2,F2       ])
    unks =     [       p0,F0,         p1,w1,         M2,F2       ]

elif case == "B2.G":
    Mq1, Fq1 = -q/12*a*a, q/2*a
    Mq2, Fq2 = -Mq1, Fq1
    #                    0                    1        2
    #                                          qqqqqqqq
    #                    |------------------------------
    #                               2a        Λ    a
    #
    u = Matrix([          0,0,              p1,0,         p2,w2  ])
    f = Matrix([         M0,F0,             Mq1,Fq1+F1,  Mq2,Fq2 ])
    unks =     [         M0,F0,             p1,F1,        p2,w2  ]

    # Elem. 1
    l = 2*a
    l2 = l*l
    l3 = l*l*l
    K = EI/l3 * Matrix(
        [
        [  4*l2 ,   -6*l ,   2*l2 ,    6*l ,   0   ,   0   ],
        [ -6*l  ,   12   ,  -6*l  ,  -12   ,   0   ,   0   ],
        [  2*l2 ,   -6*l ,   4*l2 ,    6*l ,   0   ,   0   ],
        [  6*l  ,  -12   ,   6*l  ,   12   ,   0   ,   0   ],
        [  0    ,    0   ,   0    ,    0   ,   0   ,   0   ],
        [  0    ,    0   ,   0    ,    0   ,   0   ,   0   ],
        ]
        )
    # Elem. 2
    l = a
    l2 = l*l
    l3 = l*l*l
    K += EI/l3 * Matrix(
        [
        [  0   ,   0   ,  0    ,   0   ,  0    ,   0       ],
        [  0   ,   0   ,  0    ,   0   ,  0    ,   0       ],
        [  0   ,   0   ,  4*l2 ,  -6*l ,  2*l2 ,   6*l     ],
        [  0   ,   0   , -6*l  ,  12   , -6*l  , -12       ],
        [  0   ,   0   ,  2*l2 ,  -6*l ,  4*l2 ,   6*l     ],
        [  0   ,   0   ,  6*l  , -12   ,  6*l  ,  12       ],
        ]
        )

elif case == "B2.H":
    M, c = var("M, c")
    Fc = c*w1
    K = K_2E(EI, a/2)
    #                    0              1              2
    #                                   M              M
    #                    |------------------------------
    #                                   Λ
    #                                   |  Fc = c w₁
    #                                   |
    #
    u = Matrix([        0,0,          p1,w1,          p2,w2      ])
    f = Matrix([       M0,F0,          M,-Fc,          M,0       ])
    unks =     [       M0,F0,         p1,w1,          p2,w2      ]

elif case == "AB-1":
    K = K_1E(EI, a/2)
    #                    1                             2
    #                   F/2
    #
    #                    ◆
    #                    -------------------------------
    #                    ◆             a/2             Λ
    #
    #
    #
    u = Matrix([       0,w1,                       p2,0    ])
    f = Matrix([      M1,F/2,                       0,F2   ])
    unks =     [      M1,w1,                       p2,F2   ]

elif case == "AB-2":
    K = K_1E(EI, a)
    #                    1                             2
    #                                    |
    #                                    | F
    #                                    V            -M
    #                    |------------------------------
    #                                    a             Λ
    #
    #
    #
    # Equiv. Nodal Loads due to F:
    M, F = var("M, F")
    Me1, Fe1, Me2, Fe2 = -F*a/8, F/2, F*a/8, F/2
    #
    u = Matrix([       0,0,                        p2,0      ])
    f = Matrix([  M1+Me1,F1 + Fe1,             -M+Me2,F2+Fe2 ])
    unks =     [      M1,F1,                       p2,F2     ]


eq = Eq(K*u , f)
sol = solve(eq, unks)

for s in unks:
    pprint("\n")
    for tmp in [s, sol[s]]:
        tmp = tmp.expand()
        tmp = tmp.simplify()
        pprint(tmp)
        # pprint(latex(tmp))

        if case =="5.7" and quants == True:
            tmp = sol[s]
            tmp = tmp.subs(sub_list)
            if s == M1:
                pprint("\nM1 / Nm:")
                tmp /= Newton*m
            elif s == F1:
                pprint("\nF1 / N:")
                tmp /= Newton
            elif s == p2:
                pprint("\np2 / rad (unrealistic task):")
                tmp = tmp
            elif s == F2:
                pprint("\nF2 / N:")
                tmp /= Newton
            tmp = iso_round(tmp,"1.0")
            pprint(tmp)

if case == "B2.H":
    pprint("\nw2 for ca³ = 10 EI:")
    tmp = sol[w2]
    tmp = tmp.subs(EI, c*a**3/10)
    tmp = N(tmp,3)
    pprint(tmp)