Gelöste Aufgaben/UEBH

Aufgabenstellung

Dies ist eine Näherungslösung für einen Standard-Lastfall des Euler-Bernoulli-Balkens. Hier geht es also darum, die Näherungslösung mit der analytisch richtigen Lösung zu vergleichen.

Der Euler-Bernoulli-Balken AB wird durch ein Moment M zwischen den beiden gelenkigen Lagern belastet. 

Lageplan

Gesucht ist die Biegelinie mit dem Ansatz von Ritz und zwei Trial-Funktionen.

Lösung mit Maxima

Beim Verfahren von Ritz arbeiten wir mit

• dem Prinzip vom Minimum der Potentiellen Energie und
• Ansatzfunktionen über die gesamte Länge des Balkens.

Header

Wir berechnen die Potentielle Energie U des Systems in Abhängigkeit von den generalisierten Koordinaten W_i und erhalten aus

$\displaystyle {\frac {d\,U}{d\,W_{i}}}{\stackrel {!}{=}}0$

die Gleichung für den gesuchten Koeffizienten W_i der Trial-Funktionen.

/*******************************************************/
/* MAXIMA script                                       */
/* version: wxMaxima 15.08.2                           */
/* author: Andreas Baumgart                            */
/* last updated: 2018-02-16                            */
/* ref: TMC, Labor 3                                   */
/* description: Ritz approach to EBB, load-case 5      */
/*                                                     */
/*******************************************************/

/* declare variables and functions */
declare("Π", alphabetic); /* strain energy */
assume(l>0);

Declarations

Um die Lösung dimensionslos zu machen, nutzen wir die analytische Lösung des Problems , hier die Beträge der maximalen Auslenkung des Balkens und der Verdrehung des Balkens am Momenten-Angriffspunkt hier für a = ℓ/2:

${\begin{array}{lll}{\hat {W}}&=\displaystyle {\frac {{\sqrt {3}}\,\ell ^{2}}{8\,\cdot E\,I}}\cdot M&\ldots {\text{ die maximale Auslenkung des Balkens für }}a=\ell /2,\\{\hat {\Phi }}&=\displaystyle {\frac {\ell }{4\cdot E\,I}}\cdot M&\ldots {\text{ die Verdrehung des Balkens am Momenten-Angriffspunkt für }}a=\ell /2\end{array}}$.

Mit diesen Skalierungsfaktoren muss für α=1/2 die maximale Auslenkung und die maximale Verdrehung im Momenten-Angriffspunkt nahezu eins sein.

Dimensionslose Orts-Koordinaten sind

${\begin{array}{ll}x&=\xi \cdot \ell ,\\a&=\alpha \cdot \ell .\end{array}}$.

/*anayltic solution */
analytic: w(xi) = M*l^2/(2*EI)*(xi^3/3+alpha^2*xi-2*alpha*xi+(2*xi)/3
/*                                   foeppel-function         */
                       - (if xi<alpha then 0 else xi-alpha)^2);

/* make dim'less with load case #5 */
dimless:[[W   = M*l^2/(72*sqrt(3)*EI)*q[w],
          Phi = M*l/(12*EI)*q[phi]],
         [x=xi*l,
          a=alpha*l]];

Formfunctions

Gesucht ist hier ein Ansatz mit zwei Trial Functions, also

$w(x)=W_{1}\cdot \phi _{1}(x)+W_{2}\cdot \phi _{2}(x)$

und den gesuchten Koeffizienten W₁ und W₂. Die Formfunktionen ϕ_i(x) müssen dabei alle geometrischen Randbedingungen erfüllen, also:

${\begin{array}{ll}w(0)&=0,\\w(\ell )&=0.\end{array}}$

Die Funktionen müssen nicht die Randbedingungen für Schnittmomente in A und B erfüllen - schließlich suchen wir nach einer Näherungslösung!

Um die Funktion ϕ_i(x) zu konstruieren, gibt es zwei Ansätze:

1. Wir setzen die Funktionen anschaulich aus Nullstellen x_N zusammen: dazu nehmen wir das Produkt aus (x_N-x), hier mit mindestens einer Nullstelle in x_N=0 und in x_N=ℓ. In dimensionslosen Koordinaten ist dann

   $\phi _{1}(\xi )=\xi \cdot (1-\xi ){\text{ mit }}\xi =\displaystyle {\frac {x}{\ell }}$.

   Die zweite Funktion gewinnen wir, indem wir dieser Funktion eine weitere Nullstelle - z.B. im Kraft-Angriffspunkt - hinzufügen, also mit  (α-ξ) multiplizieren:

   $\phi _{2}(\xi )=\xi \cdot (1-\xi )\cdot (\alpha -\xi )$.

   Die geforderten Nullstellen für die Erfüllung der geometrischen Randbedingungen bleiben so erhalten und wir haben eine weitere, linear von ϕ₁(x) unabhängige Ansatzfunktion. Mit dieser Formulierung könnten wir nun arbeiten. Die gesuchten Koeffizienten W₁ und 'W₂'  ergeben sich dann aus der Gleichgewichtsbedingung des Prinzip vom Minimum der Potentiellen Energie. Allerdings machen diese Koeffizienten keinen "Sinn" - wir können ihnen im Problem keine anschauliche Bedeutung zumessen. Damit wir die Koeffizienten aber auch anschaulich deuten können, skalieren wir die ϕ_i(x) so, dass

   ${\begin{array}{ll}\phi _{1}(\displaystyle {\frac {1}{2}})&=1{\text{ und}}\\\displaystyle {\frac {d\phi _{2}}{dx}}|_{x=a}&=1\end{array}}$,

   also

   $\displaystyle {{C}_{1}}=4\;\;{{C}_{2}}=-4\ell$.

   Das ist schon brauchbar: bei dieser Variante der Formulierung denken wir uns die Näherungslösung als Linearkombination einer jeweils - bzgl. x = ℓ/2 - punkt-symmetrischen und achsen-symmetrischen Trial-Function zusammengesetzt. Diesen Ansatz verfolgen wir in UEBO weiter. Wir können einen Schritt weiter zur Formulierung von anschaulich "denkbaren" Ansatzfunktionen gehen - hier im zweiten Weg:
2. Dabei gehen wir von der  - zunächst - allgemeinen Formulierung

   $w(\xi )=a_{3}\cdot \xi ^{3}+a_{2}\cdot \xi ^{2}+a_{1}\cdot \xi +a_{0}$

   mit den unbekannten a_i aus. Von diesen vier verbrauchen wir zwei für die geometrischen Randbedingungen und zwei bleiben frei für die Wahl von anschaulichen Koordinaten der Verformung. Dazu wählen wir  

   ${\begin{array}{ll}w(0)&=0\\w(\ell )&=0\\w(a)&=W\\\displaystyle {\frac {dw}{dx}}|_{x=a}&=\Phi \end{array}}$.

   

   Trial-Functions ϕ(x)

   Hier ergibt sich

   ${\begin{array}{l}\displaystyle \phi _{1}=-{\frac {\left(3\cdot {{\alpha }^{2}}-2\cdot \alpha \right)\cdot \xi +\left(1-3\cdot {{\alpha }^{2}}\right)\cdot {{\xi }^{2}}+\left(2\cdot \alpha -1\right)\cdot {{\xi }^{3}}}{{{\alpha }^{4}}-2\cdot {{\alpha }^{3}}+{{\alpha }^{2}}}}\\\displaystyle \phi _{2}={\frac {\alpha \cdot \ell \cdot \xi +\left(-\alpha -1\right)\cdot \ell \cdot {{\xi }^{2}}+\ell \cdot {{\xi }^{3}}}{{{\alpha }^{2}}-\alpha }}\end{array}}$.

   Für x=a ist die Auslenkung w(a) = W und die Steigung w'(a) = Φ - die beiden Eigenschaften sind in den Trial-Functions vollständig entkoppelt. Und mit diesen etwas komplizierten ϕs arbeiten wir im folgenden. So sehen sie aus für  α=7/10.

/* derive trial-function(s) */
trial: w(x) = sum(c[i]*x^i,i,0,3);
GBC: [subst([x=0],     subst(trial,w(x))   )= 0,
      subst([x=l],     subst(trial,w(x))   )= 0,
      subst([x=a],     subst(trial,w(x))   )= W,
      subst([x=a],diff(subst(trial,w(x)),x))= Phi];

C: makelist(c[i],i,0,3);
sol[1] : solve(GBC,C)[1];
trial: ratsimp(subst(sol[1],trial));

/* trial function printout .... */
trials : ratsimp(subst(dimless[2],[phi[1]=coeff(expand(subst(trial,w(x))),W),
                                   phi[2]=coeff(expand(subst(trial,w(x))),Phi)]));
print(trials)$
preamble: "set yrange [] reverse";
plot2d(subst([alpha=7/10],subst(trials,[phi[1],phi[2]/l])),[xi,0,1],
  [legend, "ϕ1","ϕ2/ℓ"],
  [gnuplot_preamble, preamble],
  [title, "trial-functions für α=7/10"],
  [xlabel, "ξ →"],
  [ylabel, "ϕ →"] );

Potential Energy

Für die Gleichgewichtsbedingungen setzten wir Π (aus Abschnitt Euler-Bernoulli-Balken) und A in U ein und schreiben die skalare Gleichung allgemein in Matrizenform an. Dabei müssen wir

$\displaystyle {\frac {d\phi }{x}}={\frac {d\phi }{\xi }}\cdot \underbrace {\displaystyle {\frac {d\xi }{x}}} _{\displaystyle ={\frac {1}{\ell }}}$

berücksichtigen und erhalten mit

$A=M\cdot \Phi$

das Potential in Matrix-Schreibweise:

$U=\displaystyle {\frac {1}{2}}\cdot \displaystyle {\underline {Q}}^{T}\cdot {\underline {\underline {A}}}\cdot {\underline {Q}}-{\underline {Q}}^{T}\cdot {\underline {b}}$.

In der Spaltenmatrix der unbekannten Koordinaten

${\underline {Q}}=\left({\begin{array}{c}W\\\Phi \end{array}}\right)$

stehen nun die zwei Unbekannten W und Φ. Einsetzen der Ansatzfunktion in die Formänderungsenergie und die Arbeitsfunktion liefern für A und b:

${\underline {\underline {A}}}=2\;EI\cdot {\begin{pmatrix}\displaystyle {\frac {\left(2-12\cdot \alpha +30\cdot {{\alpha }^{2}}-36\cdot {{\alpha }^{3}}+18\cdot {{\alpha }^{4}}\right)}{\left({{\alpha }^{8}}-4\cdot {{\alpha }^{7}}+6\cdot {{\alpha }^{6}}-4\cdot {{\alpha }^{5}}+{{\alpha }^{4}}\right)\cdot {\ell ^{3}}}}&\displaystyle -{\frac {\left(-2+7\cdot \alpha -9\cdot {{\alpha }^{2}}+6\cdot {{\alpha }^{3}}\right)}{\left({{\alpha }^{6}}-3\cdot {{\alpha }^{5}}+3\cdot {{\alpha }^{4}}-{{\alpha }^{3}}\right)\cdot {\ell ^{2}}}}\\\displaystyle -{\frac {\left(-2+7\cdot \alpha -9\cdot {{\alpha }^{2}}+6\cdot {{\alpha }^{3}}\right)}{\left({{\alpha }^{6}}-3\cdot {{\alpha }^{5}}+3\cdot {{\alpha }^{4}}-{{\alpha }^{3}}\right)\cdot {\ell ^{2}}}}&\displaystyle {\frac {\left(2-2\cdot \alpha +2\cdot {{\alpha }^{2}}\right)}{\left({{\alpha }^{4}}-2\cdot {{\alpha }^{3}}+{{\alpha }^{2}}\right)\cdot \ell }}\end{pmatrix}}$,

${\underline {b}}=\left({\begin{array}{c}0\\M\end{array}}\right)$.

✔ Arbeiten mit Matlab:

In numerischen Anwendungen - zum Beispiel mit Matlab - kommen wir mit der parametrischen Darstellung in α nicht durch. Wir können auch nicht erst das Potential aufschreiben und durch Ableiten die Gleichgewichtsbedingungen anschreiben. Da hilft es, sich klar zu machen, dass die Koeffizienten der System-Matrix A aus $\displaystyle a_{i,j}=EI\cdot \int _{0}^{\ell }\phi _{i}''\cdot \phi _{j}''\;dx$ kommen. Einfach wird es, wenn die Trial-Functions Polynome sind (bei uns sind sie deshalb immer Polynome). Dann bieten die Matlab-Werkzeuge zur Polynom-Manipulation ein praktisches Instrumentarium, um Ableitungen (polyder), Produkte (conv), Integrale (polyint) und Funktionswerte (polyval) zu bestimmen. Dieser Ansatz liefert z.B. $\displaystyle {{a}_{1,2}}=-{\frac {\left(-4+14\cdot \alpha -18\cdot {{\alpha }^{2}}+12\cdot {{\alpha }^{3}}\right)\cdot {\mathit {EI}}}{{{\alpha }^{6}}-3\cdot {{\alpha }^{5}}+3\cdot {{\alpha }^{4}}-{{\alpha }^{3}}}}$. Und mit dieser Überlegung ist auch sofort klar, warum die System-Matrizen für diese Näherungslösungen immer symmetrisch sein müssen! (Ist es?)

/* define potential energy of system */
PMPE : [U = Π - A,
        Π = 1/2*EI*'integrate('diff(w(x),x,2)^2,x,0,l),
        A = M*Phi(a)];

PMPE: subst(Phi(a) = subst([x=a],diff(subst(trial,w(x)),x)),
  subst(trial,
    subst(PMPE[3],subst(PMPE[2], PMPE[1]))));
PMPE : expand(ev(PMPE,nouns));

/* print system matrices */
print("A = ",  2*ratsimp(subst(dimless[2],matrix([coeff(subst(PMPE,U),W,2),coeff(coeff(subst(PMPE,U),W),Phi)/2],
                                                 [coeff(coeff(subst(PMPE,U),W),Phi)/2,coeff(subst(PMPE,U),Phi,2)]))))$
print("b = ",- ratsimp(matrix([coeff(coeff(subst(PMPE,U),W,1),Phi,0)],[coeff(coeff(subst(PMPE,U),Phi,1),W,0)])))$

Equilibrium Conditions

Diese Gleichung erfüllt die Gleichgewichtsbedingungen

$\displaystyle {\frac {dU}{dW}}{\stackrel {!}{=}}0{\text{ und }}\displaystyle {\frac {dU}{d\Phi }}{\stackrel {!}{=}}0$,

wenn

${\begin{array}{ll}\left(-12\cdot {{\alpha }^{5}}+30\cdot {{\alpha }^{4}}-32\cdot {{\alpha }^{3}}+18\cdot {{\alpha }^{2}}-4\cdot \alpha \right)\cdot \ell \cdot EI\Phi +\left(36\cdot {{\alpha }^{4}}-72\cdot {{\alpha }^{3}}+60\cdot {{\alpha }^{2}}-24\cdot \alpha +4\right)\cdot EI\cdot W&=0\\\left(-4\cdot {{\alpha }^{4}}+8\cdot {{\alpha }^{3}}-8\cdot {{\alpha }^{2}}+4\cdot \alpha \right)\cdot \ell \cdot \Phi +\left(12\cdot {{\alpha }^{3}}-18\cdot {{\alpha }^{2}}+14\cdot \alpha -4\right)\cdot EI\cdot W&=-\left({{\alpha }^{6}}-3\cdot {{\alpha }^{5}}+3\cdot {{\alpha }^{4}}-{{\alpha }^{3}}\right)\cdot {{\ell }^{2}}\cdot M\end{array}}$.

/* equilibrium condition */
equcon: ['diff(U,W) = 0,'diff(U,Phi) = 0]$
equcon: subst(PMPE,equcon)$
equcon: ev(equcon,nouns)$
equcon: ratsimp(subst(dimless[2],equcon));

Solving

Auflösen der Gleichungen nach den unbekannten Koordinaten W und Φ liefert

${\begin{array}{ll}W&\displaystyle ={\frac {\left(2\cdot \alpha -9\cdot {{\alpha }^{2}}+16\cdot {{\alpha }^{3}}-15\cdot {{\alpha }^{4}}+6\cdot {{\alpha }^{5}}\right)\cdot {{\ell }^{2}}\cdot M}{6\cdot EI}},\\\Phi &\displaystyle ={\frac {\left(1-6\cdot \alpha +15\cdot {{\alpha }^{2}}-18\cdot {{\alpha }^{3}}+9\cdot {{\alpha }^{4}}\right)\cdot \ell \cdot M}{3\cdot EI}}\end{array}}$.

Noch einfacher wird das Ergebnis mit den dimensionslosen Koordinaten

${\begin{array}{l}\displaystyle q_{W}={\frac {W}{\hat {W}}}\\\displaystyle q_{\Phi }={\frac {\Phi }{\hat {\Phi }}}\end{array}}$

zu

${\begin{array}{ll}{{q}_{W}}&=8\cdot {{3}^{\frac {5}{2}}}\cdot {{\alpha }^{5}}-20\cdot {{3}^{\frac {5}{2}}}\cdot {{\alpha }^{4}}+64\cdot {{3}^{\frac {3}{2}}}\cdot {{\alpha }^{3}}-4\cdot {{3}^{\frac {7}{2}}}\cdot {{\alpha }^{2}}+8\cdot {{3}^{\frac {3}{2}}}\cdot \alpha ,\\{{q}_{\Phi }}&=36\cdot {{\alpha }^{4}}-72\cdot {{\alpha }^{3}}+60\cdot {{\alpha }^{2}}-24\cdot \alpha +4\end{array}}$.

/* solve */
sol[2] : ratsimp(solve(equcon,[W,Phi]))[1];
/* ... or in dimensionless coordinates */
sol[3] : ratsimp(solve(subst(dimless[1],equcon),[q[w],q[phi]]))[1];
/* approximated solution */
ritz : w(xi)=ratsimp(subst(sol[2],subst(dimless[2],rhs(trial))));

Post-Processing

Parameterstudie: q_W und q_ϕ über α

Die gesuchten Koordinaten  q_W und q_Φ sind dimensionslos. Wir können sie direkt für verschiedene Werte von α auftragen.

Dies ist nicht die Biegelinie des Balkens! Man sieht aber - mit etwas Übung - direkt aus dieser Darstellung die Charakteristika der Lösung:

• die Punkt-Symmetrie der Lösung für α=1/2,
• die schlechte Approximation der Lösung für α=1/2:  q_Φ sollte dort nahe "1" sein, um die analytische Lösung gut abzubilden.

Und so sieht die normierte Biegelinie des Balkens im vergleich von Ritz-Näherung zu analytischer Lösung für verschiedene Werte von a aus:

Parameterstudie Biegelinie

Die dicken Linien gehören zu Näherung nach dem Ritz-Ansatz, die dünnen zur analytischen Lösung. Besonders für α=1/2 liefert der Ritz-Ansatz kein überzeugendes Ergebnis. Hier müssten wir mehr Trial-Functions "spendieren".

✔ Schnittmomenten-Verlauf:

Probieren Sie aus, wie gut der Ritz-Ansatz den Verlauf des Schnittmoments approximiert!

/* post-processing */
/* plot solutions */
preamble: "set yrange [] reverse";

plot2d([subst(sol[3],q[w]),subst(sol[3],q[phi])],[alpha,0,1],
  [legend, "q[w]", "q[Φ]"],
  [gnuplot_preamble, preamble],
  [xlabel, "α →"],
  [ylabel, "q[w],q[Φ] →"] );

/* abbreviation for max. displacement */
wH : subst([q[w]=1],subst(dimless[1],W));

/* plot w(x) for different αs */
/* and compare analytic with Ritz-solution */

leg : append([legend], makelist(simplode (["α = ",i,"/10"]),i,1,9,2));
toPlot : [makelist(subst([alpha=i/10],ratsimp(subst(ritz,w(xi))/wH)), i,1,9,2),
          makelist(subst([alpha=i/10],ratsimp(subst(analytic,w(xi))/wH)), i,1,9,2)];
plot2d(append(toPlot[1],toPlot[2]),[xi,0,1],
  leg,
  [color, green, blue, red, magenta, black],
  [style, [lines,3], [lines,3], [lines,3], [lines,3], [lines,3], [lines,1], [lines,1], [lines,1], [lines,1], [lines,1]],
  [gnuplot_preamble, preamble],
  [xlabel, "x/l →"],
  [ylabel, "w(x)/W →"] );

---

Links

• ...

Literature

• ...