Gelöste Aufgaben/FEAF

Aufgabenstellung

Wir erkunden hier das Thema "Schwingungen von Kontinua" und bauen dabei auf die Ergebnisse zur Berechnung

• der statischen Auslenkung mit FEM aus FEAD und
• der Schwingungen des Feder-Masse-Systems aus FEAE

auf.

Lageplan

Gesucht ist die Längsschwingung des Stabes beim Loslassen aus seiner unverformten Referenzlage. Dabei arbeiten wir mit der Methode der Finiten Elemente zum Aufstellen der Bewegungsgleichungen. Die Integeationskonstanten der Lösung passen wir an die Anfangsbedingungen

• keine Anfangs-Auslenkung
• keine Anfangs-Geschwindigkeit

an.

Lösung mit Maxima

Header

Für die mathematische Behandlung - insbesondere der Auflösung quadratischer Gleichungen - setzen wir in Maxima voraus, dass

$E>0,A>0,{\ell }_0>0$ .

Später brauchen wir für die dimensionslose Formulierung noch eine Bezugszeit t_Bez und eine Bezugslänge l_Bez.

Dazu nehmen wir eine "Anleihe" bei der analytischen Lösung des Schwingungsproblems (vgl. Aufgabe SKER):

Analytische Lösung: homohener Lösungsanteil

Analytische Lösung: partikularer Lösungsanteil

Die partielle Bewegungsgleichung des Stabes $\varrho A\cdot \ddot{u}-EA\cdot u^{″}=0$ hat Lösungen $u(x,t)=\sum _{i=1}^I{U}_i\cdot \sin ({\kappa }_i\cdot x)\cdot \cos ({\omega }_{0,i}\cdot t)$, die für ${\kappa }_i=\frac{\pi (2i-1)}{2l_0};i\ge 1$ und ${\omega }_0=\pi \sqrt{\frac{E}{\varrho {\ell }_0^2}(i\cdot (i-1)+\frac{1}{4})}$ unsere Randbedingungen erfüllen. Für die langsamste Eigenmode ist ${\omega }_{0,1}=\frac{\pi }{2{\ell }_0}\cdot \sqrt{\frac{E}{\varrho }}$ , also ist $\begin{array}{ll}{T}_1& =\frac{2\pi }{{\omega }_{0,1}}\\ & =4{\ell }_0\sqrt{\frac{\varrho }{E}}\end{array}$ und wir wählen   $t_{Bez}:=T$ .

Der partikulare Lösungsanteil ist $u_p(x)=\varrho g\frac{{\ell }_0^2}{E}\xi (1-\frac{\xi }{2});\xi =\left(\frac{x}{{\ell }_0}\right)$ . Die statische Auslenkung am unteren Ende ist demnach $\begin{array}{ll}{u}_p({\ell }_0)& =\varrho g\frac{{\ell }_0^2}{2E}\\ & =:u_s\end{array}$ . Wir wählen   ${\ell }_{Bez}:=u_s$ .

Modell-Parameter

Für die Diskretisierung wählen wir als Anzahl der Finiten Elemente

$I=2$

und damit je Element

${\ell }_i=\frac{{\ell }_0}{I}$.

Wir wählen lineare Ansatzfunktionen je Element, also

${\varphi }_1=1-\xi ;{\varphi }_2=\xi .$

Die abhängigen Koordinaten des FEM-Modells sind dann zunächst - bis zum Einarbeiten der Randbedingungen -

${\underset{\_}{Q}}_t=\left(\begin{array}{c}{U}_0(t)\\ U_1(t)\\ ⋮\\ U_I(t)\\ \end{array}\right)$.

📑 toggle code

/*******************************************************/
/* MAXIMA script                                       */
/* version: wxMaxima 16.04.2                           */
/* author: Andreas Baumgart                            */
/* last updated: 2018-01-07                            */
/* ref: FEM, PVMPE using two coordinates               */
/* description: solve by principle virt. Work          */
/*******************************************************/
 
/* settings */
workdir: "C:\\Users\\X...X\\plots\\";
 
/* assumptions */
assume(E>0,A>0,l[0]>0);
 
/* start: include these results                               */
/* from the analytical solution - see "solved problems": SKER */
sol[1] : [kappa=(2*%pi*i-%pi)/(2*l[0]),omega[0]^2=(4*%pi^2*E*i^2-4*%pi^2*E*i+%pi^2*E)/(4*l[0]^2*rho)]$
u[p]   : [u(x)=-(g*rho*x^2-2*l[0]*g*rho*x)/(2*E),
          U[i]=-32/(8*%pi^3*i^3-12*%pi^3*i^2+6*%pi^3*i-%pi^3)]$
abbrev : [omega[0,0]  =(2*%pi)/T,
          omega[0,0]^2=(%pi^2*E)/(4*l[0]^2*rho),
	  omega[0,0]=(%pi*sqrt(E))/(2*l[0]*sqrt(rho)),
	  T=(4*l[0]*sqrt(rho))/sqrt(E)]$
params : [t=T*tau,u[s]=(l[0]^2*g*rho)/(2*E)]$
/* end: include these results                               */

/* choose number of finite elements :*/
I : 2;
params: append(params, [l[i] = l[0]/I]);
 
/* coordinates */
Q[t]: makelist( U[i](t),i,1,I);
 
/* tial function */
phi : [1-xi, xi];

Equilibrium Conditions

Die Gleichgewichtsbedingung mit dem Prinzip der virtuellen Verrückungen setzen sich additiv aus

$\delta W=\sum _{i=1}^I(\delta W_i^a-\delta {\mathrm{\Pi }}_i)+\delta W_R^a$,

zusammen - also den virtuellen Arbeiten je Element zuzüglich von virtuellen Arbeiten am Rand. In unserem Beispiel ist allerdings

$\delta W_R^a=0$

Je Element sind die virtuelle Formänderungsenergie (vgl. Finite Elemente Methode)

$\delta {\mathrm{\Pi }}_i=(\delta U_{i-1},\delta U_i)\cdot {\underset{\_}{\underset{\_}{K}}}_i\cdot \left(\begin{array}{l}{U}_{i-1}(t)\\ U_i(t)\end{array}\right)$

mit der Element-Steifigkeitsmatrix

${\underset{\_}{\underset{\_}{K}}}_i=\frac{EA}{{\ell }_i}\cdot \left(\begin{array}{rr}1& -1\\ -1& 1\end{array}\right)$.

Und je Element sind die virtuelle Arbeit der D'Alembert'schen Trägheitskraft

$\delta W_i^{a,A}=-(\delta U_{i-1},\delta U_i)\cdot {\underset{\_}{\underset{\_}{M}}}_i\cdot \left(\begin{array}{l}{\ddot{U}}_{i-1}(t)\\ {\ddot{U}}_i(t)\end{array}\right)$

mit der Element-Massenmatrix

${\underset{\_}{\underset{\_}{M}}}_i=\frac{\varrho A{\ell }_i}{6}\cdot \left(\begin{array}{rr}2& 1\\ 1& 2\end{array}\right)$.

Analog folgt für die virtuelle Arbeit der Gewichtskraft

$\delta W_i^{a,g}=(\delta U_{i-1},\delta U_i)\cdot {\underset{\_}{G}}_i$

mit

${\underset{\_}{G}}_i=\frac{\varrho Ag{\ell }_i}{2}\cdot \left(\begin{array}{r}1\\ 1\end{array}\right)$.

Die System-Matrizen komponieren wir nun durch Hinzuaddieren der Anteile je Element.

Beispiel: die Gesamt-Steifigkeitsmatrix für I=4:

$\underset{\_}{\underset{\_}{K}}=\frac{EA}{{\ell }_i}\cdot \left(\begin{array}{ccccc}+1& -1& 0& 0& 0\\ -1& 1+1& -1& 0& 0\\ 0& -1& 1+1& -1& 0\\ 0& 0& -1& +1+1& -1\\ 0& 0& 0& -1& 1\end{array}\right)$

Die Beiträge der vier Elemente sind hier in rot, grün, blau und schwarz eingefärbt.

In Matrix-Schreibweise lautet die Bewegungsgleichung des Gesamt-Systems jetzt:

$\underset{\_}{\underset{\_}{M}}\cdot \underset{\_}{\ddot{Q}}+\underset{\_}{\underset{\_}{K}}\cdot \underset{\_}{Q}=\underset{\_}{G}$.

📑 toggle code

/*elemente mass matrix */
M[i] : rho*A*l[i]*makelist(
                  makelist(integrate(
                                 phi[j]*phi[k], xi,0,1) ,j,1,2),
                                                         k,1,2);
 
/*elemente stiffness matrix */
K[i] : E*A/l[i]*makelist(
                makelist(integrate(
                          diff(phi[j],xi)*diff(phi[k],xi),xi,0,1) ,j,1,2),
                                                                   k,1,2);
 
/*elemente load col-matrix */
G[i] : rho*A*l[i]*g*makelist(integrate(phi[j], xi,0,1) ,j,1,2);/* compose total system matrices */
 
M[0] : zeromatrix(I+1,I+1)$
K[0] : zeromatrix(I+1,I+1)$
G[0] : zeromatrix(I+1,  1)$
 
for e:1 thru I do
   for row:1 thru 2 do
      (for col:1 thru 2 do
         (K[0][(e-1)+row][(e-1)+col] : K[0][(e-1)+row][(e-1)+col] + K[i][row][col],
          M[0][(e-1)+row][(e-1)+col] : M[0][(e-1)+row][(e-1)+col] + M[i][row][col]),
      G[0][(e-1)+row][1] : G[0][(e-1)+row][1] + G[i][row])$
 
/* insert boundary conditions */
M[0] : submatrix(1,M[0],1);
K[0] : submatrix(1,K[0],1);
G[0] : submatrix(1,G[0]  );
 
/* control print-out */
print(M[0],"∙",transpose(diff(Q[t],t,2)),"+",K[0],"∙",transpose(Q[t]),"=",G[0])$

Solving

Die Lösung des Anfangswertproblems setzt sich aus zwei Teilen zusammen:

• der partikularen Lösung, die die Rechte Seite "G" erfüllt und
• der homogenen Lösung, die die Rechte Seite "0" erfüllt.

Die Gesamtlösung Q_t setzt sich nun - bei diesem linearen System - additiv aus partikularer Q_p und Q_h homogener Lösung zusammen:

${\underset{\_}{Q}}_t(t)={\underset{\_}{Q}}_p(t)+{\underset{\_}{Q}}_h(t)$.

📑 toggle code

/**************** total solution */
Q[t] : Q[p] + Q[h];

Right-Hand-Side Approach

Die rechte Seite der Bewegungsgleichung G ist nicht zeitabhängig - sie ist statisch. Also ist auch die Lösung Q_p statisch, wir suchen nach der Lösung des Gleichungssystems

$\underset{\_}{\underset{\_}{K}}\cdot {\underset{\_}{Q}}_p=\underset{\_}{G}$

und die ist

${\underset{\_}{Q}}_p=\left(\begin{array}{c}\frac{3g{l}_i^2\rho }{2E}\\ \frac{2g{l}_i^2\rho }{E}\end{array}\right)$

Mit u_s also

${\underset{\_}{Q}}_p=u_s\cdot \left(\begin{array}{c}\frac{3}{4}\\ 1\end{array}\right)$.

📑 toggle code

/***************************************************************************/
/* p a r t i c u l a r   s o l u t i o n                                   */
/***************************************************************************/
Q[p] : linsolve_by_lu(K[0],G[0])[1];

Homogeneous Solution

Zur Lösung der homogenen Bewegungsgleichung

$\underset{\_}{\underset{\_}{M}}\cdot \underset{\_}{\ddot{Q}}+\underset{\_}{\underset{\_}{K}}\cdot \underset{\_}{Q}=\underset{\_}{0}$

setzen wir an

${\underset{\_}{Q}}_h={\underset{\_}{\hat{Q}}}_h\cdot e^{\lambda \cdot t}$

und erhalten das Eigenwertproblem

$\underset{=:\underset{\_}{\underset{\_}{D}}}{\underbrace{({\lambda }^2\underset{\_}{\underset{\_}{M}}+\underset{\_}{\underset{\_}{K}})}}\cdot {\underset{\_}{\hat{Q}}}_h=\underset{\_}{0}$

mit

den Eigenwerten	$\lambda$ und
den Eigenvektoren	${\underset{\_}{\hat{Q}}}_h$.

Die Berechnung der Eigenwerte λ wird einfacher mit

${\lambda }^2=-\mathrm{\Lambda }$

und wir transformieren parallel die Bewegungsgleichungen auf die dimensionslose Zeit

$\tau :=\frac{t}{t_{Bez}}$ mit $t_{Bez}=T.$,

also

$\underset{\_}{\ddot{Q}}=\frac{d^2\underset{\_}{Q}}{d{\tau }^2}\cdot \underset{\frac{1}{T^2}}{\underbrace{\frac{d^2\tau }{d{t}^2}}}$.

Wir finden zwei Eigenwerte aus der Bedingung $\det (\underset{\_}{\underset{\_}{D}})=0$:

$\begin{array}{lr}{\mathrm{\Lambda }}_1=& 4.2\cdot 10^1\\ {\mathrm{\Lambda }}_2=& 5.1\cdot 10^2\\ \end{array}$

und

$\begin{array}{lr}{\lambda }_1=& j\cdot 6.4\\ {\lambda }_2=& j\cdot 23\\ \end{array}\text{ mit }j:=\sqrt{-1}$.

Für die Matrix

$\underset{\_}{\underset{\_}{D}}({\lambda }_i)={\lambda }^2\cdot \frac{1}{T^2}\cdot \underset{\_}{\underset{\_}{M}}+\underset{\_}{\underset{\_}{K}}$

stellen wir fest:

Die Matrix hat

• einen Rang (rank) von 1 - es gibt eine linear abhängige Zeile im Gleichungssystem und - entsprechend -
• einen Rangabfall (nullity)  von 1.

Die Eigenvektoren spannen nun den Nullraum (nullspace) der Matrix auf. Normmert auf die Länge 1 sind das

• für λ₁:

${\underset{\_}{\hat{Q}}}_{h,1}=\left(\begin{array}{c}0.58\\ 0.82\end{array}\right)$ und

• für λ₂:

${\underset{\_}{\hat{Q}}}_{h,2}=\left(\begin{array}{r}0.58\\ -0.82\end{array}\right)$.

Die homogene Lösung der Bewegungsgleichung lautet damit

${\underset{\_}{Q}}_h(t)=c_1\cdot \left(\begin{array}{r}0.58\\ 0.82\end{array}\right)\cdot e^{j\cdot 6.4\cdot \tau }+c_2\cdot \left(\begin{array}{r}0.58\\ -0.82\end{array}\right)\cdot e^{j\cdot 23\cdot \tau }$

mit den Integrationskonstanten c₁ und c₂. Durch die komplexen Eigenwerte sind die beiden Integrationskonstanten nun auch komplexwertig, also

$\begin{array}{l}{c}_1=c_{R,1}+j\cdot c_{I,1}\\ c_2=c_{R,2}+j\cdot c_{I,2}\end{array}$.

📑 toggle code

/***************************************************************************/
/* h o m o g e n e o u s   s o l u t i o n                                 */
/***************************************************************************/

/* direct approach will usually not work for systems with now > 5 */
/* load(eigen); eigensystem: uniteigenvectors(subst(params, invert(M[0]).K[0])); */
/* -> not employed */
 
/**************** homogeneous solution */
/* solve eigenvalue problem */
D : -Lambda*(1/T)^2*M[0] + K[0];
D : subst(params, D);
charPoly : determinant(D);
 
/* eigenvalues */
eigval : solve(subst(abbrev,charPoly)=0,Lambda);
eigval : subst(abbrev, eigval);
 
/* save for later ....*/
eigvalList : makelist(subst(eigval[j],[%lambda[j]=%i*sqrt(Lambda)]),j,1,I);
 
/* eigenvectors */
prop: [rnk = rank(ratsimp(expand(subst(eigval[1],subst(abbrev,D))))),
       nly = nullity(ratsimp(expand(subst(eigval[1],subst(abbrev,D)))))];
 
 
/* approach with "nullspace" works only for small problems .... */
if I<4 then
   (eigvec : makelist(apply (addcol, args (nullspace(
                           ratsimp(expand(subst(abbrev,subst(eigval[i],D))))))),i,1,I))
else
   (eigvec : [],
    for i:1 thru I do
       (Ared : submatrix(I,ratsimp(expand(subst(abbrev,subst(eigval[i],D/(E*A/l[0])))))),
        bred : -float(col(Ared,I)),
        Ared : float(submatrix(Ared,I)),
        ev : addrow(linsolve_by_lu(Ared,bred)[1],[1]),
        eigvec : append(eigvec,[ev])));
 
/* norm eigenvectors (to be of length "1" */
eigvec : float(ratsimp(makelist(eigvec[i]/sqrt(eigvec[i].eigvec[i]),i,1,I)));

Adapt to Initial Conditions

Diese vier reell-wertigen Integrationskonstanten bekommen wir aus vier Anfangsbedingungen für die beiden Massen, hier für I=2

$\begin{array}{l}{U}_1(0)=0\\ U_2(0)=0\\ {\dot{U}}_1(0)=0\\ {\dot{U}}_2(0)=0\end{array}$.

Wir finden

$c_{R,1}=-2.5\frac{g{\ell }_i^2\varrho }{E},c_{I,1}=0,c_{R,2}=-0.074\frac{g{\ell }_i^2\varrho }{E},c_{I,2}=0$.

und damit

${\underset{\_}{Q}}_h(t)=\left(\begin{array}{r}-0.75\\ -1.0\end{array}\right)+\left(\begin{array}{r}-0.73\\ -1.0\end{array}\right)\cdot \cos (6.4\tau )+\left(\begin{array}{r}-0.021\\ 0.03\end{array}\right)\cdot \cos (23\tau )$.

📑 toggle code

/***************************************************************************/
/* t o t a l   ( c o m p o s e d )   s o l u t i o n                       */
/***************************************************************************/
/* adapt to initial values w[1](0)=0, w[2](0)=0 */
Q[h] : sum(subst(eigval[j],subst([%lambda[j]=%i*sqrt(Lambda),i=j],
                           (c[Re,j]+%i*c[Im,j])*eigvec[j]*%e^(%lambda[j]*tau))),j,1,I);
 
/* update Q[t]*/
Q[t] : ''(Q[t]);
 
/* adapt to initial values */
intcons : flatten(makelist([c[Re,i],c[Im,i]], i, 1, I));
intcond : float(append(makelist(subst([tau=0],     Q[t]     )[i][1]=0,i,1,I),
                       makelist(subst([tau=0],diff(Q[t],tau))[i][1]=0,i,1,I)));
 
/* integration constants: */
sol[2] : ratsimp(solve(realpart(intcond),intcons)[1]);
 
Q[t] : float(realpart(expand(subst(params,subst(sol[2],Q[t]/u[s])))))$

Post-Processing

Die FEM-Lösung im Zeitbereich, angepasst an die Anfangsbedingungen, sieht nun so aus:

I=2: Auslenkung u1(t), u2(t).

I=5: Aulenkungen ui(t).

Interessant ist die Auftragung von analytischer und FEM-Lösung als Animation über der Zeit: Man erkennt, wie die FE-Lösung sowohl Form als auch Zeitverlauf der analytischer Lösung erfasst, man erkennt jedoch auch, wie die FE-Lösung - Aufgrund ihrer höheren Eigenfrequenz - der analytischen Lösung vorauseilt.

📑 toggle code

/* plot results */
 
legends: append([legend], makelist(simplode(["U[",i,"]"]),i,1,I));
plot2d(makelist(Q[t][i][1],i,1,I),[tau,0,4], legends,
    [xlabel, "τ/1 →"], [ylabel, "U/u[s] →"]);
 
/* analytiv and FEM Solution */
analytic : append(sol[1],
           [omega[0]*t = sqrt(ratsimp(subst(abbrev,(subst(sol[1],omega[0]^2*T^2)))))*tau,
            x=l[0]*xi]);
sol[3] : append([ratsimp(subst([x=xi*l[0]],subst(params,subst(u[p],u(x)/u[s]))))],
makelist(subst([i=j],subst(u[p],subst(analytic,U[i]*sin(kappa*x)*cos(omega[0]*t)))),j,1,10));
 
timeFcts : makelist(subst(eigvalList[i],%e^(%lambda[i]*tau)),i,1,I);
Q[l] : subst(params,append([Q[p]/u[s]],
            realpart(makelist(coeff(subst(sol[2],Q[h]),timeFcts[i])/u[s]*timeFcts[i],i,1,I))));
 
sol[4] : append( [[0,0]], makelist([i, Q[t][i][1]],i,1,I));
 
toPlot: [subst([xi=xi/I], sum(sol[3][i],i,1,length(sol[3]))),
         [discrete, sol[4]]];
 
NSteps : 19$
fpprintprec: 4$
for j:0 thru 10*NSteps-1 do
   (tstep : simplode([if j<10 then "00" elseif j<100 then "0" else "",j]),
    plot2d(subst([tau=j/NSteps],toPlot),[xi,0,I], [y,-0.1,2.1],
                     [title, simplode(["t/T=",float(j/NSteps+0.01)])], [legend, "analytic", "FEM"],
                     [style, [lines,1,1], [lines,1,2]], [xlabel, "ξ →"], [ylabel, "u →"],
                     [png_file,simplode([workdir,"FEM-",I,"-step-",tstep,".png"])]))$

Links

• ...

Literature

• ...