Gelöste Aufgaben/MaMa

Aufgabenstellung

Ein Stab der Länge $3a$ und Biegesteifigkeit $EI$ ist links fest eingespannt und wird bei 2/3 seiner Länger durch ein gelenkiges Lager gestützt. Zwischen dem freien rechten Rand und dem gelenkigen Lager wirkt eine konstante Streckenlast $q_{0}$.

Lageplan

Gesucht ist die analytische Lösung für den Euler-Bernoulli-Balken und die Verläufe der Schnittgrößen. Im Vordergrund der Ausarbeitung steht das Zusammenspiel von Maxima (Computeralgebra) und Matlab (Numerik) zum Aufstellen, Verwalten und Lösen des linearen Gleichungssystems.

Die Aufgabe ist einfach genug, um sie analytisch mit Maxima allein zu lösen. Hier geht es aber darum, die "Zusammenarbeit" zwischen den beiden Software-Paketen exemplarisch vorzustellen. Dieses Zusammenspiel wird dann interessant, wenn die Anzahl der Gleichungen groß wird - typischer weise > 50. Dann ist eine algebraische Lösung oft nicht mehr möglich. Maxima eignet sich dann hervorragend für das Aufstellen und Verwalten der Gleichungen. Die konsolidierten Gleichungen werden dann zur numerischen Lösung an Matlab übergeben. Wie das geht, zeigen wir im Folgenden.

Lösung mit Maxima und Matlab

Im folgenden finden Sie die Skripte für Maxima und Matlab, mit denen diese Aufgabe bearbeitet wurde.

Declarations

Wir zeichnen zunächst die Berieche, Knoten und Koordinaten aus:

Bezeichnung für Bereiche und Knoten einführen.

Für Bereich $I$ führen wir die unabhängige Koordinate $x_{I}$ ein, für Bereich $II$ die Koordinate $x_{II}$. Die zwei Bereich haben die Längen

• $\ell _{I}$ und
• $\ell _{II}$.

Um die gesuchte Lösung dimensionslos ansetzten zu können, wählen wir eine Vergleichslösung. Hier bietet sich der Kragbalken unter Streckenlast an. Diese liefert für die maximale Auslenkung von $w(x)$ und die maximlae Kippung der Querschnitte von $\phi (x)$

• $w_{ref}={\frac {q_{0}a^{4}}{8EI}}$ und
• $\varphi _{ref}={\frac {q_{0}a^{3}}{6EI}}$

Außerdem wählen wir

• $M_{ref}={\frac {q_{0}a^{2}}{2}}$ und
• $Q_{ref}=q_{0}a$

Boundary Value Problem

Das Randwertproblem wird im Feld von der Euler-Bernoulli-Differentialgleichung des Balkens beschreiben und an den Rändern durch kinematische oder Kraft-Randbedingungen.

Differentialgleichung des Euler-Bernoulli-Balkens

Die Gleichgewichtsbedingung am Balken führt auf die Differentialgleichung

$EIw^{IV}=q_{0}$

Deren allgemeine Lösung (also ohne Anpassung an Rand- und Übergangsbedingungen) ist

... für die Auslenkung des Querschnitts:	$EIw(x)={\frac {1}{24}}q_{0}x^{4}+{\frac {1}{6}}C_{3}x^{3}+{\frac {1}{2}}C_{2}x^{2}+C_{1}x+C_{0}$

mit den 4 Integrationskonstanten - und damit unseren Unbekannten je Abschnitt -

$C_{0},C_{1},C_{2}{\text{ und }}C_{3}$.

Aus Gleichgewichts- und Kinematik-Beziehungen kommen dann die weiteren Gleichungen

... für den Kippwinkel des Querschnitts:	$EI\varphi (x)=EI{\frac {dw}{dx}}=+{\frac {1}{6}}q_{0}x^{3}+{\frac {1}{2}}C_{3}x^{2}+C_{2}x+C_{1}$
... für das Biegemoment:	$M(x)=-EI{\frac {d^{2}w}{dx^{2}}}=-{\frac {1}{2}}q_{0}x^{2}-C_{3}x-C_{2}$ und
... für die Querkraft:	$Q(x)={\frac {dM}{dx}}=-q_{0}x-C_{3}$ und

In unserer Aufgabe haben wir zwei Bereiche - in B müssen wir diese allgemeinen Lösungen aneinander anstückeln. Um sie zu unterscheiden, hängen wir für jede Lösung den Index des Bereichs an, also

Bereich I $0<x_{I}<\ell _{I}$	Bereich II $0<x_{II}<\ell _{II}$
$EIw_{I}(x)={\frac {1}{24}}q_{I,0}x_{I}^{4}+{\frac {1}{6}}C_{I,3}x_{I}^{3}+{\frac {1}{2}}C_{I,2}x_{I}^{2}+C_{I,1}x+C_{I,0}$	$EIw_{II}(x)={\frac {1}{24}}q_{II,0}x_{I}^{4}+{\frac {1}{6}}C_{II,3}x_{I}^{3}+{\frac {1}{2}}C_{II,2}x_{I}^{2}+C_{II,1}x+C_{II,0}$

Und das gilt in gleicher Weise für die Gleichungen von $\varphi ,M,Q$.

Rand- und Übergangsbedingungen

Die Rand- und Übergangsbedingungen im Stab finden wir entweder aus kinematischen Überlegungen oder aus Gleichgewichtsbeziehungen mit den Schnittkräften des Balkens.

Kinematische Rand- und Übergangsbedingungen

Für die Kinematik des Systems schauen wir uns die angenommene Verformung des Systems an:

Angenommene Verformung des Stabes.

Aus dem Bild entnehmen wir, dass

• ... in A:
  • $w_{I}(0)=0$
  • $\varphi _{I}(0)=0$
• ... in B:

  

  Kippwinkel bei "B".

  • $w_{I}(\ell _{I})=0$
  • $w_{II}(0)=0$ und
  • $\varphi _{I}(\ell _{I})=\varphi _{II}(0)=\Phi _{B}$

gilt.

Kraft- und Momenten Rand- und Übergangsbedingungen

Für die Übergangsbedingungen in B brauchen wir ein Freikörperbild:

Schnittbild in B.

Die Gleichgewichtsbeziehungen erhalten wir dann zu

• ... in B:
  • $w_{I}(0)=0$ und
• ... in C:
  • $Q_{II}(\ell _{I}I)=0$
  • $M_{II}(\ell _{I}I)=0$

Dies sind insgesamt 8 Rand- bzw. Übergangsbedingungen für die 8 Integrationskonstanten

$C_{I,0},C_{I,1},C_{I,2},C_{I,3},C_{II,0},C_{II,1},C_{II,2}{\text{ und }}C_{II,3}$.

Konsolidierung der Gleichungen

Jetzt kommt Maxima ins Spiel.

Maxima-Quellcode

Im Quellcode für Maxima werden die Elemente aus dem Lösungsprozess oben direkt umgesetzt. So finden sich in "fcts" die Funktionen für $w,\varphi ,M,Q$.

Diese Funktionen können wir mit Maxima direkt in den Randbeindungen ansprechen und erhalten die 5+3=8 Gleichungen

${\begin{array}{lccl}1)&C_{I,0}/EI&=&0\\2)&C_{I,1}/EI&=&0\\3)&4C_{I,3}a^{3}/(3EI)+2C_{I,2}a^{2}/EI+2C{I,1}a/EI+C_{I,0}/EI&=&0\\4)&2C_{1,3}a^{2}/EI+2C_{I,2}a/EI+C_{1,1}/EI&=&C_{II,1}/EI\\5)&-2C_{I,3}a-C_{I,2}&=&-C_{II,2}\\6)&C_{II,0}/EI&=&0\\7)&-q_{0}a^{2}/2-C_{II,3}a-C_{II,2}&=&0\\8)&-q_{0}a-C_{II,3}&=&0\end{array}}$

Diese Gleichungen kann man einzeln mit passenden Faktoren durchmultiplizieren und man erhält das einfacher aussehende lineare Gleichungssystem

${\begin{pmatrix}1&0&0&0&0&0&0&0\\0&1&0&0&0&0&0&0\\1&2a&2{{a}^{2}}&{\frac {4{{a}^{3}}}{3}}&0&0&0&0\\0&1&2a&2{{a}^{2}}&0&-1&0&0\\0&0&-1&-\left(2a\right)&0&0&1&0\\0&0&0&0&1&0&0&0\\0&0&0&0&0&0&-1&-a\\0&0&0&0&0&0&0&-1\end{pmatrix}}\cdot {\begin{pmatrix}{C_{I,0}}\\{C_{I,1}}\\{C_{I,2}}\\{C_{I,3}}\\{C_{II,0}}\\{C_{II,1}}\\{C_{II,2}}\\{C_{II,3}}\end{pmatrix}}={\begin{pmatrix}0\\0\\0\\0\\0\\0\\{\frac {{q_{0}}{{a}^{2}}}{2}}\\{q_{0}}a\end{pmatrix}}$

Maxima hilft auch dabei, die von Null verschiedenen Elemente herauszusortieren und für Matlab zu formatieren für ${\underline {\underline {A}}}$:

   A( 1 , 1 ) =  1;
   A( 2 , 2 ) =  1; 
   A( 3 , 1 ) =  1;
   A( 3 , 2 ) =  2*a; 
   A( 3 , 3 ) =  2*a^2; 
   A( 3 , 4 ) =  (4*a^3)/3; 
   A( 4 , 2 ) =  1; 
   A( 4 , 3 ) =  2*a; 
   A( 4 , 4 ) =  2*a^2; 
   A( 4 , 6 ) =  -1; 
   A( 5 , 3 ) =  -1; 
   A( 5 , 4 ) =  -(2*a); 
   A( 5 , 7 ) =  1; 
   A( 6 , 5 ) =  1; 
   A( 7 , 7 ) =  -1; 
   A( 7 , 8 ) =  -a; 
   A( 8 , 8 ) =  -1; 

und ${\underline {b}}$:

   b(7,1) = q0*a^2/2;
   b(8,1) = q0*a;

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/* Modul Numerische Methoden der Mechanik                            */
/* Script zum Unterricht                                             */
/* Löst: Gross e.a., Technische Mechanik, Bd 2, A 3.16               */
/* Autor: Andreas Baumgart                                           */
/* Last Updated: 2024-11-13                                          */
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/

/* parameters */
params: [ℓ[1]=2*a, ℓ[2]=a, q[1,0]=0];

/* reference-values */
refs: [w[ref]=q[2,0]*a^4/8/EI,
       φ[ref]=q[2,0]*a^3/6/EI,
       M[ref]=q[2,0]*a^2/2,
       Q[ref]=q[2,0]*a];

/*****************************************/
/* the differential equation for the EBB */
ebb: [EI*diff(w[i](x),x,4)=0, w[i](x) = 1/EI*(q[i,0]*x^4/24+sum(1/j!*C[i,j]*x^j,j,0,3))];
/* ... and its solution */
fcts: [w[i](x) =    diff(subst(ebb,w[i](x)),x,0),
       φ[i](x) =    diff(subst(ebb,w[i](x)),x,1),
       M[i](x) =-EI*diff(subst(ebb,w[i](x)),x,2),
       Q[i](x) =-EI*diff(subst(ebb,w[i](x)),x,3)];

/* boundary conditions */
bvp: [subst([i=1,x=  0 ], subst(fcts,w[i](x))) = 0,
      subst([i=1,x=  0 ], subst(fcts,φ[i](x))) = 0,
      subst([i=1,x=ℓ[1]], subst(fcts,w[i](x))) = 0,
      subst([i=1,x=ℓ[1]], subst(fcts,φ[i](x))) = subst([i=2,x=  0 ], subst(fcts,φ[i](x))),
      subst([i=1,x=ℓ[1]], subst(fcts,M[i](x))) = subst([i=2,x=  0 ], subst(fcts,M[i](x))),
      subst([i=2,x=  0 ], subst(fcts,w[i](x))) = 0,
      subst([i=2,x=ℓ[2]], subst(fcts,M[i](x))) = 0,
      subst([i=2,x=ℓ[2]], subst(fcts,Q[i](x))) = 0];
      
X: flatten(makelist(makelist(C[i,j],j,0,3),i,1,2));

/* scale equations - so they look nicer ... */
scale: [EI,EI,EI,EI,1,EI,1,1];
for i:1 thru length(bvp) do
    bvp[i]: subst(params,bvp[i])*scale[i]; 

/* preprocess */
bvp: expand(subst(refs,subst(params,bvp)));
ACM: augcoefmatrix(bvp,X);
A: submatrix(ACM,9);
b:-col(ACM,9);

print(A,"*",transpose(X),"=",b)$

/* you'll need to transfere these coefficients to Matlab: */
for row:1 thru 8 do
    for col:1 thru 8 do
        if notequal(subst(a=1,A[row][col]), 0) then
            print("A(",row,",",col,") = ", A[row][col])$

/* Output for Matlab .....
   A( 1 , 1 ) =  1;
   A( 2 , 2 ) =  1; 
   A( 3 , 1 ) =  1;
   A( 3 , 2 ) =  2*obj.a; 
   A( 3 , 3 ) =  2*obj.a^2; 
   A( 3 , 4 ) =  (4*obj.a^3)/3; 
   A( 4 , 2 ) =  1; 
   A( 4 , 3 ) =  2*obj.a; 
   A( 4 , 4 ) =  2*obj.a^2; 
   A( 4 , 6 ) =  -1; 
   A( 5 , 3 ) =  -1; 
   A( 5 , 4 ) =  -(2*obj.a); 
   A( 5 , 7 ) =  1; 
   A( 6 , 5 ) =  1; 
   A( 7 , 7 ) =  -1; 
   A( 7 , 8 ) =  -obj.a; 
   A( 8 , 8 ) =  -1; 
*/

/* solve                           */
/* we could also solve in Maxima : */
sol: solve(bvp,X)[1];

/* postprocess - an example: */
toPlot: expand(subst([x=ξ*a],subst(refs,subst(params,subst(sol,makelist(subst([i=j],subst(fcts, w[i](x)/w[ref])),j,1,2))))));
/* plot dimensionless displacements of beams in secsions I und II */
print("reference displacement w[ref] = ", subst(refs,w[ref]))$
plot2d([[parametric,ξ,toPlot[1],[ξ,0,2]],[parametric,2+ξ,toPlot[2],[ξ,0,1]]],[xlabel,"x/a->"],[ylabel,"w/w[ref]->"], [legend, "w_1", "w_2"]);

/* these are the functions needed in postprocessing in Matlab : */
for k:1 thru 4 do
    (print(" for ", lhs(fcts[k])/lhs(refs[k]),":"),
     print(expand(subst([x=ξ*a],subst(params,makelist(subst([i=j],subst(refs,subst(fcts, lhs(fcts[k])/lhs(refs[k])))),j,1,2))))))$
/* note that x=xi*a has been introduced as the dimensionless coordinate */

Matlab^©-Quellcode

Genau so wird der Code dann in Matlab eingefügt - hier mit der Änderung, dass a und q0 Properties der EBB-Class sind und entsprechend innerhalb der Klasse als obj.a bzw obj.q angesprochen werden müssen, laso

        :
        :
        function A = sysmat(obj)
            %populate
            %   system matrix A

            % create matrix
            A   = zeros(8,8);
            % populate matrix (only non-zero elements)
            A( 1 , 1 ) =  1;
            A( 2 , 2 ) =  1; 
            A( 3 , 1 ) =  1;
            A( 3 , 2 ) =  2*obj.a; 
            A( 3 , 3 ) =  2*obj.a^2; 
            A( 3 , 4 ) =  (4*obj.a^3)/3; 
            A( 4 , 2 ) =  1; 
            A( 4 , 3 ) =  2*obj.a; 
            A( 4 , 4 ) =  2*obj.a^2; 
            A( 4 , 6 ) =  -1; 
            A( 5 , 3 ) =  -1; 
            A( 5 , 4 ) =  -(2*obj.a); 
            A( 5 , 7 ) =  1; 
            A( 6 , 5 ) =  1; 
            A( 7 , 7 ) =  -1; 
            A( 7 , 8 ) =  -obj.a; 
            A( 8 , 8 ) =  -1; 
        end
                
        function b = rhs(obj)
            %populate
            %   right-hand side

            % cretae matrix
            b   = zeros(8,1);
            % populate matrix
            b(7,1) = obj.q*obj.a^2/2;
            b(8,1) = obj.q*obj.a;
        end
        :
        :

Folder Structure

Die Datei-Struktur links zeigt das Skript MaMa.m sowie Classes und Functions in den jeweiligen Ordnern. Die Excel-Datei hält alle System-Parameter.

Den komplette Quellcode zu diesem Programm können Sie über dieses ZIP-File rechts herunterladen.

archive

download compressed archive →

Solving

Die Gleichungen könnten wir in Maxima mit dem Befehl

solve(bvp,X)[1]

lösen - das Ergebnis ist dann

$C_{I,0}=0,C_{I,1}=0,C_{I,2}=-q_{0}a^{2}/4,C_{I,3}=3q_{0}a/8,C_{II,0}=0,C_{II,1}=q_{0}a^{3}/4,C_{II,2}=q_{0}a^{2}/2,C_{II,3}=-q_{0}a$

Das gleiche passiert in der Matlab-Funktion linsolve in der solveAB.m:

function x = solveAb(sys)
%   solve 
%   linear system of equations A*x = b

A = sys.sysmat();
b = sys.rhs();

x = linsolve(A,b);
end

Wir erhalten die Lösung

X =

   1.0e+03 *

         0
         0
   -0.5000
    0.7500
         0
    0.5000
    1.0000
   -2.0000

${\displaystyle }$ ${\displaystyle }$

Post-Processing

Mit der algebraischen Lösung tragen wir $w(x)/w_{ref}$ auf und finden

Dimensionslose Auslenkung aus der algebraischen Lösung mit Maxima.

Die gleiche Lösung- und zusätzliche Graphen - kommen aus der numerischen Lösung mit Matlab:

Dimensionslose Auslenkung aus der numerischen Lösung mit Maxima.

Die Lösungen stimmen - natürlich - überein. In Fällen, bei denen das lineare Gleichungssystem so wenige Unbekannte wie hier hat, würde man also die Lösung direkt in Maxima ermitteln. Wenn es komplizierter wird, muss man oft in die Numerik gehen, um in endlicher Zeit eine Lösung zu bekommen.

---

Links

• ...

Literature

• ...