Gelöste Aufgaben/FEAF

Aus numpedia
Zur Navigation springen Zur Suche springen


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,0>0 .

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

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

ϱAu¨EAu=0

hat Lösungen

u(x,t)=i=1IUisin(κix)cos(ω0,it),

die für

κi=π(2i1)2l0;i1

und

ω0=πEϱ02(i(i1)+14)

unsere Randbedingungen erfüllen.

Für die langsamste Eigenmode ist

ω0,1=π20Eϱ , also ist

T1=2πω0,1=40ϱE

und wir wählen  

tBez:=T .

Der partikulare Lösungsanteil ist

up(x)=ϱg02Eξ(1ξ2);ξ=(x0) .

Die statische Auslenkung am unteren Ende ist demnach

up(0)=ϱg022E=:us .

Wir wählen  

Bez:=us .

Modell-Parameter

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

I=2

und damit je Element

i=0I.

Wir wählen lineare Ansatzfunktionen je Element, also

ϕ1=1ξ;ϕ2=ξ.

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

Q_t=(U0(t)U1(t)UI(t)).

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

δW=i=1I(δWiaδΠi)+δWRa,

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

δWRa=0

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

δΠi=(δUi1,δUi)K__i(Ui1(t)Ui(t))

mit der Element-Steifigkeitsmatrix

K__i=EAi(1111).

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

δWia,A=(δUi1,δUi)M__i(U¨i1(t)U¨i(t))

mit der Element-Massenmatrix

M__i=ϱAi6(2112).

Analog folgt für die virtuelle Arbeit der Gewichtskraft

δWia,g=(δUi1,δUi)G_i

mit

G_i=ϱAgi2(11).

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

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

K__=EAi(+1100011+1100011+110001+1+1100011)

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:

M__Q¨_+K__Q_=G_.

/*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 Qt setzt sich nun - bei diesem linearen System - additiv aus partikularer Qp und Qh homogener Lösung zusammen:

Q_t(t)=Q_p(t)+Q_h(t).

/**************** 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 Qp statisch, wir suchen nach der Lösung des Gleichungssystems

K__Q_p=G_

und die ist

Q_p=(3gli2ρ2E2gli2ρE)

Mit us also

Q_p=us(341).

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

M__Q¨_+K__Q_=0_

setzen wir an

Q_h=Q^_heλt

und erhalten das Eigenwertproblem

(λ2M__+K__)=:D__Q^_h=0_

mit

den Eigenwertenλ und
den EigenvektorenQ^_h.

Die Berechnung der Eigenwerte λ wird einfacher mit

λ2=Λ

und wir transformieren parallel die Bewegungsgleichungen auf die dimensionslose Zeit

τ:=ttBez mit tBez=T.,

also

Q¨_=d2Q_dτ2d2τdt21T2.

Wir finden zwei Eigenwerte aus der Bedingung det(D__)=0:

Λ1=4.2101Λ2=5.1102

und

λ1=j6.4λ2=j23 mit j:=1.

Für die Matrix

D__(λi)=λ21T2M__+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 λ1:

Q^_h,1=(0.580.82) und

  • für λ2:

Q^_h,2=(0.580.82).

Die homogene Lösung der Bewegungsgleichung lautet damit

Q_h(t)=c1(0.580.82)ej6.4τ+c2(0.580.82)ej23τ

mit den Integrationskonstanten c1 und c2. Durch die komplexen Eigenwerte sind die beiden Integrationskonstanten nun auch komplexwertig, also

c1=cR,1+jcI,1c2=cR,2+jcI,2.

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

U1(0)=0U2(0)=0U˙1(0)=0U˙2(0)=0.

Wir finden

cR,1=2.5gi2ϱE,cI,1=0,cR,2=0.074gi2ϱE,cI,2=0.

und damit

Q_h(t)=(0.751.0)+(0.731.0)cos(6.4τ)+(0.0210.03)cos(23τ).

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

I=2: Animation der Verschiebung des Endpuntkes. I=5: Animation der Verschiebung des Endpunktes.

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

  • ...