Gelöste Aufgaben/UEBO

Aufgabenstellung

Diese Problemstellung liefert einen Näherungsansatz für eine Standardlösung zum Euler-Bernoulli-Balken.

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

Lageplan

Gesucht ist eine Lösung für die Biegelinie mit dem Ansatz von Ritz und zwei Trial-Funktionen.

(Weg "1" wie in UEBH beschrieben.)

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

$\frac{dU}{d{W}_i}\stackrel{!}{=}0$

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

📑 toggle code

/*******************************************************/
/* MAXIMA script                                       */
/* version: wxMaxima 15.08.2                           */
/* author: Andreas Baumgart                            */
/* last updated: 2019-10-13                            */
/* ref: TMC, Labor 4                                   */
/* description: Ritz approach to EBB, load-case 5      */
/*                                                     */
/*******************************************************/

/* declare variables and functions */
declare("Π", alphabetic); /* strain energy */
assume(ℓ>0);
dimless:[x=xi*ℓ,
         a=alpha*ℓ];

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 für a = ℓ und der Verdrehung des Balkens am Momenten-Angriffspunkt für a = ℓ/2:

$W^{max}=\frac{M{\ell }^2}{9\sqrt{3}\cdot EI}$	die maximale Auslenkung des Balkens für a=ℓ
${\mathrm{\Phi }}^M=\frac{M\ell }{12\cdot EI}$	die Verdrehung des Balkens am Momenten-Angriffspunkt für a=ℓ/2

Dimensionslose Orts-Koordinaten sind

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

📑 toggle code

/*analytic solution vgl. Lexikon/Euler-Bernoulli-Blaken */
analytic: w(xi) = M*ℓ^2/(6*EI)*(xi^3+xi*(2-6*alpha+3*alpha^2)
/*                                   foeppel-function         */
                       - 3*(if xi<alpha then 0 else xi-alpha)^2);
sectionI : M*ℓ^2/(6*EI)*(xi^3+xi*(2-6*alpha+3*alpha^2));
/* make dim'less with load case #5 */
Phi[C] : subst([xi=1/2],diff(subst([alpha=1/2],sectionI),xi)/ℓ);
W[max] : -subst([alpha=1],subst(solve(
                           diff(subst([alpha=1],sectionI),xi)/ℓ=0,xi)[2],sectionI));

Formfunctions

Bei der Suche nach passenden Trial-Functions ϕ lassen wir uns ebenfalls von der analytischen Lösung des Problems "inspirieren":

$w_a=\frac{M{\ell }^2}{3EI}((3\cdot {\alpha }^2-6\cdot \alpha +2)\cdot \xi +{\xi }^3-3<\xi -\alpha {>}^2)$

Der Funktionsverlauf von w^a hat zwei charakteristische Ausprägungen:

... für a=0	... für a=ℓ/2
Diese Lösung - mit dem angreifenden Moment in A - hat eine starke symmetrische Komponente bzgl der Stab-Mitte.	Diese Lösung - mit dem angreifenden Moment in der Stab-Mitte - ist punktsymmetrisch zum Momenten-Angriffspunkt.

Und so wählen wir unsere Trial-Functions als

${\varphi }_1=c_1\cdot \xi \cdot (1-\xi )\text{ und }{\varphi }_2=c_2\cdot \xi \cdot (1-\xi )\cdot (\alpha -\xi )$.

Trial-Functions

Für α=7∙ℓ/10 sehen sie so aus;

Die Koeffizienten c₁ und c₂ haben wir dabei so gewählt, dass

$\begin{array}{ll}{\varphi }_1(\frac{1}{2})& =1\\ {{\varphi }_2}^{\prime }(\alpha )& =1\end{array}$.

Mit den neuen, gesuchten Wichtungsfaktoren q_w und q_ϕ ist die Ansatzfunktion zur Lösung mit dem Verfahren von Rayleigh-Ritz damit

$w(\xi )=\frac{M{\ell }^2}{3EI}(q_w\left(4\frac{\xi -{\xi }^2}{3\sqrt{3}}\right)+q_{\varphi }(-{\xi }^3+(1+\alpha ){\xi }^2-\alpha \xi ))$

Aufgrund der gewählten Skalierungsfaktoren erwarten wir als Ergebnis näherungsweise

• für α=½: q_w ≈ 0 und q_ϕ ≈ 1,
• für α= 0: q_w ≈ 1 und q_ϕ ≈ 0.

📑 toggle code

/* derive trial-function(s) */
phi : [c[1]*xi*(1-xi), c[2]*xi*(1-xi)*(alpha-xi)];
GBC: [subst([xi = 1/2 ],      phi[1]      )= 1,
      subst([xi = 1/2 ], diff(phi[2],xi)/ℓ)= 1];
sol[1] : solve(GBC,[c[1],c[2]])[1];
phi: ratsimp(subst(sol[1],phi));

ansatz: [w(xi) = q[w]*W[max]*phi[1] + q[p]*Phi[C]*phi[2]];

preamble: "set yrange [] reverse";
plot2d(subst([alpha=7/10],[phi[1],phi[2]/ℓ]),[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

$\frac{d\varphi }{x}=\frac{d\varphi }{\xi }\cdot \underset{=\frac{1}{\ell }}{\underbrace{\frac{d\xi }{x}}}$

berücksichtigen und erhalten mit der Arbeitsfunktion des Moments

$A=M\cdot w^{\prime }{|}_{\xi =\alpha }$

das Potential in Matrix-Schreibweise:

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

wobei

$\underset{\_}{Q}=\left(\begin{array}{c}{q}_w\\ q_{\varphi }\end{array}\right)$.

Einsetzen der Ansatzfunktion in die Formänderungsenergie und die Arbeitsfunktion liefert für die Matrizen A und b:

$\underset{\_}{\underset{\_}{A}}=\left(\begin{array}{cc}\frac{128}{81}& \frac{16}{3^{\frac{5}{2}}}-\frac{32\cdot \alpha }{3^{\frac{5}{2}}}\\ \frac{16}{3^{\frac{5}{2}}}-\frac{32\cdot \alpha }{3^{\frac{5}{2}}}& \frac{8\cdot {\alpha }^2}{3}-\frac{8\cdot \alpha }{3}+\frac{8}{3}\end{array}\right)$,

$\underset{\_}{b}=\left(\begin{array}{c}\frac{8}{3^{\frac{3}{2}}}-\frac{16\cdot \alpha }{3^{\frac{3}{2}}}\\ 2\cdot \alpha -2\cdot {\alpha }^2\end{array}\right)$.

📑 toggle code

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

PMPE: subst([Phi(a) = subst([xi=alpha],diff(subst(ansatz,w(xi))/ℓ,xi))],
      subst(ansatz,
      subst(PMPE[3],subst(PMPE[2], PMPE[1]))));
PMPE : expand(ev(PMPE,nouns));

/* unknowns */
Q : [q[w],q[p]];

Equilibrium Conditions

Diese Gleichung erfüllt die Gleichgewichtsbedingungen

$\frac{dU}{d{Q}_i}\stackrel{!}{=}0$,

wenn

$\left(\begin{array}{cc}\frac{128}{81}& \frac{16}{3^{\frac{5}{2}}}-\frac{32\cdot \alpha }{3^{\frac{5}{2}}}\\ \frac{16}{3^{\frac{5}{2}}}-\frac{32\cdot \alpha }{3^{\frac{5}{2}}}& \frac{8\cdot {\alpha }^2}{3}-\frac{8\cdot \alpha }{3}+\frac{8}{3}\end{array}\right)\cdot \left(\begin{array}{c}{q}_w\\ q_p\end{array}\right)=\left(\begin{array}{c}\frac{8}{3^{\frac{3}{2}}}-\frac{16\cdot \alpha }{3^{\frac{3}{2}}}\\ 2\cdot \alpha -2\cdot {\alpha }^2\end{array}\right)$.

📑 toggle code

/* equilibrium condition */
eom : makelist(expand(diff(subst(PMPE,U)/(M^2*ℓ/(6*EI)),Q[i])),i,1,2);

A :   funmake('matrix, makelist(makelist(coeff(eom[i],Q[j]),j,1,2),i,1,2));
b : - funmake('matrix, makelist([subst(makelist(Q[j]=0,j,1,2),eom[i])],i,1,2));

print(A,"∙",transpose(Q),"=",b)$

Solving

Auflösen der Gleichungen nach den unbekannten Koordinaten q_w und q_ϕ liefert

$\begin{array}{ll}{q}_w& =-\frac{3^{\frac{3}{2}}}{8}(-2+7\alpha -9{\alpha }^2+6{\alpha }^3)\\ q_{\varphi }& =-\frac{1}{2}(1-6\alpha +6{\alpha }^2)\end{array}$.

Damit ist die gesuchte Näherungs-Lösung

$w(\xi )=\frac{M{\ell }^2}{6EI}((3\cdot {\alpha }^2-6\cdot \alpha +2)\cdot \xi +(-9\cdot {\alpha }^2+12\cdot \alpha -3)\cdot {\xi }^2+(6\cdot {\alpha }^2-6\cdot \alpha +1)\cdot {\xi }^3)$.

📑 toggle code

/* 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

Verlauf der Koordinaten q_w, q_ϕ

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

Wir sehen:

• für α=½: die Lösung wird - wie erwartet - nur durch ϕ₂ beschreiben - also q_w ≈ 0 und q_ϕ ≈ 1; allerdings ist die Qualität der Lösung mit q_ϕ = 1/4 sehr schlecht - hier drückt der Sprung in der Momenten-Kennlinie der analytischen Lösung auf das Ergebnis (s.u.).
• für α= 0: die Lösung wird - wie erwartet - primär durch ϕ₁ beschreiben, also q_w ≈ 1 und q_ϕ ≈ 0. Hier zeigt die Lösung mit q_w = 1.3 und q_ϕ = -0.5 einen recht großen Lösungs-Anteil der punktsymmetrischen Trial-Function.

Parameterstudie Biegelinie

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:

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

Vergleich: analytische / numerische Lösung für den Biegemomenten-Verlauf.

📑 toggle code

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

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

/* 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))/W[max])), i,1,9,2),
          makelist(subst([alpha=i/10],ratsimp(subst(analytic,w(xi))/W[max])), 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/ℓ →"],
  [ylabel, "w(x)/W →"] );

Post-Processing - Nachtrag

Wieso die Näherungslösung - besonders für α=½ - so schlecht ist, erkennt man beim Auftragen der Biegemomente im Stab für

• die analytische Lösung

$\frac{M(x)}{M}=\{\begin{array}{cl}-\xi & \dots \text{ für }\xi <\frac{1}{2}\\ 1-\xi & \dots \text{ sonst. }\end{array}$ und

• die numerische Lösung

$\frac{M(x)}{M}=\frac{1}{4}(2\xi -1)$.

📑 toggle code

/* Vergleich der Momenten-Kennlinien */
displacements : M*ℓ^2/(6*EI)*[[(xi^3+xi*(2-6*alpha+3*alpha^2)), (xi^3+xi*(2-6*alpha+3*alpha^2) - 3*(xi-alpha)^2)],
                             [((3*alpha^2-6*alpha+2)*xi+(-9*alpha^2+12*alpha-3)*xi^2+(6*alpha^2-6*alpha+1)*xi^3)]];
bendingmoments: subst([xi=t],subst([alpha=1/2],-EI*diff(displacements,xi,2)/ℓ^2/M));

plot2d([[parametric, t, bendingmoments[1][1], [t, 0, 1/2]],
        [parametric, t, bendingmoments[1][2], [t, 1/2, 1]],
        [parametric, t, bendingmoments[2][1], [t,  0 , 1]]],
        [color, red, red, blue],
        [style, [lines,1], [lines,1], [lines,3]],
        [legend, "analytic, section I", "analytic, section II", "Rayleigh-Ritz"],
        [xlabel, "x/ℓ →"],
        [ylabel, "M(x)/M →"] );

---

Links

• ...

Literature

• ...