Gelöste Aufgaben/FEAD: Unterschied zwischen den Versionen

Aus numpedia
Zur Navigation springen Zur Suche springen
 
(6 dazwischenliegende Versionen desselben Benutzers werden nicht angezeigt)
Zeile 9: Zeile 9:


<onlyinclude>
<onlyinclude>
[[Datei:Screenshot 20210111-063733~2.png|100px|left|mini|Caption]]
[[Datei:FEAB-01.png|mini|186x186px|Lageplan|alternativtext=|links]]
[[Datei:FEAB-01.png|mini|175x175px|Lageplan]]
Gesucht ist die Näherungslösung für die statische Auslenkung der Stab-Querschnitte mit dem Prinzip der Virtuellen Verrückungen und dem Ansatz nach der [[Randwertprobleme/Methoden zur Lösung von Randwertproblemen/Finite Elemente Methode|Methode der Finiten Elemente]]. Dazu wir unterteilen die Struktur in Elemente und setzen [[Randwertprobleme/Methoden zur Lösung von Randwertproblemen/Finite Elemente Methode/FEM: Trial-Functions für lineare Ansatz-Polynome|Lineare Trialfunctions]] für die Verformung in den Elemeten an.
Gesucht ist die Näherungslösung für die statische Auslenkung der Stab-Querschnitte mit dem Prinzip der Virtuellen Verrückungen und dem Ansatz nach der [[Randwertprobleme/Methoden zur Lösung von Randwertproblemen/Finite Elemente Methode|Methode der Finiten Elemente]]. Dazu wir unterteilen die Struktur in Elemente und setzen [[Randwertprobleme/Methoden zur Lösung von Randwertproblemen/Finite Elemente Methode/FEM: Trial-Functions für lineare Ansatz-Polynome|Lineare Trialfunctions]] für die Verformung in den Elemeten an.
</onlyinclude>


==Lösung mit Matlab==
==Lösung mit Matlab==
Lorem Ipsum ....
<!-------------------------------------------------------------------------------->{{MyCodeBlock|title=Header
 
|text=Die Lösung basiert auf Maxima 16.04.2 - wir interessieren uns also vor allem für die Struktur der Lösung.
==tmp==
<!-------------------------------------------------------------------------------->Die Lösung basiert auf Maxima 16.04.2 - wir interessieren uns also vor allem für die Struktur der Lösung.{{MyCodeBlock|title=Header
|text=Text
|code=
|code=
<syntaxhighlight lang="lisp" line start=1>
<syntaxhighlight lang="lisp" line start=1>
1+1
/*******************************************************/
/* MAXIMA script                                      */
/* version: wxMaxima 16.04.2                          */
/* author: Andreas Baumgart                            */
/* last updated: 2017-10-10                            */
/* ref: FENV step 4 im Prozess: FE-Ansätze            */
/* description: mit dem PvV werden die Bewegungsgl.    */
/*          für einen Stab unter Gewichtskraft erstellt*/
/*******************************************************/
</syntaxhighlight>
</syntaxhighlight>
}}
}}
==tmp==
<!-------------------------------------------------------------------------------->{{MyCodeBlock|title=Declarations
<!-------------------------------------------------------------------------------->Für Maxima brauchen wir einige Deklarationen.
|text=Für Maxima brauchen wir einige Deklarationen.


Die Anzahl der Finiten-Elementen ''I'' wählen wir zu
Die Anzahl der Finiten-Elementen ''I'' wählen wir zu


<math>I=3</math>
::<math>I=3</math>


Wir wählen Elemente gleicher Länge, also ist die Elementlänge ''l<sub>i</sub>'' im Element ''i''
Wir wählen Elemente gleicher Länge, also ist die Elementlänge ''l<sub>i</sub>'' im Element ''i''


<math>\begin{array}{ll}\ell_i &= \ell_e\\ & = \ell_0 /I \end{array}</math>.{{MyCodeBlock|title=Declarations
::<math>\begin{array}{ll}\ell_i &= \ell_e\\ & = \ell_0 /I \end{array}</math>.
|text=Text
|code=
|code=
<syntaxhighlight lang="lisp" line start=1>
<syntaxhighlight lang="lisp" line start=1>
1+1
/*******************************************************/
/* declare variational variables - see 6.3 Identifiers */
declare("δW", alphabetic);
declare("δA", alphabetic);
declare("δΠ", alphabetic);
declare("δQ", alphabetic);
declare("δu", alphabetic);
 
/*******************************************************/
/* parameter: Number of Elements */
I : 3; /* change as required */
</syntaxhighlight>
</syntaxhighlight>
}}
}}
==tmp==
<!-------------------------------------------------------------------------------->Wir können die Trial-Functions als


<math>\displaystyle u(x) = \sum_{i=0}^I U_i\cdot\phi_i(x)</math>
<!-------------------------------------------------------------------------------->{{MyCodeBlock|title=Formfunctions
|text=
Wir können die Trial-Functions als
 
::<math>\displaystyle u(x) = \sum_{i=0}^I U_i\cdot\phi_i(x)</math>


mit
mit


<math>\displaystyle \phi_i(x) = \left\{ \begin{array}{cl} \xi_{i-1} & \text{ wobei } x = \ell_e\cdot\left((i-1) + \xi_{i-1}\right),\; 0< \xi_{i-1}<1 \\ \left(1- \xi_{i}\right)  & \text{ wobei } x = \ell_e\cdot\left(i+ \xi_{i}\right),\; 0< \xi_{i}<1 \end{array}\right.</math>.
::<math>\displaystyle \phi_i(x) = \left\{ \begin{array}{cl} \xi_{i-1} & \text{ wobei } x = \ell_e\cdot\left((i-1) + \xi_{i-1}\right),\; 0< \xi_{i-1}<1 \\ \left(1- \xi_{i}\right)  & \text{ wobei } x = \ell_e\cdot\left(i+ \xi_{i}\right),\; 0< \xi_{i}<1 \end{array}\right.</math>.
 
Für ''i=2'' sieht die Trial-Function dann so aus:[[Datei:FEAD-trialfct.png|mini|Trial-Function für ''i=2''.]]<pre>


Für ''i=2'' sieht die Trial-Function dann so aus:
[[Datei:FEAD-trialfct.png|mini|Trial-Function für ''i=2''.]]<br clear="all"/>
<!-- plot commands
plot2d([[parametric, 1+t, t, [t,0,1]],[parametric, 2+t, 1-t, [t,0,1]]],
plot2d([[parametric, 1+t, t, [t,0,1]],[parametric, 2+t, 1-t, [t,0,1]]],


Zeile 56: Zeile 73:


             [xlabel, "x/l ->"], [ylabel, "ϕ2/1"])
             [xlabel, "x/l ->"], [ylabel, "ϕ2/1"])
-->
Für die Methode der Finiten Elemente ist allerdings die Sichtweise je Elemement anschaulicher, also


</pre>Für die Methode der Finiten Elemente ist allerdings die Sichtweise je Elemement anschaulicher, also
::<math>\displaystyle u_i(x_i) = \left\{ \begin{array}{c} \displaystyle U_{i-1} \cdot \left(1- \left(\frac{x_i}{\ell_i} \right)\right) + U_{i} \cdot \left(\frac{x_i}{\ell_i}\right) \\ 0 \text{ sonst }\end{array}\right.</math>,
 
<math>\displaystyle u_i(x_i) = \left\{ \begin{array}{c} \displaystyle U_{i-1} \cdot \left(1- \left(\frac{x_i}{\ell_i} \right)\right) + U_{i} \cdot \left(\frac{x_i}{\ell_i}\right) \\ 0 \text{ sonst }\end{array}\right.</math>,


Und für Element ''i=2'' wird der Verlauf der Querschnitts-Verschiebung also durch
Und für Element ''i=2'' wird der Verlauf der Querschnitts-Verschiebung also durch


<math>\displaystyle u_2(x_2) = \left\{ \begin{array}{c} \displaystyle U_{1} \cdot \left(1- \left(\frac{x_2}{\ell_2} \right)\right) + U_{2} \cdot \left(\frac{x_2}{\ell_2}\right) \\ 0 \text{ sonst }\end{array}\right.</math>.
::<math>\displaystyle u_2(x_2) = \left\{ \begin{array}{c} \displaystyle U_{1} \cdot \left(1- \left(\frac{x_2}{\ell_2} \right)\right) + U_{2} \cdot \left(\frac{x_2}{\ell_2}\right) \\ 0 \text{ sonst }\end{array}\right.</math>.


entsprechend
entsprechend


<math>\displaystyle \delta u_2(x_2) = \left\{ \begin{array}{c} \displaystyle \delta U_{1} \cdot \left(1- \left(\frac{x_2}{\ell_2} \right)\right) + \delta U_{2} \cdot \left(\frac{x_2}{\ell_2}\right) \\ 0 \text{ sonst }\end{array}\right.</math>
::<math>\displaystyle \delta u_2(x_2) = \left\{ \begin{array}{c} \displaystyle \delta U_{1} \cdot \left(1- \left(\frac{x_2}{\ell_2} \right)\right) + \delta U_{2} \cdot \left(\frac{x_2}{\ell_2}\right) \\ 0 \text{ sonst }\end{array}\right.</math>


beschreiben.
beschreiben.
Zeile 73: Zeile 90:
Bei 3 Finiten Elementen sind die System-Koordinaten dann
Bei 3 Finiten Elementen sind die System-Koordinaten dann


<math>\underline{Q} = \left(\begin{array}{l}U_0\\U_1\\U_2\\U_3\end{array}\right)</math>{{MyCodeBlock|title=Formfunctions
::<math>\underline{Q} = \left(\begin{array}{l}U_0\\U_1\\U_2\\U_3\end{array}\right)</math>
|text=Text
|code=
|code=
<syntaxhighlight lang="lisp" line start=1>
<syntaxhighlight lang="lisp" line start=1>
1+1
/* coordinates and their variations */
Q : makelist( U[i],i,0,I);
δQ : makelist(δU[i],i,0,I);
 
/* trial functions */
Phi : [1-(x/l[i]),(x/l[i])];
 
/* Ansatz for element i */
u[i]: sum( U[i+j-2]*Phi[j],j,1,2);
δu[i]: sum(δU[i+j-2]*Phi[j],j,1,2);
</syntaxhighlight>
</syntaxhighlight>
}}
}}
==tmp==
<!-------------------------------------------------------------------------------->Die gesamte virtuelle Arbeit am System können wir nun als Summe der virtuellen Einzelarbeiten je Element hinschreiben, also


<math>\begin{array}{ll}\delta W & = \delta W^a - \delta\Pi\\ &= \sum_{i=1}^I \left( \delta W^a_i - \delta\Pi_i \right) \end{array}</math>.
<!-------------------------------------------------------------------------------->
{{MyCodeBlock|title=Element-wise contributions to ''δW''
|text=
Die gesamte virtuelle Arbeit am System können wir nun als Summe der virtuellen Einzelarbeiten je Element hinschreiben, also
 
::<math>\begin{array}{ll}\delta W & = \delta W^a - \delta\Pi\\ &= \sum_{i=1}^I \left( \delta W^a_i - \delta\Pi_i \right) \end{array}</math>.


Für jedes Element erhalten wir die virtuelle Formänderungsenergie
Für jedes Element erhalten wir die virtuelle Formänderungsenergie


<math>\begin{array}{ll} \displaystyle \delta \Pi_i &=  \int_{0}^{\ell_i} EA \cdot u_i' \cdot \delta u_i' \; dx_i\\ \displaystyle            &= \left(\delta U_{i-1},\delta U_{i}\right)\cdot                            \displaystyle  \frac{EA}{\ell_i}\cdot                            \left(\begin{array}{rr}1&-1\\-1&1\end{array}\right)\cdot                            \left(\begin{array}{l}U_{i-1}\\U_{i}\end{array}\right) \end{array}                            </math>
::<math>\begin{array}{ll} \displaystyle \delta \Pi_i &=  \int_{0}^{\ell_i} EA \cdot u_i' \cdot \delta u_i' \; dx_i\\ \displaystyle            &= \left(\delta U_{i-1},\delta U_{i}\right)\cdot                            \displaystyle  \frac{EA}{\ell_i}\cdot                            \left(\begin{array}{rr}1&-1\\-1&1\end{array}\right)\cdot                            \left(\begin{array}{l}U_{i-1}\\U_{i}\end{array}\right) \end{array}                            </math>


und die virtuelle Arbeit der äußeren, eingeprägten Gewichtskraft
und die virtuelle Arbeit der äußeren, eingeprägten Gewichtskraft


<math>\begin{array}{ll} \displaystyle \delta W^a_i &=  \int_0^{\ell_i} \varrho \cdot A \cdot g \cdot \delta u_i \; dx_i\\ \displaystyle              &= \frac{1}{2} \cdot \left(\delta U_{i-1},\delta U_{i}\right)\cdot                              \varrho \cdot A \cdot g \cdot \ell_i \cdot            \left(\begin{array}{l}1\\1\end{array}\right) \end{array}                        </math>.
::<math>\begin{array}{ll} \displaystyle \delta W^a_i &=  \int_0^{\ell_i} \varrho \cdot A \cdot g \cdot \delta u_i \; dx_i\\ \displaystyle              &= \frac{1}{2} \cdot \left(\delta U_{i-1},\delta U_{i}\right)\cdot                              \varrho \cdot A \cdot g \cdot \ell_i \cdot            \left(\begin{array}{l}1\\1\end{array}\right) \end{array}                        </math>.
 
{{MyCodeBlock|title=Element-wise contributions to ''δW''
|text=Text
|code=
|code=
<syntaxhighlight lang="lisp" line start=1>
<syntaxhighlight lang="lisp" line start=1>
1+1
/* Equilibrium */
δΠ[i] : E*A*integrate(diff(u[i],x)*diff(δu[i],x),x,0,l[i]);
δA[i] : integrate(rho*g*A*δu[i], x,0,l[i]);
 
K[i] : funmake('matrix, makelist(makelist(coeff(coeff(expand( δΠ[i]), δU[i+j]),U[i+k]), j,-1,0), k,-1,0));
P[i] : funmake('matrix,          makelist(    [coeff(expand( δA[i]), δU[i+j])]      , j,-1,0)        );
</syntaxhighlight>
</syntaxhighlight>
}}
}}
==tmp==
 
<!-------------------------------------------------------------------------------->Für jedes Element müssen wir nun die virtuellen Arbeiten zum Gesamt-Gleichungssystem zusammenaddieren.
<!-------------------------------------------------------------------------------->{{MyCodeBlock|title=Equlilibrium Conditions
|text=
Für jedes Element müssen wir nun die virtuellen Arbeiten zum Gesamt-Gleichungssystem zusammenaddieren.


Die Gleichgewichtsbedingung
Die Gleichgewichtsbedingung


<math>\begin{array}{ll}\delta W &= \delta W^a - \delta \Pi\\ &\stackrel{!}{=}0\end{array}</math>
::<math>\begin{array}{ll}\delta W &= \delta W^a - \delta \Pi\\ &\stackrel{!}{=}0\end{array}</math>


liefern ein Gleichungssystem, in das wir nun die kinematische Randbedingung ''U<sub>0</sub>=0, δU<sub>0</sub>=0'' durch Streichen der ersten Zeile des Gleichungssystems und der ersten Spalte der Gesamt-Steifigkeitsmatrix einarbeiten.
liefern ein Gleichungssystem, in das wir nun die kinematische Randbedingung ''U<sub>0</sub>=0, δU<sub>0</sub>=0'' durch Streichen der ersten Zeile des Gleichungssystems und der ersten Spalte der Gesamt-Steifigkeitsmatrix einarbeiten.
Zeile 111: Zeile 143:
Wir wählen nun noch für jedes Element die gleiche Element-Länge ''ℓ<sub>i</sub> = ℓ<sub>e</sub>'', mit ''ℓ<sub>e</sub> = ℓ<sub>0</sub>/3'' und es bleibt
Wir wählen nun noch für jedes Element die gleiche Element-Länge ''ℓ<sub>i</sub> = ℓ<sub>e</sub>'', mit ''ℓ<sub>e</sub> = ℓ<sub>0</sub>/3'' und es bleibt


<math>\displaystyle \frac{E A}{\ell_i}\begin{pmatrix}2 & -1 & 0\\ -1 & 2 & -1\\ 0 & -1 & 1\end{pmatrix}\cdot\left(\begin{array}{l}U_1\\U_2\\U_3\end{array}\right) = \varrho A \ell_i g \cdot\left(\begin{array}{l}1\\1\\\frac{1}{2}\end{array}\right)</math>.{{MyCodeBlock|title=Equlilibrium Conditions
::<math>\displaystyle \frac{E A}{\ell_i}\begin{pmatrix}2 & -1 & 0\\ -1 & 2 & -1\\ 0 & -1 & 1\end{pmatrix}\cdot\left(\begin{array}{l}U_1\\U_2\\U_3\end{array}\right) = \varrho A \ell_i g \cdot\left(\begin{array}{l}1\\1\\\frac{1}{2}\end{array}\right)</math>.
|text=Text
|code=
|code=
<syntaxhighlight lang="lisp" line start=1>
<syntaxhighlight lang="lisp" line start=1>
1+1
/* Assembling the System Equations */
 
K[0] : zeromatrix(length(Q),length(Q))$
P[0] : zeromatrix(length(Q),    1    )$
 
for j:1 thru I do
  (for row:1 thru 2 do
      (P[0][j-1+row][  1    ]: P[0][j-1+row][  1    ]+P[i][row][ 1 ],
        for col:1 thru 2 do
          K[0][j-1+row][j-1+col]: K[0][j-1+row][j-1+col]+K[i][row][col]));
 
/* Incorporate Boundary Condition U[0]=0 */
K[0] : submatrix(1,K[0],1);
P[0] : submatrix(1,P[0]);
</syntaxhighlight>
</syntaxhighlight>
}}
}}
==tmp==
<!-------------------------------------------------------------------------------->Die Gleichgewichtsbedingungen lösen wir nun und erhalten


<math>\displaystyle {{U}_{0}}=0,{{U}_{1}}=\frac{5g\,{{\ell}_{i}^{2}}\varrho}{2E},{{U}_{2}}=\frac{4g\,{{\ell}_{i}^{2}}\varrho}{E},{{U}_{3}}=\frac{9g\,{{\ell}_{i}^{2}}\varrho}{2E}</math>{{MyCodeBlock|title=Solving
<!-------------------------------------------------------------------------------->{{MyCodeBlock|title=Solving
|text=Text
|text=
Die Gleichgewichtsbedingungen lösen wir nun und erhalten
 
::<math>\displaystyle {{U}_{0}}=0,{{U}_{1}}=\frac{5g\,{{\ell}_{i}^{2}}\varrho}{2E},{{U}_{2}}=\frac{4g\,{{\ell}_{i}^{2}}\varrho}{E},{{U}_{3}}=\frac{9g\,{{\ell}_{i}^{2}}\varrho}{2E}</math>
|code=
|code=
<syntaxhighlight lang="lisp" line start=1>
<syntaxhighlight lang="lisp" line start=1>
1+1
/* solution */
sol[0]: linsolve_by_lu(K[0], P[0]);
sol[1]: append([U[0]=0],makelist(U[i]=sol[0][1][i][1],i,1,I))
</syntaxhighlight>
</syntaxhighlight>
}}
}}
==tmp==
 
<!-------------------------------------------------------------------------------->Wir tragen die Ergebnisse für die numerische Näherungslösung gegen die exakte Lösung auf. Dabei setzen wir l<sub>i</sub>=l<sub>0</sub>/3. Die Verschiebungen sind elementweise für die Elemente e = 1,2,3 aufgetragen:[[Datei:FEAD-plot.png|mini|Verschiebung u(x) - analytische und FE-Lösung]]
<!-------------------------------------------------------------------------------->
{{MyCodeBlock|title=Post-processing
{{MyCodeBlock|title=Post-processing
|text=Text
|text=
Wir tragen die Ergebnisse für die numerische Näherungslösung gegen die exakte Lösung auf. Dabei setzen wir l<sub>i</sub>=l<sub>0</sub>/3. Die Verschiebungen sind elementweise für die Elemente e = 1,2,3 aufgetragen:[[Datei:FEAD-plot.png|mini|Verschiebung u(x) - analytische und FE-Lösung]]
|code=
|code=
<syntaxhighlight lang="lisp" line start=1>
<syntaxhighlight lang="lisp" line start=1>
1+1
/* plot results*/
/* define piecewise, normalized functions u[0]    */
u[0]: subst(sol[1],makelist(subst([i=j],subst([x=xi*l[i]],u[i])),j,1,I))$
/* normalize */
u[0]: expand(subst([l[i]=l[0]/I],u[0])/(rho*g*l[0]^2/E));
 
/* analytic solution */
analytic : subst([xi=xi/I],(xi)*(1-(xi)/2))$
toPlot: append(        [[parametric,    t, analytic, [t, 0, I]]],
              makelist([parametric, i-1+t,  u[0][i], [t, 0, 1]],i,1,I))$
/* use t as independent variable for parametric plot */
toPlot: subst([xi=t], toPlot)$
 
/* plot */
legende : append([legend, "analytic"],makelist(simplode(["e = ",i]),i,1,I))$
plot2d(toPlot, legende,
              [title, sconcat("number of elements I = ",I)],
              [xlabel, "x/l[e] ->"], [ylabel, "u(x)/(rho*g*l/E)"],
              [gnuplot_preamble,"set key bottom"])$
</syntaxhighlight>
</syntaxhighlight>
}}
}}
Zeile 147: Zeile 213:


'''Links'''
'''Links'''
*...
*[[Gelöste Aufgaben/FEAB]]
*[[Gelöste Aufgaben/FEAC]]


'''Literature'''
'''Literature'''
*...
*[[Sources/Literatur#Strang2008|Strang 2008]]

Aktuelle Version vom 9. März 2021, 11:17 Uhr

Aufgabenstellung

Statt Formfunktionen über die ganze Stablänge anzusetzten wie in FEAB gehen wir jetzt nach der Methode der Finiten Elemente vor.


Lageplan

Gesucht ist die Näherungslösung für die statische Auslenkung der Stab-Querschnitte mit dem Prinzip der Virtuellen Verrückungen und dem Ansatz nach der Methode der Finiten Elemente. Dazu wir unterteilen die Struktur in Elemente und setzen Lineare Trialfunctions für die Verformung in den Elemeten an.


Lösung mit Matlab

Header

Die Lösung basiert auf Maxima 16.04.2 - wir interessieren uns also vor allem für die Struktur der Lösung.


/*******************************************************/
/* MAXIMA script                                       */
/* version: wxMaxima 16.04.2                           */
/* author: Andreas Baumgart                            */
/* last updated: 2017-10-10                            */
/* ref: FENV step 4 im Prozess: FE-Ansätze             */
/* description: mit dem PvV werden die Bewegungsgl.    */
/*          für einen Stab unter Gewichtskraft erstellt*/
/*******************************************************/



Declarations

Für Maxima brauchen wir einige Deklarationen.

Die Anzahl der Finiten-Elementen I wählen wir zu

Wir wählen Elemente gleicher Länge, also ist die Elementlänge li im Element i

.

/*******************************************************/
/* declare variational variables - see 6.3 Identifiers */
declare("δW", alphabetic);
declare("δA", alphabetic);
declare("δΠ", alphabetic);
declare("δQ", alphabetic);
declare("δu", alphabetic);

/*******************************************************/
/* parameter: Number of Elements */
I : 3; /* change as required */




Formfunctions

Wir können die Trial-Functions als

mit

.

Für i=2 sieht die Trial-Function dann so aus:

Trial-Function für i=2.


Für die Methode der Finiten Elemente ist allerdings die Sichtweise je Elemement anschaulicher, also

,

Und für Element i=2 wird der Verlauf der Querschnitts-Verschiebung also durch

.

entsprechend

beschreiben.

Bei 3 Finiten Elementen sind die System-Koordinaten dann


/* coordinates and their variations */
 Q : makelist( U[i],i,0,I);
δQ : makelist(δU[i],i,0,I);

/* trial functions */
Phi : [1-(x/l[i]),(x/l[i])];

/* Ansatz for element i */
 u[i]: sum( U[i+j-2]*Phi[j],j,1,2);
δu[i]: sum(δU[i+j-2]*Phi[j],j,1,2);




Element-wise contributions to δW

Die gesamte virtuelle Arbeit am System können wir nun als Summe der virtuellen Einzelarbeiten je Element hinschreiben, also

.

Für jedes Element erhalten wir die virtuelle Formänderungsenergie

und die virtuelle Arbeit der äußeren, eingeprägten Gewichtskraft

.

/* Equilibrium */
δΠ[i] : E*A*integrate(diff(u[i],x)*diff(δu[i],x),x,0,l[i]);
δA[i] : integrate(rho*g*A*δu[i], x,0,l[i]);

K[i] : funmake('matrix, makelist(makelist(coeff(coeff(expand( δΠ[i]), δU[i+j]),U[i+k]), j,-1,0), k,-1,0));
P[i] : funmake('matrix,          makelist(     [coeff(expand( δA[i]), δU[i+j])]       , j,-1,0)         );




Equlilibrium Conditions

Für jedes Element müssen wir nun die virtuellen Arbeiten zum Gesamt-Gleichungssystem zusammenaddieren.

Die Gleichgewichtsbedingung

liefern ein Gleichungssystem, in das wir nun die kinematische Randbedingung U0=0, δU0=0 durch Streichen der ersten Zeile des Gleichungssystems und der ersten Spalte der Gesamt-Steifigkeitsmatrix einarbeiten.

Wir wählen nun noch für jedes Element die gleiche Element-Länge i = ℓe, mit e = ℓ0/3 und es bleibt

.

/* Assembling the System Equations */

K[0] : zeromatrix(length(Q),length(Q))$
P[0] : zeromatrix(length(Q),    1    )$

for j:1 thru I do
   (for row:1 thru 2 do
       (P[0][j-1+row][  1    ]: P[0][j-1+row][  1    ]+P[i][row][ 1 ],
        for col:1 thru 2 do
           K[0][j-1+row][j-1+col]: K[0][j-1+row][j-1+col]+K[i][row][col]));

/* Incorporate Boundary Condition U[0]=0 */
K[0] : submatrix(1,K[0],1);
P[0] : submatrix(1,P[0]);




Solving

Die Gleichgewichtsbedingungen lösen wir nun und erhalten


/* solution */
sol[0]: linsolve_by_lu(K[0], P[0]);
sol[1]: append([U[0]=0],makelist(U[i]=sol[0][1][i][1],i,1,I))




Post-processing

Wir tragen die Ergebnisse für die numerische Näherungslösung gegen die exakte Lösung auf. Dabei setzen wir li=l0/3. Die Verschiebungen sind elementweise für die Elemente e = 1,2,3 aufgetragen:

Verschiebung u(x) - analytische und FE-Lösung

/* plot results*/
/* define piecewise, normalized functions u[0]     */
u[0]: subst(sol[1],makelist(subst([i=j],subst([x=xi*l[i]],u[i])),j,1,I))$
/* normalize */
u[0]: expand(subst([l[i]=l[0]/I],u[0])/(rho*g*l[0]^2/E));

/* analytic solution */
analytic : subst([xi=xi/I],(xi)*(1-(xi)/2))$
toPlot: append(        [[parametric,     t, analytic, [t, 0, I]]], 
               makelist([parametric, i-1+t,  u[0][i], [t, 0, 1]],i,1,I))$
/* use t as independent variable for parametric plot */
toPlot: subst([xi=t], toPlot)$

/* plot */
legende : append([legend, "analytic"],makelist(simplode(["e = ",i]),i,1,I))$
plot2d(toPlot, legende,
               [title, sconcat("number of elements I = ",I)],
               [xlabel, "x/l[e] ->"], [ylabel, "u(x)/(rho*g*l/E)"],
               [gnuplot_preamble,"set key bottom"])$




Spannungen im Stab:
Tragen Sie auch die Spannungen im Stab über die Stablänge an! Berechnen Sie die Spannungen auf Basis der Dehnung 




Links

Literature