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

  ${\underset{\_}{\underset{\_}{K}}}_i=\frac{E{I}_i}{{\ell }_i^3}\left(\begin{array}{cccc}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{array}\right)$

• und die Element-"Rechte Seite"

  ${\underset{\_}{P}}_i=\rho A_{i}g{\ell }_i\cdot \left(\begin{array}{c}\frac{1}{2}\\ \frac{{\ell }_i}{12}\\ \frac{1}{2}\\ -\frac{{\ell }_i}{12}\end{array}\right)$.

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

${\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

${\tilde{A}}_1=\frac{(\frac{13\cdot h_0}{7}+\frac{h_1}{7})\cdot b}{2},{\tilde{A}}_2=\frac{(\frac{10\cdot h_0}{7}+\frac{4\cdot h_1}{7})\cdot b}{2},{\tilde{A}}_3=\frac{(\frac{4\cdot h_0}{7}+\frac{10\cdot h_1}{7})\cdot b}{2}$

und

${\tilde{I}}_1=\frac{{(\frac{13\cdot h_0}{7}+\frac{h_1}{7})}^3\cdot b}{96},{\tilde{I}}_2=\frac{{(\frac{10\cdot h_0}{7}+\frac{4\cdot h_1}{7})}^3\cdot b}{96},{\tilde{I}}_3=\frac{{(\frac{4\cdot h_0}{7}+\frac{10\cdot h_1}{7})}^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 \left(\begin{array}{cccccc}\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{array}\right)\cdot \underset{\_}{Q}=mg\cdot \left(\begin{array}{c}\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{array}\right)$.

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

$\hat{W}=\frac{q_{ref}{\ell }_0^4}{8E{I}_{ref}}$

und

$q_{ref}=\frac{mg}{{\ell }_0},E{I}_{ref}=Eb\frac{{\left(\frac{h_0+h_1}{2}\right)}^3}{12}$

zu

$\begin{array}{ccl}\frac{W_0}{\hat{W}}& =& 0,\\ \frac{{\mathrm{\Phi }}_0}{\hat{W}}& =& 0,\\ \frac{W_1}{\hat{W}}& =& 0.016,\\ \frac{{\mathrm{\Phi }}_1}{\hat{W}}& =& \frac{0.2093}{{\ell }_0},\\ \frac{W_2}{\hat{W}}& =& 0.1317,\\ \frac{{\mathrm{\Phi }}_2}{\hat{W}}& =& \frac{0.5545}{{\ell }_0},\\ \frac{W_3}{\hat{W}}& =& 0.5937,\\ \frac{{\mathrm{\Phi }}_3}{\hat{W}}& =& \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

• ...