Macaulays methode

Macaulays methode is een methode waarmee de verplaatsingsfunctie van een ligger van meerdere segmenten kan worden bepaald in één functie [Mac19]. Daarbij is het niet nodig om een balk onder te verdelen in segmenten met randvoorwaarden (afhankelijk van verdeelde belastingen, puntlasten, opleggingen, etc.) en blijft de notatie direct herkenbaar bij integratie/differentiatie (in tegenstelling tot de Dirac delta / Heaviside functies). Er bestaat een implementatie voor deze methode voor balken, maar een uitbreiding naar 2d-constructies zoals portalen ontbreekt.

Voorbeeld Macaulays methode zonder uitbreidingen

Laten we de volgende constructie bekijken:

Fig. 1 Voorbeeld constructie

Definieer belasting

Deze belasting (inclusief oplegreacties) kan als volgt worden beschreven volgens Macaulay:

\[q\left( x \right)=\underbrace{{{A}_{v}}{{\left\langle x \right\rangle }^{-1}}+{{B}_{v}}{{\left\langle x-15 \right\rangle }^{-1}}}_{\text{oplegreacties}}\underbrace{+10{{\left\langle x \right\rangle }^{0}}-10{{\left\langle x-4 \right\rangle }^{0}}}_{\text{verdeelde belasting}}\underbrace{+35{{\left\langle x-9 \right\rangle }^{-1}}}_{\text{puntlast}}\]

Dit kan in sympy worden genoteerd als:

import sympy as sym
sym.init_printing()

x = sym.symbols('x')
EI = sym.symbols('EI')
C_1, C_2, C_3, C_4 = sym.symbols('C_1 C_2 C_3 C_4') 
A_v, B_v = sym.symbols('A_v B_v')

Q = 10
L = 15
a = 4
b = 5
F = 35

q =  A_v * sym.SingularityFunction(x,0,-1) + Q * sym.SingularityFunction(x,0,0) - Q * sym.SingularityFunction(x,a,0) + F * sym.SingularityFunction(x, a+b, -1) + B_v * sym.SingularityFunction(x,L,-1)
display(q)

Definieer randvoorwaarden

Voor dit probleem gelden de volgende randvoorwaarden:

\[\begin{split}\begin{align} & w\left( 0 \right)=0 \\ & w\left( 15 \right)=0 \\ & M\left( 0 \right)=0 \\ & M\left( 15 \right)=0 \\ & V\left( {{0}^{-}} \right)=0 \\ & V\left( {{15}^{+}} \right)=0 \\ \end{align}\end{split}\]

Dit kan in sympy als volgt worden gedefinieerd:

V = -sym.integrate(q, x) + C_1
M = sym.integrate(V, x) + C_2
kappa = M / EI
phi = sym.integrate(kappa, x) + C_3
w = -sym.integrate(phi, x) + C_4
display(w)

De dwarskracht net links van \(0\) en net rechts van \(15\) wordt gemodelleerd met een ‘dummy’-variabele van \(\cfrac{1}{\infty}\).

oo = sym.Dummy('oo', prime=True)
very_small = 1/oo

Nu kunnen in SymPy de vergelijkingen voor de randvoorwaarden worden opgesteld:

Eq1 = sym.Eq(w.subs(x,0),0)
Eq2 = sym.Eq(w.subs(x,L),0)
Eq3 = sym.Eq(M.subs(x,0),0)
Eq4 = sym.Eq(M.subs(x,L),0)
Eq5 = sym.Eq(V.subs(x,0-very_small),0).subs(oo,sym.oo)
Eq6 = sym.Eq(V.subs(x,L+very_small),0).subs(oo,sym.oo)
display(Eq1, Eq2, Eq3, Eq4, Eq5, Eq6)

Los op voor randvoorwaarden

Oplossen naar de randvoorwaarden geeft:

\[w\left( x \right)=\frac{26692x\underbrace{-292{{\left\langle x \right\rangle }^{3}}-158{{\left\langle x-15 \right\rangle }^{3}}}_{\text{oplegreacties}}\underbrace{+15{{\left\langle x \right\rangle }^{4}}-15{{\left\langle x-4 \right\rangle }^{4}}}_{\text{verdeelde belasting}}\underbrace{+210{{\left\langle x-9 \right\rangle }^{3}}}_{\text{puntlast}}}{36EI}\]

welke oplossing ook met code kan worden gevonden:

sol = sym.solve((Eq1,Eq2,Eq3,Eq4,Eq5,Eq6),(C_1,C_2,C_3,C_4,A_v,B_v))
display(sol)

display(w.subs(sol).factor(EI))

Dit kan ook worden omgeschreven naar een functie per domein:

\(\displaystyle \begin{cases} \frac{79 x^{3} - 3555 x^{2} + 40781 x - 78465}{18 EI} & \text{for}\: x > 9 \\\frac{- 13 x^{3} - 360 x^{2} + 7633 x - 960}{9 EI} & \text{for}\: x > 4 \\\frac{x \left(15 x^{3} - 292 x^{2} + 26692\right)}{36 EI} & \text{for}\: x > 0\end{cases}\)

wat met de volgende code kan worden gevonden:

display(sym.simplify(w.subs(sol).rewrite(sym.Piecewise)))

Dit resultaat kan ook worden geplot voor een \(EI\) van \(10000\):

import numpy as np
import matplotlib.pyplot as plt

%config InlineBackend.figure_formats = ['svg']

w_numpy = sym.lambdify(x, w.subs(sol).subs(EI,10000).rewrite(sym.Piecewise))

x_numpy = np.linspace(0,15,100)
plt.figure()
plt.plot(x_numpy,w_numpy(x_numpy))
plt.xlabel('$x$')
plt.ylabel('$w$');
ax = plt.gca()
ax.spines['right'].set_color('none')
ax.spines['top'].set_color('none')
ax.spines['bottom'].set_position('zero')
ax.spines['left'].set_position('zero')
ax.invert_yaxis()

The dwarskrachtverdeling kan ook worden gevonden:

\[V\left( x \right)=\underbrace{\frac{146{{\left\langle x \right\rangle }^{0}}}{3}+\frac{79{{\left\langle x-15 \right\rangle }^{0}}}{3}}_{\text{oplegreacties}}\underbrace{-10{{\left\langle x \right\rangle }^{1}}+10{{\left\langle x-4 \right\rangle }^{1}}}_{\text{verdeelde belasting}}\underbrace{-35{{\left\langle x-9 \right\rangle }^{0}}}_{\text{puntlast}}\]

met de volgende code:

display(V.subs(sol))

En ook deze formule kan worden omgeschreven naar een formule per domein:

\[\begin{split}\displaystyle \begin{cases} - \frac{79}{3} & \text{for}\: x > 9 \\\frac{26}{3} & \text{for}\: x > 4 \\\frac{146}{3} - 10 x & \text{for}\: x > 0 \end{cases}\end{split}\]

display(sym.simplify(V.subs(sol).rewrite(sym.Piecewise)))

en kan ook worden geplot:

V_numpy = sym.lambdify(x, V.subs(sol).rewrite(sym.Piecewise))

x_numpy = np.linspace(0,15.01,10000)
plt.figure()
plt.plot(x_numpy,V_numpy(x_numpy))
plt.xlabel('$x$')
plt.ylabel('$V$');
ax = plt.gca()
ax.spines['right'].set_color('none')
ax.spines['top'].set_color('none')
ax.spines['bottom'].set_position('zero')
ax.spines['left'].set_position('zero')
ax.invert_yaxis()

Voorbeeld Macaulays methode in SymPy

De methode van Macaulay is ook geïmplementeerd in SymPy en kan daar worden toegepast zonder vergelijkingen op te moeten stellen.

Let op: het assenstelsel is in SymPy gedefinieerd met de verticale as positief naar boven.

import sympy as sym
from sympy.physics.continuum_mechanics.beam import Beam
sym.init_printing()

E, I = sym.symbols('E, I')
b = Beam(15, E, I)
A_v = b.apply_support(0,type='pin')
B_v = b.apply_support(15,type='roller')
b.apply_load(-10, 0, 0, 4)
b.apply_load(-35, 9, -1)
b.solve_for_reaction_loads(A_v, B_v)
b.deflection().simplify()

b.plot_deflection({E : 10000, I : 1});

b.plot_shear_force();