Gelöste Aufgaben/UEBL

Aufgabenstellung

Die Problemstellung ist identisch zu UEBI (und UEBJ, UEBK):

Der Euler-Bernoulli-Balken AB wird durch seine Gewichtskraft belastet. Er ist in A fest eingespannt und hat eine konstante Breite b sowie eine zwischen A und B linear veränderliche Höhe h.

Lageplan

Gesucht ist eine Näherungslösung mit der Methode der Finiten Elementen, die nur Elemente mit stückweise konstanten Querschnitts-Eigenschaften zulassen.

Gegeben sind für den Balken:

• Länge ℓ, Breite b,
• E-Modul E, Dichte ρ und
• die Höhe h₀=b und h₁ jeweils in A und B; dazwischen ist die Höhe linear veränderlich.

Lösung mit Maxima

Header

Unsere FEM-Formulierung mit kubischen Ansatz-Polynomen für Euler-Bernoulli-Balken ist auf Stäbe mit elementweise konstantem Querschnitt beschränkt: E und I müssen je Element konstant sein.

Für dieses Problem können wir die Formulierung also - nominell - nicht anwenden. Wir müssten im Abschnitt Virtuelle Formänderungsenergie die Integration für linear-veränderliche Querschnittsflächen ansetzen und durchführen.

Das machen wir hier nicht.

FE-Modell

Wir schauen uns an, wie die Lösung für das Ersatz-Modell aussieht, bei dem wir uns den Originalstab durch einen stufenweise abgesetzten Stab ersetzt denken - hier für zwei Elemente.

Die Querschnittseigenschaften je Element seien die des Original-Stabes im Element-Mittelpunkt.

/*******************************************************/
/* MAXIMA script                                       */
/* version: wxMaxima 15.08.2                           */
/* author: Andreas Baumgart                            */
/* last updated: 2019-09-01                            */
/* ref: TMC                                            */
/* description: EBB with I(x), A(x)                    */
/*              solved as Initial Value Problem        */
/*******************************************************/

Declarations

Wir benutzen die Parameter-Deklarationen aus UEBI und laden außerdem die analytische Lösung von dort.

Einziger Unterschied: um die Element-Länge ℓ_i in Element "i" eindeutig von der Gesamt-Länge des Stabes zu unterscheiden, wählen wir für die Gesamt-Länge den Parameter ℓ₀ statt ℓ.

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

/* system parameters                                        */
params: [q(xi) = A(xi)*rho*g,
         A(xi) = b*h(xi),
         I(xi) = b*h(xi)^3/12,
         h(xi) = H[0]*(1-xi)+ H[1]*xi];
params: append(params,
               solve((H[0]+H[1])/2*b*ℓ[0]*rho=m, rho));

geometry : [alpha=1/2];

dimless: [x = xi*ℓ, H[0]=b, H[1]=alpha*b];

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

/* load analytic solutions as discrete tables from UEBI */
file_search_maxima: append(file_search_maxima, ["C:/Users/XXXX/OneDrive/Confluence Sources/UEBI/"]);
batch("analyticSolFromUEBI.data")$

FEM-Formulation

Wir verwenden aus Abschnitt FEM-Formulierung für den Euler-Bernoulli-Balken

• die Standard Element-Steifigkeitsmatrix

  $\displaystyle {\underline {\underline {K}}}_{i}={\frac {EI_{i}}{\ell _{i}^{3}}}\;{\begin{pmatrix}12&6\cdot {{\ell }_{i}}&-12&6\cdot {{\ell }_{i}}\\6\cdot {{\ell }_{i}}&4\cdot {{\ell }_{i}^{2}}&-6\cdot {{\ell }_{i}}&2\cdot {{\ell }_{i}^{2}}\\-12&-6\cdot {{\ell }_{i}}&12&-6\cdot {{\ell }_{i}}\\6\cdot {{\ell }_{i}}&2\cdot {{\ell }_{i}^{2}}&-6\cdot {{\ell }_{i}}&4\cdot {{\ell }_{i}^{2}}\end{pmatrix}}$

• und die Element-"Rechte Seite"

  ${\underline {P}}_{i}=\rho \,A_{i}\,g\,\ell _{i}\cdot {\begin{pmatrix}\displaystyle {\frac {1}{2}}\\\displaystyle {\frac {{\ell }_{i}}{12}}\\\displaystyle {\frac {1}{2}}\\-\displaystyle {\frac {{\ell }_{i}}{12}}\end{pmatrix}}$.

Hier dürfen also der (mittlere) Querschnittsflächeninhalt A_i, das (mittlere) Flächenmoment I_i und die Länge ℓ_i je Element unterschiedliche Werte annehmen.

Ab hier haben wir einige Freiheiten und die füllen wir folgendermaßen:

• Die Anzahl der Finiten Element ist 3.
• Die Element-Längen wählen wir so, dass sie sich in Stab-Längsrichtung verdoppeln, also $\ell _{i+1}=2\cdot \ell _{i}$

Dann sind mit

$\lambda _{i}:=\ell _{i}/\ell _{0}$

die relativen Element-Längen

$\displaystyle {{\lambda }_{1}}={\frac {1}{7}},{{\lambda }_{2}}={\frac {2}{7}},{{\lambda }_{3}}={\frac {4}{7}}$.

Für diese drei Element bestimmen wir die "mid-point"-Querschnittseigenschaften zu

$\displaystyle {\tilde {A}}_{1}={\frac {\left({\frac {13\cdot {{h}_{0}}}{7}}+{\frac {{h}_{1}}{7}}\right)\cdot b}{2}},{\tilde {A}}_{2}={\frac {\left({\frac {10\cdot {{h}_{0}}}{7}}+{\frac {4\cdot {{h}_{1}}}{7}}\right)\cdot b}{2}},{\tilde {A}}_{3}={\frac {\left({\frac {4\cdot {{h}_{0}}}{7}}+{\frac {10\cdot {{h}_{1}}}{7}}\right)\cdot b}{2}}$

und

$\displaystyle {\tilde {I}}_{1}={\frac {{{\left({\frac {13\cdot {{h}_{0}}}{7}}+{\frac {{h}_{1}}{7}}\right)}^{3}}\cdot b}{96}},{\tilde {I}}_{2}={\frac {{{\left({\frac {10\cdot {{h}_{0}}}{7}}+{\frac {4\cdot {{h}_{1}}}{7}}\right)}^{3}}\cdot b}{96}},{\tilde {I}}_{3}={\frac {{{\left({\frac {4\cdot {{h}_{0}}}{7}}+{\frac {10\cdot {{h}_{1}}}{7}}\right)}^{3}}\cdot b}{96}}$

Einsetzen und komponieren der Gesamt-Systemmatrix - mit Einarbeiten der geometrischen Randbedingungen - liefert

${\frac {{{b}^{4}}\cdot E}{{\ell }_{0}^{3}}}\cdot {\begin{pmatrix}{\frac {21411}{64}}&-{\frac {16227\cdot {{\ell }_{0}}}{896}}&-27&{\frac {27\cdot {{\ell }_{0}}}{7}}&0&0\\-{\frac {16227\cdot {{\ell }_{0}}}{896}}&{\frac {8865\cdot {{\ell }_{0}^{2}}}{3136}}&-{\frac {27\cdot {{\ell }_{0}}}{7}}&{\frac {18\cdot {{\ell }_{0}^{2}}}{49}}&0&0\\-27&-{\frac {27\cdot {{\ell }_{0}}}{7}}&{\frac {14553}{512}}&-{\frac {6183\cdot {{\ell }_{0}}}{1792}}&-{\frac {729}{512}}&{\frac {729\cdot {{\ell }_{0}}}{1792}}\\{\frac {27\cdot {{\ell }_{0}}}{7}}&{\frac {18\cdot {{\ell }_{0}^{2}}}{49}}&-{\frac {6183\cdot {{\ell }_{0}}}{1792}}&{\frac {1395\cdot {{\ell }_{0}^{2}}}{1568}}&-{\frac {729\cdot {{\ell }_{0}}}{1792}}&{\frac {243\cdot {{\ell }_{0}^{2}}}{3136}}\\0&0&-{\frac {729}{512}}&-{\frac {729\cdot {{\ell }_{0}}}{1792}}&{\frac {729}{512}}&-{\frac {729\cdot {{\ell }_{0}}}{1792}}\\0&0&{\frac {729\cdot {{\ell }_{0}}}{1792}}&{\frac {243\cdot {{\ell }_{0}^{2}}}{3136}}&-{\frac {729\cdot {{\ell }_{0}}}{1792}}&{\frac {243\cdot {{\ell }_{0}^{2}}}{1568}}\end{pmatrix}}\cdot {\underline {Q}}=m\,g\cdot {\begin{pmatrix}{\frac {25}{98}}\\{\frac {23\cdot {{\ell }_{0}}}{4116}}\\{\frac {20}{49}}\\{\frac {16\cdot {{\ell }_{0}}}{1029}}\\{\frac {12}{49}}\\-{\frac {8\cdot {{\ell }_{0}}}{343}}\end{pmatrix}}$.

/*****************************************************************/
/* FEM-Formulation                                               */
NOE : 3;     /* number of elements */
NON : NOE+1; /* number of nodes    */

/* element matrices */
K[i] : (E*I[i]/ℓ[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[i]*g*ℓ[i]*matrix([   1/2  ],
                              [ ℓ[i]/12],
                              [   1/2  ],
                              [-ℓ[i]/12]);

/* element lengths                    */
lengths : solve(append([sum(λ[i],i,1,NOE)=1], makelist(λ[i]=2*λ[i-1],i,2,NOE)), makelist(λ[i],i,1,NOE))[1];
xCoord  : subst(lengths, makelist(sum(λ[i],i,1,j),j,0,NOE));

/* element cross sectional properties      */
height : makelist(subst([xi=xCoord[i]],subst(params,h(xi))),i,1,NON);

/** discrete proeprties at element mid-point-sections  */
Amp : makelist(  (height[i]+height[i+1])/2       * b, i,1,NOE);
Imp : makelist( ((height[i]+height[i+1])/2)^3/12 * b, i,1,NOE);

/* define all element matrices  */
for e:1 thru NOE do (
  K[e] : ratsimp(subst(geometry,subst(lengths,subst(dimless,subst([ℓ[e]=ℓ[0]*λ[e], A[e]=Amp[e], I[e]=Imp[e]], subst([i=e], K[i])))))),
  P[e] : ratsimp(subst(geometry,subst(lengths,subst(dimless,subst([ℓ[e]=ℓ[0]*λ[e], A[e]=Amp[e], I[e]=Imp[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])));

/* incorporate bounday conditions */

AA : submatrix(1,2,K[0],1,2);
bb : submatrix(1,2,P[0]);

Solving

Das lineare Gleichungssystem lösen wir mit der Referenz-Auslenkung

$\displaystyle {\hat {W}}={\frac {q_{ref}\,\ell _{0}^{4}}{8\,E\,I_{ref}}}$

und

$\displaystyle q_{ref}={\frac {m\,g}{\ell _{0}}},E\,I_{ref}=E\,b\,{\frac {\left(\displaystyle {\frac {h_{0}+h_{1}}{2}}\right)^{3}}{12}}$

zu

${\begin{array}{ccl}\displaystyle {\frac {W_{0}}{\hat {W}}}&=&0,\\\displaystyle {\frac {\Phi _{0}}{\hat {W}}}&=&0,\\\displaystyle {\frac {W_{1}}{\hat {W}}}&=&0.016,\\\displaystyle {\frac {\Phi _{1}}{\hat {W}}}&=&\displaystyle {\frac {0.2093}{\ell _{0}}},\\\displaystyle {\frac {W_{2}}{\hat {W}}}&=&0.1317,\\\displaystyle {\frac {\Phi _{2}}{\hat {W}}}&=&\displaystyle {\frac {0.5545}{\ell _{0}}},\\\displaystyle {\frac {W_{3}}{\hat {W}}}&=&0.5937,\\\displaystyle {\frac {\Phi _{3}}{\hat {W}}}&=&\displaystyle {\frac {0.8932}{\ell _{0}}}\end{array}}$

mit den Bezeichnungen nach Abschnitt FEM: Trial-Functions für kubische Ansatz-Polynome.

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

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

Post-Processing

Die Ergebnisse schauen wir uns als dimensionsloser Darstellung an, wobei wir wie in UEBI die Standard-Lösungen für den Balken unter konstanter Streckenlast als Referenz-Lösung ansetzen.

Wir finden

• ... für die Auslenkung w:
  

  

  Biegelinie w(x)

• ... für den Kippwinkel ϕ:
  

  

  Kippwinkel ϕ(x)

• ... für das Biegemoment M:
  

  

  Biegemomentenverlauf M(x)

• ... für die Querkraft Q:
  

  

  Querkraftverlauf 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) ->"];

/* diagrams 
analyticSolution : [WA,PhiA,MA,QA]$ /* from file .... */                                                              */
for diag:1 thru 4 do (
    /* pick a reference value to plot dimensionless */
    Ref : [W[ref], ϕ[ref], M[ref], Q[ref]][diag],
    toPlot : ratsimp(makelist([parametric, xCoord[e]*(1-t)+xCoord[e+1]*t,
        subst(geometry,subst(dimless,subst(reference,subst([xi=t],subst(lengths,subst(nodalSol,subst([i=e],subst([ℓ[i]=λ[i]*ℓ[0]],
                        /* factor is 1 or -EI - depending on property to plot */
                        [1,1,-E*Imp[e],-E*Imp[e]][diag]*
                        sum( [W[i-1],Phi[i-1],W[i],Phi[i]][j]*diff(phi[j],xi,diag-1)/ℓ[i]^(diag-1),j,1,4)
                        )))/Ref))))), [t,0,1]],e,1,NOE)),
    /* also plot the analytical solution */
    toPlot: append([[discrete, analyticSolution[diag]]],toPlot),
    preamble: if diag<=2 then "set yrange [] reverse" else "set yrange []",
    plot2d(toPlot, [legend, "analytic", "Element 1", "Element 2", "Element 3", "Element 4", "Element 5"],
                   [gnuplot_preamble, preamble],
                   [xlabel, "x/ℓ ->"],
                   [ylabel, textlabels[diag]]) )$

---

Links

• UEBI (analytischer Lösung), UEBJ, UEBK

Literature

• ...