Gelöste Aufgaben/Kw60

Aufgabenstellung

Der Euler-Bernoulli-Balken ABC (Länge ℓ₀, Biegesteifigkeikeit EI, Querschnittsflläche A) ist am rechten Rand fest eingespannt und wird durch das Gewicht

• aufgrund seiner eigenen Masse m₀ sowie
• aufgrund von zwei Personen der Massen m_A bzw. m_B

belastet.

Lageplan

Gesucht ist eine Näherungslösung mit der Methode der Finiten Elemente und zwei Elementen. Dabei sollen die Verläufe der Biegelinie und die Schnittlasten mit der analytischen Lösung verlglichen werden.

Es ist ℓ_B =  ℓ₀ /3, m_A =  m₀ /2 und m_B =  m₀ /2.

Lösung mit Maxima

Header

Das Modell erstellen wir mit der FEM-Formulierung für den Euler-Bernoulli-Balken mit zwei Elementen.

/*******************************************************/
/* MAXIMA script                                       */
/* version: wxMaxima 15.08.2                           */
/* author: Andreas Baumgart                            */
/* last updated: 2019-09-01                            */
/* ref: TMC                                            */
/* description: simple EBB with two section            */
/*              solved as FE Problem                   */
/*******************************************************/

Declarations

Zum späteren Plotten verwenden wir Referenzgrößen, die wir der analytischen Standard-Lösung eines Kragbalkens entnehmen.

Dies sind

$\begin{array}{ccl}{W}_{{r}{e}{f}}& =& \frac{q_0\cdot {\ell }_0^4}{8\cdot I_0\cdot E},\\ {\mathrm{\Phi }}_{{r}{e}{f}}& =& \frac{q_0\cdot {\ell }_0^3}{6\cdot I_0\cdot E},\\ M_{{r}{e}{f}}& =& \frac{m_0\cdot {\ell }_0\cdot g}{2},\\ Q_{{r}{e}{f}}& =& m_0\cdot g,\\ q_0& =& \frac{m_0\cdot g}{{\ell }_0}\end{array}$.

Wir verwenden dimensionslose Größen als Abkürzungen:

${\ell }_1={\lambda }_1{\ell }_0,{\ell }_B={\lambda }_2{\ell }_0,{\lambda }_1=1-{\lambda }_2,m_A={\theta }_A{m}_0,m_B={\theta }_B{m}_0$.

/************************************************************/
declare( "ℓ", alphabetic);
declare( "λ", alphabetic);
declare( "ϕ", alphabetic);

/* system parameters                                        */
params: [q[0] = A[0]*rho*g];
params: append(params,
               solve(A[0]*ℓ[0]*rho=m[0], rho));

dimless: [ℓ[1]=λ[1]*ℓ[0],ℓ[2]=λ[2]*ℓ[0],λ[1]=1-λ[2], m[A]=theta[A]*m[0], m[B]=theta[B]*m[0]];

geometry: [λ[2] = 1/3, theta[A]=1/2, theta[B]=1/2];

/* make equations of motion dim'less with load case #1 from Gross e.a., same as UEBI */
reference : [W[ref] = q[0]*ℓ[0]^4/(8*E*I[0]), ϕ[ref] = q[0]*ℓ[0]^3/(6*E*I[0]), M[ref] = m[0]*g*ℓ[0]/2, Q[ref] = m[0]*g, q[0] = m[0]*g/ℓ[0]];

Analytic Solution

Die analytische Lösung erzeugen wir aus Superpositionen der eingeprägten Lösungen zu den drei Gewichten auf die Struktur.

Dazu müssen wir ein bisschen umdenken: in den Standard-Lösungen ist der Balken jeweils links fest eingespannt - hier rechts.

Wir nehmen das komplementäre Koordinatensystem

$\overline{\xi }=1-\xi$

und lesen die Lösung ab:

$\begin{array}{lccl}EIw(\overline{\xi })& =& & \frac{q_0{\ell }_0^4\cdot (6\cdot {\overline{\xi }}^2-4\cdot {\overline{\xi }}^3+{\overline{\xi }}^4)}{24}\\ & & +& \frac{m_{A}g{\ell }_0^3\cdot (3\cdot {\overline{\xi }}^2-{\overline{\xi }}^3)}{6}\\ & & +& \{\begin{array}{ll}\frac{m_{B}g{\ell }_0^3\cdot (3\cdot {\lambda }_2\cdot {\overline{\xi }}^2-{\overline{\xi }}^3)}{6}& \text{ für }\overline{\xi }\le {\lambda }_2\\ \frac{m_{B}g{\ell }_0^3\cdot (3\cdot {\lambda }_2\cdot {\overline{\xi }}^2-{\overline{\xi }}^3+{(\overline{\xi }-{\lambda }_2)}^3)}{6}& \text{ sonst}\end{array}\end{array}$.

/**************************************************************************/
/* analytic solution                                                      */
/* ... from superposition of analytical solutions for cantileverd beam    */
analySol : [q[0]*ℓ[0]^4/(24*E*I[0])*(6*xi^2-4*xi^3+xi^4),
            [m[B]*g*ℓ[0]^3/(6*E*I[0])*(3*xi^2*λ[2]-xi^3 +     0),
             m[B]*g*ℓ[0]^3/(6*E*I[0])*(3*xi^2*λ[2]-xi^3 + (xi-λ[2])^3)],
            m[A]*g*ℓ[0]^3/(6*E*I[0])*(3*xi^2* 1 -xi^3)];
/* obs: zwo solutions for load case 2 (weight in "B")                     */

FE Formulation

Die System-Matrix setzen wir aus den Element-Steifigkeitsmatrizen aus FEM: Trial-Functions für kubische Ansatz-Polynome zusammen und arbeiten die geometrischen Randbedingungen ein, indem wir die Zeilen und Spalten bzgl der Verschiebung und Verdrehung in C (W₂, Φ₂) streichen.

Das Gleichungssystem lautet dann

$\frac{EI}{{\ell }_0^3}\cdot \left(\begin{array}{cccc}-\frac{12}{{\lambda }_2^3-3\cdot {\lambda }_2^2+3\cdot {\lambda }_2-1}& \frac{6\cdot {\ell }_0}{{\lambda }_2^2-2\cdot {\lambda }_2+1}& \frac{12}{{\lambda }_2^3-3\cdot {\lambda }_2^2+3\cdot {\lambda }_2-1}& \frac{6\cdot {\ell }_0}{{\lambda }_2^2-2\cdot {\lambda }_2+1}\\ \frac{6\cdot {\ell }_0}{{\lambda }_2^2-2\cdot {\lambda }_2+1}& -\frac{4\cdot {\ell }_0^2}{{\lambda }_2-1}& -\frac{6\cdot {\ell }_0}{{\lambda }_2^2-2\cdot {\lambda }_2+1}& -\frac{2\cdot {\ell }_0^2}{{\lambda }_2-1}\\ \frac{12}{{\lambda }_2^3-3\cdot {\lambda }_2^2+3\cdot {\lambda }_2-1}& -\frac{6\cdot {\ell }_0}{{\lambda }_2^2-2\cdot {\lambda }_2+1}& -\frac{12-36\cdot {\lambda }_2+36\cdot {\lambda }_2^2}{{\lambda }_2^6-3\cdot {\lambda }_2^5+3\cdot {\lambda }_2^4-{\lambda }_2^3}& -\frac{12\cdot {\ell }_0\cdot {\lambda }_2-6\cdot {\ell }_0}{{\lambda }_2^4-2\cdot {\lambda }_2^3+{\lambda }_2^2}\\ \frac{6\cdot {\ell }_0}{{\lambda }_2^2-2\cdot {\lambda }_2+1}& -\frac{2\cdot {\ell }_0^2}{{\lambda }_2-1}& -\frac{12\cdot {\ell }_0\cdot {\lambda }_2-6\cdot {\ell }_0}{{\lambda }_2^4-2\cdot {\lambda }_2^3+{\lambda }_2^2}& -\frac{4\cdot {\ell }_0^2}{{\lambda }_2^2-{\lambda }_2}\end{array}\right)\cdot \underset{\_}{Q}=m_{0}g\cdot \left(\begin{array}{c}\frac{1-{\lambda }_2+2\cdot {\theta }_A}{2}\\ \frac{{\ell }_0-2\cdot {\ell }_0\cdot {\lambda }_2+{\ell }_0\cdot {\lambda }_2^2}{12}\\ \frac{1+2\cdot {\theta }_B}{2}\\ \frac{2\cdot {\ell }_0\cdot {\lambda }_2-{\ell }_0}{12}\end{array}\right)$

mit

$\underset{\_}{Q}=\left(\begin{array}{c}{W}_0\\ {\mathrm{\Phi }}_0\\ W_1\\ {\mathrm{\Phi }}_1\\ \end{array}\right)$.

/**************************************************************************/
/* FE Formulation                                                         */
/* element matrices */
K[i] : (E*I[0]/ℓ[i]^3)*matrix([12,6*ℓ[i],-12,6*ℓ[i]],
                              [6*ℓ[i],4*ℓ[i]^2,-6*ℓ[i],2*ℓ[i]^2],
                              [-12,-6*ℓ[i],12,-6*ℓ[i]],
                              [6*ℓ[i],2*ℓ[i]^2,-6*ℓ[i],4*ℓ[i]^2]);

P[i] : rho*A[0]*g*ℓ[i]*matrix([   1/2  ],
                              [ ℓ[i]/12],
                              [   1/2  ],
                              [-ℓ[i]/12]);

/* use one element per section                                             */
NOE : 2;     /* number of elements */
NON : NOE+1; /* number of nodes    */

xCoord: subst(geometry,[0,1-λ[2],1]);

/* define both element matrices  */
for e:1 thru NOE do (
  K[e] : ratsimp(subst(dimless,subst([ℓ[e]=ℓ[0]*λ[e]], subst([i=e], K[i])))),
  P[e] : ratsimp(subst(dimless,subst([ℓ[e]=ℓ[0]*λ[e]], subst([i=e], P[i])))));

/* compose system matrix */
K[0] : zeromatrix(2*NON,2*NON);
P[0] : zeromatrix(2*NON,1);

for ele:1 thru NOE do (
   for row:1 thru 4 do (
      P[0][2*(ele-1)+row,1] : P[0][2*(ele-1)+row,1] + P[ele][row,1],
      for col:1 thru 4 do (
         K[0][2*(ele-1)+row,2*(ele-1)+col] : K[0][2*(ele-1)+row,2*(ele-1)+col] + K[ele][row,col])));

P[0] : ratsimp(subst(params,P[0]));

/* add discrete loads an points A and B */
P[0][1,1] : P[0][1,1]+m[A]*g;
P[0][3,1] : P[0][3,1]+m[B]*g;


/* incorporate bounday conditions */
AA : submatrix(5,6,K[0],5,6);
bb : submatrix(5,6,P[0]);

Solving

Lösen des linearen Gleichungssystems liefert

$\begin{array}{ccl}\frac{W_0}{W_{ref}}& =& -\frac{-3-8\cdot {\theta }_A+(4\cdot {\lambda }_2^3-12\cdot {\lambda }_2^2)\cdot {\theta }_B}{3}\\ \frac{{\mathrm{\Phi }}_0}{W_{ref}}& =& -\frac{4+12\cdot {\theta }_A+12\cdot {\lambda }_2^2\cdot {\theta }_B}{3\cdot {\ell }_0}\\ \frac{W_1}{W_{ref}}& =& \frac{6\cdot {\lambda }_2^2-4\cdot {\lambda }_2^3+{\lambda }_2^4+(12\cdot {\lambda }_2^2-4\cdot {\lambda }_2^3)\cdot {\theta }_A+8\cdot {\lambda }_2^3\cdot {\theta }_B}{3}\\ \frac{{\mathrm{\Phi }}_1}{W_{ref}}& =& -\frac{12\cdot {\lambda }_2-12\cdot {\lambda }_2^2+4\cdot {\lambda }_2^3+(24\cdot {\lambda }_2-12\cdot {\lambda }_2^2)\cdot {\theta }_A+12\cdot {\lambda }_2^2\cdot {\theta }_B}{3\cdot {\ell }_0}\\ \frac{W_2}{W_{ref}}& =& 0\\ \frac{{\mathrm{\Phi }}_2}{W_{ref}}& =& 0\end{array}$.

/***************************************************************************/
/* solve                                                                   */
sol: ratsimp(subst(dimless,linsolve_by_lu(AA,bb)))[1];

nodalSol: ratsimp(subst(dimless,subst(reference,flatten(append(
       makelist(
         [ W [i-1] = sol[2*i-1][1],
          Phi[i-1] = sol[2*i  ][1]],i,1,NOE),
         [ W [  2] =  0 ,
          Phi[  2] =  0 ])))));

Post-Processing

Und die Ergebnisse der analytischen und FE-Rechnungen plotten wir ...

... für w(x):

Biegelinie w(x)

... für Φ(x):

Kippung w'(x)

... für M(x):

Biegemoment M(x)

... für Q(x):

Querkraft Q(x)

/***************************************************************************/
/* plot results                                                            */

/* Trial-Fucntions                                                         */
phi : [       (xi-1)^2*(2*xi+1),
         ℓ[i]* xi     *(  xi-1)^2,
        -      xi^2   *(2*xi-3),
         ℓ[i]* xi^2   *(  xi-1)];

textlabels : ["<- w(x)/w[ref]", "<- w'(x)/ϕ[ref]", "M(x)/(m*g*ℓ/2) ->", "Q(x)/(m g) ->"];

for D:0 thru 3 do (
    /* pick a reference value to plot dimensionless */
    Ref : [W[ref], ϕ[ref], M[ref], Q[ref]][D+1],
    Fac : [   1  ,    1  ,-E*I[0],-E*I[0]][D+1],
    toPlotFEM : makelist([parametric,
                      xCoord[e]*(1-t)+xCoord[e+1]*t,
                      expand(subst([xi=t],subst(geometry,subst(reference,subst(nodalSol,subst(dimless, subst([i=e],
                              Fac*diff(sum( [W[i-1],Phi[i-1],W[i],Phi[i]][j]*phi[j],j,1,4),xi,D)/ℓ[i]^D/Ref
                        ))))))),
                      [t,0,1]],e,1,NOE),
    toPlotAly : subst(geometry,subst([xi=t],
             [[parametric, 1-xi,(-1)^D*ratsimp(Fac*diff(subst(geometry,subst(reference,subst(dimless,(analySol[1]+analySol[2][1]+analySol[3])/Ref))),xi,D)/ℓ[0]^D), [t,0,λ[2]]],
              [parametric, 1-xi,(-1)^D*ratsimp(Fac*diff(subst(geometry,subst(reference,subst(dimless,(analySol[1]+analySol[2][2]+analySol[3])/Ref))),xi,D)/ℓ[0]^D), [t,λ[2],1]]])),
    toPlot: append(toPlotAly,toPlotFEM),
    preamble: if D<=1 then "set yrange [] reverse" else "set yrange []",
    plot2d(toPlot, [legend, "analytic", "analytic", "FEM", "FEM"],
               [gnuplot_preamble, preamble],
               [style, [lines, 2],[lines, 2],[lines, 1],[lines, 1]],
               [xlabel, "x/ℓ ->"],
               [ylabel, textlabels[D+1]]));

✔ Achtung: Koordinatentransformation:

Wegen der Koordinaten-Transformation zwischen ξ und ihrer komplementären Koordinate in der analytischen Lösung müssen wir den Richtungssinn für ϕ und Q umdrehen!

---

Links

• ...

Literature

• ...