Gelöste Aufgaben/UEBL: Unterschied zwischen den Versionen

Aus numpedia
Zur Navigation springen Zur Suche springen
KKeine Bearbeitungszusammenfassung
Keine Bearbeitungszusammenfassung
 
(5 dazwischenliegende Versionen desselben Benutzers werden nicht angezeigt)
Zeile 1: Zeile 1:
j
[[Category:Gelöste Aufgaben]]
[[Datei:UEBL-33.png|mini|Biegemomentenverlauf ''M(x)'']]
[[Category:Numerische Lösung]]
[[Datei:UEBL-34.png|mini|Querkraftverlauf ''Q(x)'']]
[[Category:Randwertproblem]]
[[Datei:UEBL-32.png|mini|Kippwinkel ''ϕ(x)'']]
[[Category:Biege-Belastung]][[Category:Euler-Bernoulli-Balken]]
[[Datei:UEBL-31.png|mini|Biegelinie w(x)]]
[[Category:Finite-Elemente-Methode‎]]
[[Category:Maxima‎]]


==Aufgabenstellung==
Die Problemstellung ist identisch zu [[Gelöste Aufgaben/UEBI|UEBI]] (und [[Gelöste Aufgaben/UEBJ|UEBJ]], [[Gelöste Aufgaben/UEBK|UEBK]]):


[[Datei:UEBL-11.png|mini|FE-Modell]]
Der [[Sources/Lexikon/Euler-Bernoulli-Balken|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''.<onlyinclude>
[[Datei:UEBF-01.png|200px|left|mini|Lageplan]]
Gesucht ist eine Näherungslösung mit der Methode der Finiten Elementen, die nur Elemente mit stückweise konstanten Querschnitts-Eigenschaften zulassen.
</onlyinclude>
Gegeben sind für den Balken:
 
* Länge ''ℓ'', Breite ''b,''
* E-Modul ''E'', Dichte ''ρ'' und
* die Höhe ''h''<sub>0</sub>''=b'' und ''h<sub>1</sub>'' jeweils in ''A'' und ''B''; dazwischen ist die Höhe linear veränderlich.
 
== Lösung mit Maxima ==
<!-------------------------------------------------------------------------------->
{{MyCodeBlock|title=Header
|text=
Unsere [[Randwertprobleme/Methoden zur Lösung von Randwertproblemen/FEM: Trial-Functions für kubische Ansatz-Polynome|FEM-Formulierung mit kubischen Ansatz-Polynomen]] für [[Sources/Lexikon/Euler-Bernoulli-Balken|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 [[Sources/Lexikon/Virtuelle Formänderungsenergie|Virtuelle Formänderungsenergie]] die Integration für linear-veränderliche Querschnittsflächen ansetzen und durchführen.
 
Das machen wir hier nicht.[[Datei:UEBL-11.png|mini|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.
|code=
<syntaxhighlight lang="lisp" line start=1>
/*******************************************************/
/* 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        */
/*******************************************************/
</syntaxhighlight>
}}
 
<!-------------------------------------------------------------------------------->
{{MyCodeBlock|title=Declarations
|text=
Wir benutzen die Parameter-Deklarationen aus UEBI und laden außerdem die analytische Lösung von dort.
 
Einziger Unterschied: um die Element-Länge ''ℓ<sub>i</sub>'' in Element "i" eindeutig von der Gesamt-Länge des Stabes zu unterscheiden, wählen wir für die Gesamt-Länge den Parameter ''ℓ<sub>0</sub>'' statt ''ℓ.''
|code=
<syntaxhighlight lang="lisp" line start=1>
/************************************************************/
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")$
</syntaxhighlight>
}}
 
<!-------------------------------------------------------------------------------->
{{MyCodeBlock|title=FEM-Formulation
|text=
Wir verwenden aus Abschnitt FEM-Formulierung für den Euler-Bernoulli-Balken
 
<ul>
<li>die Standard Element-Steifigkeitsmatrix
::<math>\displaystyle \underline{\underline{K}}_i = \frac{E I_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}</math></li>
<li>und die Element-"Rechte Seite"
::<math>\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}</math>.</li>
</ul>
 
Hier dürfen also der (mittlere) Querschnittsflächeninhalt ''A<sub>i</sub>'', das (mittlere) Flächenmoment ''I<sub>i</sub>'' und die Länge ''ℓ<sub>i</sub>'' je Element unterschiedliche Werte annehmen.
 
Ab hier haben wir einige Freiheiten und die füllen wir folgendermaßen:
 
<ul>
<li>Die Anzahl der Finiten Element ist 3.</li>
<li>Die Element-Längen wählen wir so, dass sie sich in Stab-Längsrichtung verdoppeln, also
<math>\ell_{i+1} = 2 \cdot \ell_i</math></li>
</ul>
 
Dann sind mit
 
::<math>\lambda_i := \ell_i/\ell_0</math>
 
die relativen Element-Längen
 
::<math>\displaystyle {{\lambda}_{1}}=\frac{1}{7},{{\lambda}_{2}}=\frac{2}{7},{{\lambda}_{3}}=\frac{4}{7}</math>.
 
Für diese drei Element bestimmen wir die "mid-point"-Querschnittseigenschaften zu
 
::<math>\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}</math>
 
und
 
::<math>\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}</math>
 
Einsetzen und komponieren der Gesamt-Systemmatrix - mit Einarbeiten der geometrischen Randbedingungen - liefert
 
::<math>\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}</math>.
|code=
<syntaxhighlight lang="lisp" line start=1>
 
/*****************************************************************/
/* 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]);
</syntaxhighlight>
}}
 
<!-------------------------------------------------------------------------------->
{{MyCodeBlock|title=Solving
|text=
Das lineare Gleichungssystem lösen wir mit der Referenz-Auslenkung
 
::<math>\displaystyle \hat{W} = \frac{q_{ref}\,\ell_0^4}{8\,E\,I_{ref}}</math>
 
und
 
::<math>\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}</math>
 
zu
 
::<math>\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}</math>
 
mit den Bezeichnungen nach Abschnitt [[Randwertprobleme/Methoden zur Lösung von Randwertproblemen/FEM: Trial-Functions für kubische Ansatz-Polynome|FEM: Trial-Functions für kubische Ansatz-Polynome]].
|code=
<syntaxhighlight lang="lisp" line start=1>
/***************************************************************************/
/* 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))))));
</syntaxhighlight>
}}
 
<!-------------------------------------------------------------------------------->
{{MyCodeBlock|title=Post-Processing
|text=
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
 
<ul>
<li>... für die Auslenkung ''w'':<br/>[[Datei:UEBL-31.png|mini|none|Biegelinie w(x)]]</li>
<li>... für den Kippwinkel ''ϕ'':<br/>[[Datei:UEBL-32.png|mini|none|Kippwinkel ''ϕ(x)'']]</li>
<li>... für das Biegemoment ''M'':<br/>[[Datei:UEBL-33.png|mini|none|Biegemomentenverlauf ''M(x)'']]</li>
<li>... für die Querkraft ''Q'':<br/>[[Datei:UEBL-34.png|mini|none|Querkraftverlauf ''Q(x)'']]</li>
</ul>
|code=
<syntaxhighlight lang="lisp" line start=1>
/***************************************************************************/
/* 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]]) )$
</syntaxhighlight>
}}
 
<hr/>
'''Links'''
* [[Gelöste Aufgaben/UEBI|UEBI]] (analytischer Lösung), [[Gelöste Aufgaben/UEBJ|UEBJ]], [[Gelöste Aufgaben/UEBK|UEBK]]
 
'''Literature'''
* ...

Aktuelle Version vom 19. April 2021, 05:49 Uhr


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 h0=b und h1 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 0 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
  • und die Element-"Rechte Seite"
    .

Hier dürfen also der (mittlere) Querschnittsflächeninhalt Ai, das (mittlere) Flächenmoment Ii 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

Dann sind mit

die relativen Element-Längen

.

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

und

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

.

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

und

zu

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

Literature

  • ...