Gelöste Aufgaben/Kw53: Unterschied zwischen den Versionen

Aus numpedia
Zur Navigation springen Zur Suche springen
Keine Bearbeitungszusammenfassung
Keine Bearbeitungszusammenfassung
 
(8 dazwischenliegende Versionen desselben Benutzers werden nicht angezeigt)
Zeile 14: Zeile 14:
<onlyinclude>
<onlyinclude>
[[Datei:Kw50-01.png|300px|left|mini|Lageplan (wie [[Gelöste Aufgaben/Kw50|Kw50]])]]
[[Datei:Kw50-01.png|300px|left|mini|Lageplan (wie [[Gelöste Aufgaben/Kw50|Kw50]])]]
Geben Sie die Lösung für ein Euler-Bernoulli-Modell der Brücke mit dem Verfahren von [[Randwertprobleme/Methoden zur Lösung von Randwertproblemen/Verfahren von Rayleigh-Ritz (EBB)|Rayleigh-Ritz (EBB)]] an - hier mit [https://de.wikipedia.org/wiki/Lagrange-Multiplikator Lagrange-Multiplikator] für die geometrische Zwangsbedingung.
Geben Sie die Lösung für ein Euler-Bernoulli-Modell der Brücke mit dem Verfahren von [[Randwertprobleme/Methoden zur Lösung von Randwertproblemen/Verfahren von Rayleigh-Ritz (EBB)|Rayleigh-Ritz (EBB)]] an - hier ohne[https://de.wikipedia.org/wiki/Lagrange-Multiplikator Lagrange-Multiplikator] für die geometrische Zwangsbedingung.


Dies ist eine Näherungslösung zu [[Gelöste Aufgaben/Kw50|Kw50]].
Dies ist eine Näherungslösung zu [[Gelöste Aufgaben/Kw50|Kw50]].
Zeile 20: Zeile 20:
Ermitteln Sie die genäherten Verläufe der Schnittgrößen und Verschiebungen im Balken für diese Parameter:
Ermitteln Sie die genäherten Verläufe der Schnittgrößen und Verschiebungen im Balken für diese Parameter:


<math>\begin{array}{ll}K_C =&\displaystyle  5 \frac{E\,I}{\ell_0}\\m_A =&\displaystyle \frac{m_B}{5}  \end{array}</math>
::<math>\begin{array}{ll}K_C =&\displaystyle  5 \frac{E\,I}{\ell_0}\\m_A =&\displaystyle \frac{m_B}{5}  \end{array}</math>


== Lösung mit Maxima ==
== Lösung mit Maxima ==
Zeile 28: Zeile 28:


<!-------------------------------------------------------------------------------->
<!-------------------------------------------------------------------------------->
{{MyCodeBlock|title=Header
{{MyCodeBlock|title=Header
|text=
|text=
Für die Lösung nutzen wir hier das Ritz-Verfahren mit einer Variante:
Für die Lösung nutzen wir hier das Ritz-Verfahren. Die geometrischen Zwangsbedingungen arbeiten wir direkt in die Trial-Functions ein.
 
* die geometrische Zwangsbedingung für die Punkte ''A, B'' durch das Seil erfassen wir durch einen [https://de.wikipedia.org/wiki/Lagrange-Multiplikator Lagrange-Multiplikator].
|code=
|code=
<syntaxhighlight lang="lisp" line start=1>
<syntaxhighlight lang="lisp" line start=1>
Zeile 41: Zeile 38:
/* author: Andreas Baumgart                            */
/* author: Andreas Baumgart                            */
/* last updated: 2019-02-12                            */
/* last updated: 2019-02-12                            */
/* ref: TM-C, Labor 1, dimensionless representation    */
/* ref: TM-C, Brigde-Problem                          */
/* description: finds the rayleigh-ritz with Lagragian */
/* description: finds the plain rayleigh-ritz solution */
/*              Multiplyers for lab problem #3         */
/*              for lab problem #3                     */
/*******************************************************/
/*******************************************************/
</syntaxhighlight>
</syntaxhighlight>
Zeile 49: Zeile 46:


<!-------------------------------------------------------------------------------->
<!-------------------------------------------------------------------------------->
{{MyCodeBlock|title=Declarations
{{MyCodeBlock|title=Declarations
|text=
|text=
Wir arbeiten mit den selben Parametern und Bezugslängen, wie in [[Gelöste Aufgaben/Kw50|Kw50]].
Wir arbeiten mit den selben Parametern und Bezugslängen, wie in [[Gelöste Aufgaben/Kw50|Kw50]].
Insbesondere gilt auch hier wieder
::<math>\displaystyle W_B = -\frac{W_A}{\sqrt{3}}</math>.
|code=
|code=
<syntaxhighlight lang="lisp" line start=1>
<syntaxhighlight lang="lisp" line start=1>
/* declare variational variables - see 6.3 Identifiers */
/* declare variational variables - see 6.3 Identifiers */
declare("Δs", alphabetic);
declare("Δs", alphabetic);
declare( "ϕ", alphabetic); /* = dw/dx*/
declare( "ϕ", alphabetic);  
declare( "Π", alphabetic); /* elastic potential */
declare( "Π", alphabetic); /* elastic potential */
declare( "ℓ", alphabetic);
declare( "ℓ", alphabetic);
declare( "λ", alphabetic);
declare( "Λ", alphabetic);


assume(ℓ[0]>0);
assume(ℓ[0]>0);
Zeile 78: Zeile 68:
geometry: [alpha[A] = 30*%pi/180,
geometry: [alpha[A] = 30*%pi/180,
           alpha[B] = 60*%pi/180,
           alpha[B] = 60*%pi/180,
           ℓ[0]    = ℓ[1]+ℓ[2],
           ℓ[1]    = ℓ[0]-ℓ[2],
           Δs[A]    = W[A]*sin(alpha[A]),
           Δs[A]    = W[A]*cos(%pi/2-alpha[A]),
           Δs[B]    = W[B]*sin(alpha[B]),
           Δs[B]    = W[B]*cos(%pi/2-alpha[B]),
           tan(alpha[B]) = H/ℓ[2],
           tan(alpha[B])= H/ℓ[2],
           tan(alpha[A]) = H/ℓ[0],
           tan(alpha[A])= H/ℓ[0],
           xi[1]    = ℓ[1]/ℓ[0],
           xi[1]    = ℓ[1]/ℓ[0],
           xi[2]    = ℓ[2]/ℓ[0]];
           xi[2]    = ℓ[2]/ℓ[0],
geometry: ratsimp(solve(geometry,[alpha[A],alpha[B],ℓ[1],ℓ[2],Δs[A],Δs[B],H,xi[1],xi[2]])[1]);
          0        = Δs[A]+Δs[B]];
geometry: ratsimp(solve(geometry,[alpha[A],alpha[B],ℓ[1],ℓ[2],Δs[A],Δs[B],H,xi[1],xi[2], W[B]])[1]);


/* reference length selected:                        */
/* reference length selected:                        */
Zeile 93: Zeile 84:


<!-------------------------------------------------------------------------------->
<!-------------------------------------------------------------------------------->
{{MyCodeBlock|title=Formfunctions
{{MyCodeBlock|title=Formfunctions
|text=
|text=
Nach "Ritz" wählen wir Trial-Functions über die gesamte Stablänge und müssen uns überlegen, welchen Aufwand wir dafür treiben wollen.
Nach "Ritz" wählen wir zwei Trial-Functions über die gesamte Stablänge. Als zugehörige, gesuchte Koordinaten wählen wir die Auslenkung in ''A'' und die Verdrehung in ''C'', also
 
Intuitiv wählen wir  für jeden Punkt ''A, B, C'' jeweils eine Koordinate, hier
 
::<math>\begin{array}{l}w(0) = W_A\\w(\ell_1) = W_B\\\displaystyle \frac{dw}{dx}|_{x=\ell} = \Phi_C\end{array}</math>.


und brauchen - zusammen der geometrischen Randbedingung ''w()=0'' ein Polynom mit vier freien Parametern - also ein Polynom dritten Grades.
::<math>\underline{Q} = \left(\begin{array}{c}W_A\\\Phi_C \end{array}\right)</math>.


Mit dem Ansatz für die Formfunktion
Die Trial-Functions müssen dann diesen geoemtrischen Zwangsbedingungen genügen:


::<math>\displaystyle \tilde{w}(x) = \sum_{i=0}^3 C_i\cdot x^i</math>
::<math>\begin{array}{l}w(0) = W_A\\w(\ell_1) = W_B \mbox{ mit } \displaystyle W_B = -\frac{W_A}{\sqrt{3}}\\\displaystyle \frac{dw}{dx}|_{x=\ell} = \Phi_C\end{array}</math>.


kommt aus den Bedingungen oben dann
Mit einem Polynom 3-ten Grades als Ansatz ist also unser Näherungsansatz für die Auslenkung dann


::<math>\displaystyle \tilde{w}( \xi) = \sum_{i=1}^3 Q_i \cdot \phi_i(\xi)</math>
::<math>\displaystyle \tilde{w}(\xi) = \sum_{i=1}^2 Q_i \cdot \phi_i(\xi)</math>


mit
mit den Trial-Functions


::<math>\underline{Q} = \left(\begin{array}{c}W_A\\W_B\\\Phi_C \end{array}\right)</math>.[[Datei:Kw51.png|mini|Trial-Functions]]Und so sehen sie aus, unsere drei [[Sources/Lexikon/Trial-Function|Trial-Functions]]:
::<math>\begin{array}{lll}\phi_1 = &-&\displaystyle \frac{1}{2}\left(\left( {{3}^{\frac{5}{2}}}+3\right) \, {{\xi}^{3}}+\left( -2\; {{3}^{\frac{5}{2}}}-8\right) \, {{\xi}^{2}}+\left( {{3}^{\frac{5}{2}}}+7\right)  \xi-2\right)\\\phi_2 =&&\ell_o\cdot \left( 3\, {{\xi}^{3}}-5 \, {{\xi}^{2}}+2\, \xi\right)\end{array}</math>.


Für die Formfunktion gilt aber immer die geometrische Zwangsbedingung:
[[Datei:Kw53-11.png|mini|Trial Functions]]Und so sehen sie aus, unsere zwei [[Sources/Lexikon/Trial-Function|Trial-Functions]]:
 
::<math>\displaystyle W_B = -\frac{W_A}{\sqrt{3}}</math>.
|code=
|code=
<syntaxhighlight lang="lisp" line start=1>
<syntaxhighlight lang="lisp" line start=1>
/* coordinates                                       */
/* coordinates                                       */
Q : [W[A],W[B],Phi[C]];
Q : [W[A],Phi[C]];
/* raw trial function                                */
v(x) := sum(C[i]*x^i,i,0,3);
v(x) := sum(C[i]*x^i,i,0,3);
 
/* gemetric constraints */
/* formfunctions                                      */
const: [subst([x= 0  ],      v(x)  ) = W[A],
const: [subst([x= 0  ],      v(x)  ) = W[A],
         subst([x=ℓ[1]],      v(x)  ) = W[B],
         subst([x=ℓ[1]],      v(x)  ) = W[B],
Zeile 132: Zeile 114:
         subst([x=ℓ[0]],      v(x)  ) =  0  ];
         subst([x=ℓ[0]],      v(x)  ) =  0  ];
trials : expand(subst(solve(subst(geometry,const), makelist(C[i],i,0,3))[1],v(x)));
trials : expand(subst(solve(subst(geometry,const), makelist(C[i],i,0,3))[1],v(x)));
phi : makelist(ratsimp(coeff(trials,Q[i])),i,1,3);
phi : makelist(ratsimp(coeff(trials,Q[i])),i,1,2);
 
trials: w(x) = sum(Q[i]*phi[i],i,1,2);
print('phi = expand(subst([x=ℓ[0]*xi],phi)))$
 
/* plot trial functions                              */
plot2d(ratsimp(subst([x = xi*ℓ[0]],phi)*[1,1/ℓ[0]]), [xi,0,1],
            [legend, "W[A]","Φ[C]"], [xlabel, "x/ℓ →"], [ylabel, "ϕ[i] →"]);
</syntaxhighlight>
</syntaxhighlight>
}}
}}


<!-------------------------------------------------------------------------------->
<!-------------------------------------------------------------------------------->
{{MyCodeBlock|title=Equilibrium Conditions
 
{{MyCodeBlock|title=Potentials
|text=
|text=
Die Potentiale aus Elastischer Energien und Arbeitsfunktion der Gewichtskräfte sind
Die Potentiale aus Elastischen Energien und Arbeitsfunktion der Gewichtskräfte sind


::<math>\begin{array}{lll} U  =& &\displaystyle \frac{1}{2}\cdot \int_0^\ell E I\; w''^2 \; dx + \frac{1}{2}\cdot K_C \cdot \Phi_C^2\\    &-&\displaystyle \int_0^\ell q_0 \; w \; dx - m_A\,g\; W_A \end{array}</math>.
::<math>\begin{array}{lll} U  =& &\displaystyle \frac{1}{2}\cdot \int_0^\ell E I\; w''^2 \; dx + \frac{1}{2}\cdot K_C \cdot \Phi_C^2\\    &-&\displaystyle \int_0^\ell q_0 \; w \; dx - m_A\,g\; W_A \end{array}</math>.


Nach dem Prinzip vom Minimum der Potentiellen Energie ist das System im Gleichgewicht, wenn das Potential ''U'' des Systems
Einsetzten der Trial-Functions liefert
 
::<math>\displaystyle \frac{dU}{dQ_i} = 0 \text{ für } Q_i \in \left( W_A, W_B, \Phi_C\right)</math>


erfüllt. Und sich zusätzlich an die Zwangsbedingung hält!
::<math>U= \displaystyle \frac{1}{2} \left({W_A},{{\Phi}_C}\right) \underbrace{\begin{pmatrix}\frac{\left( 8 {{3}^{\frac{5}{2}}}+262\right) \, \mathit{EI}}{{{\ell}_{0}^{3}}} & -\frac{\left( 5 {{3}^{\frac{5}{2}}}+17\right) \, \mathit{EI}}{{{\ell}_{0}^{2}}}\\ -\frac{\left( 5 {{3}^{\frac{5}{2}}}+17\right) \, \mathit{EI}}{{{\ell}_{0}^{2}}} & \frac{28 \mathit{EI}+{\ell_0}\, {K_C}}{{\ell_0}}\end{pmatrix}}_{\displaystyle :=\underline{\underline{A}}} \begin{pmatrix}{W_A}\\ {{\Phi}_C}\end{pmatrix} - \left({W_A},{{\Phi}_C}\right) \underbrace{\begin{pmatrix}{m_A} g-\frac{{{3}^{\frac{3}{2}}}\, {q_0}\, {\ell_0}}{8}+\frac{5 {q_0}\, {\ell_0}}{24}\\ \frac{{q_0}\, {{\ell}_{0}^{2}}}{12}\end{pmatrix}}_{\displaystyle :=\underline{b}}</math>.
|code=
|code=
<syntaxhighlight lang="lisp" line start=1>
<syntaxhighlight lang="lisp" line start=1>
Zeile 156: Zeile 144:
PMPE : [Π[P] = 1/2*integrate(EI*'diff(w(x),x,2)^2, x,0,ℓ[0]) + 1/2*K[C]*Phi[C]^2,
PMPE : [Π[P] = 1/2*integrate(EI*'diff(w(x),x,2)^2, x,0,ℓ[0]) + 1/2*K[C]*Phi[C]^2,
         A[P] =  integrate(q[0]*w(x), x,0,ℓ[0]) + m[A]*g*W[A]];
         A[P] =  integrate(q[0]*w(x), x,0,ℓ[0]) + m[A]*g*W[A]];
PMPE: subst(trials, PMPE);
PMPE: ev(PMPE,nouns);
/* Potential                                          */
U: expand(subst(PMPE,Π[P] - A[P]));
</syntaxhighlight>
</syntaxhighlight>
}}
}}


<!-------------------------------------------------------------------------------->
<!-------------------------------------------------------------------------------->
{{MyCodeBlock|title=Geometric Constraints
 
{{MyCodeBlock|title=Equilibrium Conditions
|text=
|text=
Das können wir mit dem Konzept der Lagrange-Multiplikatoren erfassen. Dafür schreiben wir
Nach dem [[Sources/Lexikon/Minimum Prinzipe|Minimum Prinzip]] hat das Potential ''U'' ein Minimum - das Sytem ist im Gleichgewicht, wenn
 
::<math>\Lambda (\underline{Q})  = U + \lambda\cdot \underbrace{\left( \Delta s_A + \Delta s_B\right)}_{\displaystyle \equiv \sqrt{3}\, W_B+W_A} </math>


mit dem Lagrange-Multiplikator ''λ''. Die neuen Gleichgewichtsbedingungen lauten nun
::<math>\underline{\underline{A}} \cdot \underline{Q} = \underline{b}</math>.


::<math>\displaystyle \frac{d\Lambda}{dQ_i} = 0 \text{ für alle } Q_i </math>
In den System-Matrizen stehen hier


und wir erhalten die vier Gleichungen
::<math>\underline{\underline{A}}\,=\,\begin{pmatrix}\displaystyle \frac{386.7\cdot \mathit{EI}}{{{\ell}_{0}^{3}}} & \displaystyle -\frac{94.9\cdot \mathit{EI}}{{{\ell}_{0}^{2}}}\\ \displaystyle -\frac{94.9\cdot \mathit{EI}}{{{\ell}_{0}^{2}}} & \displaystyle \frac{{{\ell}_{0}}\cdot {{K}_{C}}+28\cdot \mathit{EI}}{{{\ell}_{0}}}\end{pmatrix}
</math> und


::<math>\begin{array}{cc} \displaystyle \frac{\lambda }{2}-\frac{49 {m_B} g}{120}-\frac{17 {{\Phi}_C}\, EI}{{\ell_{0}^{2}}}-\frac{108 {W_B}\, EI}{{\ell_{0}^{3}}}+\frac{19 {W_A}\, EI}{{\ell_{0}^{3}}}&=0\\ \displaystyle \frac{\sqrt{3} \lambda }{2}-\frac{9 {m_B} g}{8}+\frac{135 {{\Phi}_C}\, EI}{{\ell_{0}^{2}}}+\frac{729 {W_B}\, EI}{{\ell_{0}^{3}}}-\frac{108 {W_A}\, EI}{{\ell_{0}^{3}}}&=0\\
::<math>\underline{b}\,=\,\begin{pmatrix}{{m}_{A}}\,g-0.4\cdot {{q}_{0}}\cdot {{\ell}_{0}}\\ 0.08\cdot {{q}_{0}}\cdot {{\ell}_{0}^{2}}\end{pmatrix}</math>.
\displaystyle -\frac{{\ell_0}\, {m_B} g}{12}+\frac{33 {{\Phi}_C}\, EI}{{\ell_0}}+\frac{135 {W_B}\, EI}{{\ell_{0}^{2}}}-\frac{17 {W_A}\, EI}{{\ell_{0}^{2}}}&=0\\
\displaystyle \frac{\sqrt{3}\, {W_B}}{2}+\frac{{W_A}}{2}&=0 \end{array}</math>.
|code=
|code=
<syntaxhighlight lang="lisp" line start=1>
<syntaxhighlight lang="lisp" line start=1>
/* gemetric constraints                              */
trials: w(x) = sum(Q[i]*phi[i],i,1,3);
print('phi = expand(subst([x=ℓ[0]*xi],phi)))$
plot2d(ratsimp(subst([x = xi*ℓ[0]],phi)*[1,1,1/ℓ[0]]), [xi,0,1],
            [legend, "W[A]","W[B]","Φ[C]"], [xlabel, "x/ℓ →"], [ylabel, "ϕ[i] →"]);
PMPE: subst(trials, PMPE);
PMPE: ev(PMPE,nouns);
U: expand(subst(PMPE,Π[P] - A[P]));
/* Lagrange-Function                                  */
Λ : U + subst(geometry,λ*(Δs[A]+Δs[B]));
/* Equilibrium Conditions                            */
/* Equilibrium Conditions                            */
eom : subst(params,makelist( diff(Λ,Q[i])=0, i,1,4));
A[M] : ratsimp(funmake('matrix,2*[
                    [coeff(U,Q[1],2), coeff(coeff(U,Q[1],1),Q[2],1)/2],
                    [coeff(coeff(U,Q[2],1),Q[1],1)/2, coeff(U,Q[2],2)]
                    ]));
rest : expand(U-1/2*Q.A[M].transpose(Q));
b[M] : funmake('matrix,-makelist([coeff(rest,Q[i])],i,1,2));
</syntaxhighlight>
</syntaxhighlight>
}}
}}


<!-------------------------------------------------------------------------------->
<!-------------------------------------------------------------------------------->
{{MyCodeBlock|title=Solving
{{MyCodeBlock|title=Solving
|text=
|text=
Das Lösen des Gleichungssystems liefert dann
Das Lösen des Gleichungssystems liefert dann


::<math>\displaystyle \left(\begin{array}{c} W_A\\W_B\\\Phi_C\\\lambda \end{array}\right) = \frac{m_B\, g\, \ell_0^3}{3 \; E I}  \left(\begin{array}{l}-3.78 10^{-5}\\+2.18 10^{-5}\\\displaystyle +0.00747 \frac{1}{\ell_0}\\\displaystyle +2.71 \frac{EI}{\ell_0^3} \end{array} \right)</math>.
::<math>\displaystyle \left(\begin{array}{c} W_A\\\Phi_C\end{array}\right) = \frac{m_B\, g\, \ell_0^3}{3 \; E I}  \left(\begin{array}{l} -3.78 {10}^{-5}\\\displaystyle 0.00747 \frac{1}{\ell} \end{array} \right)</math>.
|code=
|code=
<syntaxhighlight lang="lisp" line start=1>
<syntaxhighlight lang="lisp" line start=1>
/* Solving                                            */
/* Solving                                            */
sol: float(solve(eom,Q)[1]);
sol: linsolve_by_lu(subst(params,A[M]),subst(params,b[M]))[1];
sol : expand(float(sol));
sol : makelist(Q[i]=sol[i][1],i,1,2);
</syntaxhighlight>
</syntaxhighlight>
}}
}}
Zeile 215: Zeile 200:
Und die Ergebnisse können wir uns anschauen ...
Und die Ergebnisse können wir uns anschauen ...
==== ... für w(x): ====
==== ... für w(x): ====
[[Datei:Kw51-11.png|mini|Auslenkung ''w(x)''|alternativtext=|ohne]]
[[Datei:Kw53-21.png|mini|Auslenkung ''w(x)''|alternativtext=|ohne]]


<!-- Die Ergebnis-Plots sind identisch mit denen aus Kw51 -->
==== ... für ''Φ(x)'': ====
==== ... für ''Φ(x)'': ====
[[Datei:Kw51-12.png|mini|Kippwinkel ''Φ(x)''|alternativtext=|ohne]]
[[Datei:Kw51-12.png|mini|Kippwinkel ''Φ(x)''|alternativtext=|ohne]]
Zeile 228: Zeile 214:
<syntaxhighlight lang="lisp" line start=1>
<syntaxhighlight lang="lisp" line start=1>
/* Post-Processing                                    */
/* Post-Processing                                    */
w : subst([x=xi*ℓ[0]],subst(geometry,subst(sol, sum(Q[j]*phi[j],j,1,3))));
w : subst([x=xi*ℓ[0]],subst(geometry,subst(sol, sum(Q[j]*phi[j],j,1,2))));


fcts: [        w            ,
fcts: [        w            ,

Aktuelle Version vom 31. März 2021, 06:38 Uhr


Aufgabenstellung

Eine Brücke ABC der Masse mB und homogener Biegesteifigkeit EI ist in C gelenkig gelagert und in A sowie B mit einem Seil verbunden. Das undehnbare Seil wird dabei über eine kleine Rolle (Radius r ≪ ℓ) in D haftungsfrei geführt. In Punkt C ist die Brücke über eine Drehfeder der Steifigkeit KC mit dem Lager verbunden. In A steht eine Person der Masse mA.


Lageplan (wie Kw50)

Geben Sie die Lösung für ein Euler-Bernoulli-Modell der Brücke mit dem Verfahren von Rayleigh-Ritz (EBB) an - hier ohneLagrange-Multiplikator für die geometrische Zwangsbedingung.

Dies ist eine Näherungslösung zu Kw50.

Ermitteln Sie die genäherten Verläufe der Schnittgrößen und Verschiebungen im Balken für diese Parameter:

Lösung mit Maxima

In dieser Aufgabe berechnen wir eine Näherungslösung nach dem Verfahren von Rayleigh-Ritz (EBB) zu Kw50.

Alle Überlegungen zur Geometrie des Systems übernehmen wir.

Header

Für die Lösung nutzen wir hier das Ritz-Verfahren. Die geometrischen Zwangsbedingungen arbeiten wir direkt in die Trial-Functions ein.


/*******************************************************/
/* MAXIMA script                                       */
/* version: wxMaxima 18.10.1                           */
/* author: Andreas Baumgart                            */
/* last updated: 2019-02-12                            */
/* ref: TM-C, Brigde-Problem                           */
/* description: finds the plain rayleigh-ritz solution */
/*              for lab problem #3                     */
/*******************************************************/




Declarations

Wir arbeiten mit den selben Parametern und Bezugslängen, wie in Kw50.


/* declare variational variables - see 6.3 Identifiers */
declare("Δs", alphabetic);
declare( "ϕ", alphabetic); 
declare( "Π", alphabetic); /* elastic potential */
declare( "ℓ", alphabetic);

assume(ℓ[0]>0);

/* system parameters                                  */
params: [K[C]     = kappa*EI/ℓ[0],
         q[0]     = m[B]*g/ℓ[0],
         m[A]     = theta*m[B],
         theta    = 1/5,
         kappa    = 5];

geometry: [alpha[A] = 30*%pi/180,
           alpha[B] = 60*%pi/180,
           ℓ[1]     = ℓ[0]-ℓ[2],
           Δs[A]    = W[A]*cos(%pi/2-alpha[A]),
           Δs[B]    = W[B]*cos(%pi/2-alpha[B]),
           tan(alpha[B])= H/ℓ[2],
           tan(alpha[A])= H/ℓ[0],
           xi[1]    = ℓ[1]/ℓ[0],
           xi[2]    = ℓ[2]/ℓ[0],
           0        = Δs[A]+Δs[B]];
geometry: ratsimp(solve(geometry,[alpha[A],alpha[B],ℓ[1],ℓ[2],Δs[A],Δs[B],H,xi[1],xi[2], W[B]])[1]);

/* reference length selected:                         */
dimless : ℓ[Bez] = 1/3*m[B]*g*ℓ[0]^3/(EI); /*cantilevered*/




Formfunctions

Nach "Ritz" wählen wir zwei Trial-Functions über die gesamte Stablänge. Als zugehörige, gesuchte Koordinaten wählen wir die Auslenkung in A und die Verdrehung in C, also

.

Die Trial-Functions müssen dann diesen geoemtrischen Zwangsbedingungen genügen:

.

Mit einem Polynom 3-ten Grades als Ansatz ist also unser Näherungsansatz für die Auslenkung dann

mit den Trial-Functions

.
Trial Functions

Und so sehen sie aus, unsere zwei Trial-Functions:


/* coordinates                                       */
Q : [W[A],Phi[C]];
v(x) := sum(C[i]*x^i,i,0,3);
/* gemetric constraints */
const: [subst([x= 0  ],      v(x)  ) = W[A],
        subst([x=ℓ[1]],      v(x)  ) = W[B],
        subst([x=ℓ[0]], diff(v(x),x))= Phi[C],
        subst([x=ℓ[0]],      v(x)  ) =  0  ];
trials : expand(subst(solve(subst(geometry,const), makelist(C[i],i,0,3))[1],v(x)));
phi : makelist(ratsimp(coeff(trials,Q[i])),i,1,2);

trials: w(x) = sum(Q[i]*phi[i],i,1,2);
print('phi = expand(subst([x=ℓ[0]*xi],phi)))$

/* plot trial functions                               */
plot2d(ratsimp(subst([x = xi*ℓ[0]],phi)*[1,1/ℓ[0]]), [xi,0,1],
            [legend, "W[A]","Φ[C]"], [xlabel, "x/ℓ →"], [ylabel, "ϕ[i] →"]);




Potentials

Die Potentiale aus Elastischen Energien und Arbeitsfunktion der Gewichtskräfte sind

.

Einsetzten der Trial-Functions liefert

.

/******************************************************/
/* Boundary Value Problem Formulation                 */
/* elastic and gravitational potential                */

PMPE : [Π[P] = 1/2*integrate(EI*'diff(w(x),x,2)^2, x,0,ℓ[0]) + 1/2*K[C]*Phi[C]^2,
        A[P] =  integrate(q[0]*w(x), x,0,ℓ[0]) + m[A]*g*W[A]];

PMPE: subst(trials, PMPE);
PMPE: ev(PMPE,nouns);

/* Potential                                          */
U: expand(subst(PMPE,Π[P] - A[P]));




Equilibrium Conditions

Nach dem Minimum Prinzip hat das Potential U ein Minimum - das Sytem ist im Gleichgewicht, wenn

.

In den System-Matrizen stehen hier

und
.

/* Equilibrium Conditions                             */
A[M] : ratsimp(funmake('matrix,2*[
                    [coeff(U,Q[1],2), coeff(coeff(U,Q[1],1),Q[2],1)/2],
                    [coeff(coeff(U,Q[2],1),Q[1],1)/2, coeff(U,Q[2],2)]
                    ]));
rest : expand(U-1/2*Q.A[M].transpose(Q));
b[M] : funmake('matrix,-makelist([coeff(rest,Q[i])],i,1,2));




Solving

Das Lösen des Gleichungssystems liefert dann

.

/* Solving                                            */
sol: linsolve_by_lu(subst(params,A[M]),subst(params,b[M]))[1];
sol : expand(float(sol));
sol : makelist(Q[i]=sol[i][1],i,1,2);




Post-Processing

Und die Ergebnisse können wir uns anschauen ...

... für w(x):

Auslenkung w(x)

... für Φ(x):

Kippwinkel Φ(x)

... für M(x):

Biegemoment M(x)

... für Q(x):

Querkraft Q(x)

/* Post-Processing                                    */
w : subst([x=xi*ℓ[0]],subst(geometry,subst(sol, sum(Q[j]*phi[j],j,1,2))));

fcts: [         w             ,
           diff(w,xi  )/ℓ[0]  ,
       -EI*diff(w,xi,2)/ℓ[0]^2,
       -EI*diff(w,xi,3)/ℓ[0]^3];
fcts: float(subst(geometry,expand(fcts)))$ 
facts: [1/ℓ[Bez], ℓ[0]/ℓ[Bez], 1/(m[B]*g*ℓ[0]), 1/(m[B]*g)];
 
textlabels : ["← w(x)/ℓ[Bez]", "← w'(x)/(ℓ[Bez]/ℓ[0]) →", "M(x)/(m[B]*g*ℓ) →", "Q(x)/(m[B]g →"];
for i: 1 thru 4 do(
  f : expand(subst(dimless,facts[i]*fcts[i])),
  preamble: if i<=2 then "set yrange [] reverse" else "set yrange []",
  plot2d(f, [xi,0,1], [legend, false],
                      [gnuplot_preamble, preamble],
                      [xlabel, "x/ℓ →"],
                      [ylabel, textlabels[i]]))$





Links

  • Aufgabe Kw50 (analytische Lösung dieser Aufgabe)
  • Aufgabe Kw52 (Lösung dieser Aufgabe mit dem Ansatz von Rayleigh-Ritz und Lagrange-Multiplikator)
  • Aufgabe Kw53 (Lösung dieser Aufgabe mit dem Ansatz von Rayleigh-Ritz)

Literature

  • ...