Gelöste Aufgaben/FEAG

Aufgabenstellung

Analog zu FEAF untersuchen wir hier die Schwingungen eines Kontinuums beim Loslassen aus der entspannten Rugelage. Hier nicht mit einem Dehnstab, sondern einem Euler-Bernoulli-Balken.

Lageplan

Gesucht ist die Schwingung eines Euler-Bernoulli-Balkens beim Loslassen aus der Ruhelage. Wir gehen nach dem Standardrezept der Finite Elemente Methode vor, arbeiten also mit dem Prinzip der virtuellen Arbeit.

Lösung mit Matlab^®

Interessant ist hier, dass - im Gegensatz zu Stablängsschwingungen - die Eigenfrequenz nicht ein gerades Vielfaches der untersten Eigenfrequenz ist. Falls Sie ein Saiteninstrument spielen, verstehen Sie sofort, warum das wichtig ist.

Maxima können wir hier nicht gut gebrauchen: die Gleichungen werden zu umfangreich. Wir arbeiten also mehr mit numerischen Verfahren, da ist Matlab geeigneter.

Header

Im Programm arbeiten wir mit einer dimensionslosen Formulierung - wir brauchen dafür eine Bezugszeit t_Bez und eine Bezugslänge l_Bez.

Bezugsgrößen wählen

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

Analytische Lösung: homohener Lösungsanteil

Analytische Lösung: partikularer Lösungsanteil

Die homogene Bewegungsgleichung des Euler-Bernoulli-Balkens $\varrho \,A\cdot {\ddot {w}}+E\,I\cdot {w}^{IV}=0$ hat Lösungen vom Typ $\displaystyle w(x,t)=\sum _{i=1}^{\infty }W_{i}\cdot e^{\displaystyle \kappa _{i}\cdot x}\cdot e^{\displaystyle \omega _{0,i}\cdot t}$ Wir bekommen zu jedem ω0 vier κ ${\begin{array}{l}\displaystyle \kappa _{1}=+{\frac {j\,\alpha }{\ell }},\\\displaystyle \kappa _{2}=-{\frac {\alpha }{\ell }},\\\displaystyle \kappa _{3}=-{\frac {j\,\alpha }{\ell }},\\\displaystyle \kappa _{4}=+{\frac {\alpha }{\ell }}\end{array}}\;\;{\text{ mit }}\;\;\displaystyle \alpha ^{4}={\frac {{{\ell }^{4}}\,A\,\rho }{E\,I}}\cdot {{\omega }_{0}^{2}}$ Das α ist eine praktische Abkürzung, hinter der wir ω0 verstecken. Anders als beim Dehnstab (SKER) finden wir hier keine analytische Beziehung, sondern nur die numerischen Beziehungen, für die unsere Randbedingungen erfüllt sind: ${\begin{array}{c}\alpha _{1}=1.875\ldots ,\\\alpha _{2}=4.694\ldots ,\\\alpha _{3}=7.854\ldots ,\\\vdots \end{array}}$ Die langsamste Eigenmode gehört zu α1 mit der Periodendauer $\displaystyle T^{*}=0.5688\ldots \pi \ell ^{2}{\sqrt {\frac {\varrho \,A}{E\,I}}}$. Also wählen wir $t_{Bez}:=T^{*}$.

Die zugeordnete inhomogene Bewegungsgleichung des Euler-Bernoulli-Balkens $E\,I\cdot {w}^{IV}=\varrho \,A\cdot g$ hat die partikulare Lösung $\displaystyle E\,I\,{w}_{p}\left(x\right)=A\,g\,\varrho \,\ell ^{4}\cdot \left({\frac {{\xi }^{4}}{24}}-{\frac {{\xi }^{3}}{6}}+{\frac {{\xi }^{2}}{4}}\right)$ Die statische Auslenkung am unteren Ende ist demnach $\displaystyle {{w}_{s}}={\frac {{{\ell }^{4}}\,A\,g\,\rho }{8\,E\,I}}$. Wir wählen $\ell _{Bez}:=w_{s}$.

System-Parameter des FEM-Modells

Für die Diskretisierung wählen wir als Anzahl der Finiten Elemente (Number Of Elements):

$\displaystyle I=4\;\;\;({\text{ im Matlab-Code: }}NOE)$

und damit je Element

$\displaystyle \ell _{i}={\frac {\ell _{0}}{I}}$.

Wir wählen die kubische Ansatzfunktionen aus dem Abschnitt Finite Elemente Methode je Element, also

${\begin{array}{crc}\phi _{1}=&\left(\xi -1\right)^{2}\cdot \left(1+2\cdot \xi \right)\\\phi _{2}=&\left(\xi -1\right)^{2}\cdot \xi \cdot \ell _{i}\\\phi _{3}=&\left(\xi ^{2}\cdot \left(3-2\cdot \xi \right)\right)\\\phi _{4}=&\left(\left(\xi -1\right)\cdot \xi ^{2}\cdot \ell _{i}\right)\end{array}}$

Für die numerisch Implementierung stört in dieser Darstellung das Element-spezifische ℓ_i - die Elementlänge. Wir behelfen uns, indem wir die Ansatzfunktionen schreiben als

$\displaystyle \underbrace {\left({\begin{array}{c}\phi _{1}\\\phi _{2}\\\phi _{3}\\\phi _{4}\end{array}}\right)} _{\displaystyle =:{\underline {\phi }}}=\underbrace {\left({\begin{array}{cccc}1&0&0&0\\0&\ell _{i}&0&0\\0&0&1&0\\0&0&0&\ell _{i}\\\end{array}}\right)} _{\displaystyle =:{\underline {\underline {d}}}}\cdot \underbrace {\left({\begin{array}{c}\left(\xi -1\right)^{2}\cdot \left(1+2\cdot \xi \right)\\\left(\xi -1\right)^{2}\cdot \xi \\\left(\xi ^{2}\cdot \left(3-2\cdot \xi \right)\right)\\\left(\left(\xi -1\right)\cdot \xi ^{2}\right)\end{array}}\right)} _{\displaystyle =:{\underline {\tilde {\phi }}}}$,

wobei wir die Ansatz-Polynome in ${\underline {\tilde {\phi }}}$ verpacken.

Das schaut unnötig komplex aus - allerdings stecken wir die Diagonal-Matrizen d (für jedes Element eine) in eine Matlab^®-Variable, so dass sie dort nicht weiter auffällt.

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

${\underline {Q}}_{t}=\left({\begin{array}{c}W_{0}(t)\\\Phi _{0}(t)\\W_{1}(t)\\\vdots \\\Phi _{I}(t)\\W_{I}(t)\\\end{array}}\right)$.

Matlab

Classes

In dieser CODE-Sektion kommt ein Beispiele zum Aufbau des Matlab^®-Programms mit Klassen (Objektbasierte Programmierung)

1. FEM-Section

   classdef FEM_Section  % section section parameters properties rho; E; beta; h; b; l; end methods  % Constructor function self = FEM_Section(parHashMap) self.rho = parHashMap('rho'); self.E = parHashMap('E'); self.beta = parHashMap('beta'); self.h = parHashMap('h'); self.b = parHashMap('b'); self.l = parHashMap('l'); end  % end Constructor  % start bending stiffness EI function EI = bendingStiffness(self) EI = self.E*1/12.*self.b*self.h^3; end  % end bending stiffness  % start cross-sectional area function A = cossecArea(self) A = self.b*self.h; end  % end bending stiffness  % start mass function m = mass(self) m = self.rho*self.b*self.h*self.l; end  % end mass end end

2. FEM-Element-Model

   classdef FEM_Element_Modell  % modell parameters properties type;  % element type noe;  % number of elements of this type noc;  % number of coordinates per node phi;  % trial-functions M;  % element-mass matrix K;  % element-stiffness matrix G;  % element-load matrix (g) end methods  % Constructor function self = FEM_Element_Modell(parHashMap) self.type = parHashMap('Type'); self.noe = parHashMap('NOE'); if self.type == 'EBB'  % Euler-Bernoulli-Beam Elements self.phi = [[ 2,-3, 0, 1];... [ 1,-2, 1, 0];... [-2, 3, 0, 0];... [ 1,-1, 0, 0]]; self.noc = 2; elseif self.type == 'ER'  % Extensible Rod type self.phi = [[-1, 1]; ... [ 1, 0]]; self.noc = 1; end [self.M,self.K,self.G] = self.ElementMassMatrix(); end  % end Constructor function [M, K, G] = ElementMassMatrix(self) m = length(self.phi); M = zeros(m,m); K = zeros(m,m); G = zeros(m,1); for row = 1:m for col = 1:m M(row,col) = diff(... polyval(... polyint(... conv(self.phi(row,:),self.phi(col,:))),... [0,1])); if self.type == 'EBB' % -> second derivs K(row,col) = diff(... polyval(... polyint(... conv(polyder(polyder(self.phi(row,:))),polyder(polyder(self.phi(col,:))))),... [0,1])); elseif self.type == 'ER' % -> first derivs K(row,col) = diff(... polyval(... polyint(... conv(polyder(self.phi(row,:)),polyder(self.phi(col,:)))),... [0,1])); end end G(row,1) = diff(... polyval(... polyint(... self.phi(row,:)),... [0,1])); end end end end

3. FEM-Container

   classdef FEM_Container  % holds all system matrics properties  % matrices M;  % mass matrix D;  % damping matrix K;  % stiffness matrix G;  % gravitational Loading d;  % scalting matrix with l[i]-Element  % rod length l;  % length  % model properties NOE; % number of elements NON; % number of nodes NOC  % number of coordinates per node NOQ  % number of coordinates in total g;  % 9.81 m/s^2  % boundary conditions free; disabled;  % functions trials;  % solution samples; % number of samples when plotting with points end methods  % Constructor function self = FEM_Container(bezug,system,boundaries,element)  % declarations NOE = sum(element.noe); % number of all elements in all sections NON = NOE+1;  % number of nodes NOC = element.noc;  % number of coordinates per node NOQ = NON*NOC;  % number of coordintes self.NOE = NOE; self.NON = NON; self.NOC = NOC; self.NOQ = NOQ;  % number of samples when plotting self. samples = 21;  % global params self.g = 9.81/(bezug('l_ref')/bezug('t_ref')^2);  % FE-lengths self.l = [0]; for sec = 1: length(system) for ele = 1: element.noe(sec) self.l = [self.l, self.l(length(self.l))+system(sec).l/element.noe(sec)]; end end  % initialize system matrices [self.M, self.K, self.G, self.d] = self.systemMatrices(system, element);  % collect boundary conditions free = linspace(1,1, element.noc*(sum(element.noe)+1)); touch = [1]; for sec = 1: length(element.noe) touch = [touch,touch(length(touch))+element.noc*element.noe(sec)]; end for i= 1:min(length(touch),length(boundaries)) nodeNo= touch(i); free(nodeNo:nodeNo+element.noc-1) = boundaries(i).free; end disabled = ~free; self.free = (1:length(free)).*free; self.free(self.free==0)=[]; disabled = (1:length(disabled)).*disabled; self.disabled = (1:length(disabled)).*disabled; self.disabled(self.disabled==0)=[];  % trial functions self.trials = element.phi; end  % end Constructor  % start system matrices function [M,K,G,d] = systemMatrices(self, system, element) M = zeros(self.NOQ,self.NOQ); K = zeros(self.NOQ,self.NOQ); G = zeros(self.NOQ, 1 ); d = zeros(4, 4, self.NOE ); g = self.g; e = 0;% counter for all elements of all sections for section = 1: length(element.noe) L = system(section).l/element.noe(section); m = system(section).rho*system(section).cossecArea()*L;  % beinding stiffness b = system(section).bendingStiffness()/L^3; for i = 1:element.noe(section) e = e+1; if element.type=='EBB' d(:,:,e) = diag([1,L,1,L]); % diagnoal matrix with l[i]-factors elseif element.type=='ER' d(:,:,e) = diag([1,1]); end  % j = 2*(e-1)+1; M(j:j+3,j:j+3) = M(j:j+3,j:j+3)+m* d(:,:,e)*element.M*d(:,:,e); K(j:j+3,j:j+3) = K(j:j+3,j:j+3)+b* d(:,:,e)*element.K*d(:,:,e); G(j:j+3,1) = G(j:j+3,1) +m*g*d(:,:,e)*element.G; end end end  % end system matrices  % start boudaryconditions function [M,K,G] = constrained(self) M = self.M(self.free,self.free); K = self.K(self.free,self.free); G = self.G(self.free, 1 ); end  % end system matrices  % start eigensystem function [V,D] = eigensystem(self)  % ODE: A q + B q = 0  % = - = - -  % define system matrices [A,B,G] = self.constrained(); %S = linsolve(A,-B); S = inv(A)*(-B);  % eigensystem [V,D] = eig(S); end  % end eigensystem  % start particular solution function Q = particularSolution(self)  % ODE: A q + B q = G  % = - = - -  % define system matrices [A,B,G] = self.constrained(); Q = linsolve(B,G); end  % end particular solution end end

Eingabe-Parameter aus Excel

0) Bezugsgrößen.

2) Modell

1) System

2) Randbedingungen

Equilibrium Conditions

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

$\displaystyle \delta W=\sum _{i=1}^{I}\left(\delta W_{i}^{a}-\delta \Pi _{i}\right)+\delta W_{R}^{a}$,

zusammen - also den virtuellen Arbeiten je Element zuzüglich von virtuellen Arbeiten am Rand. In unserem Beispiel haben wir allerdings keine eingeprägten, äußeren Lasten am Rand, also ist

$\delta W_{R}^{a}=0$.

Je Element erfassen wir die virtuellen Arbeiten durch Element-Matrizen. So ist z.B. für die virtuelle Formänderungsenergie (vgl. Finite Elemente Methode)

$\delta \Pi _{i}=\left(\delta W_{i-1},\delta \Phi _{i-1},\delta W_{i},\delta \Phi _{i}\right)\cdot {\underline {\underline {K}}}_{i}\cdot \left({\begin{array}{c}W_{i-1}(t)\\\Phi _{i-1}(t)\\W_{i}(t)\\\Phi _{i}(t)\end{array}}\right)$

mit der Element-Steifigkeitsmatrix ${\underline {\underline {K}}}_{i}$.

Aus der FEM-Formulierung für den Euler-Bernoulli-Balken könnten wir die Element-Matrizen der virtuellen Arbeit der D'Alembert'schen Trägheitskraft sowie der Virtuelle Formänderungsenergie herauskopieren. Auch die "rechte Seite" mit der virtuelle Arbeit der Gewichtskraft aus

$\displaystyle \delta W_{i}^{a,g}=\int _{\ell _{i}}\varrho \;A\;g\cdot \delta w(x)\;dx_{i}$,

könnten wir als Element-Lastmatrix (Spaltenmatrix) abspeichern.

$\displaystyle {\underline {G}}_{i}={\frac {\varrho \,A\,g\,\ell _{i}}{2}}\cdot \left({\begin{array}{c}+1\\\displaystyle +{\frac {\ell _{i}}{6}}\\+1\\\displaystyle -{\frac {\ell _{i}}{6}}\end{array}}\right)$.

Dass es auch anders geht,sehen Sie in der Definition der Matlab^®-Klasse FEM_Element_Modell, bei der die Element direkt aus den Trial-Functions - hier den ${\underline {\tilde {\phi }}}$ -

hergeleitet. werden:

                 :
          K(row,col) = diff(...
                            polyval(...
                                 polyint(...
                                     conv(polyder(self.phi(row,:)),polyder(self.phi(col,:)))),...
                                     [0,1]));
                 :

Komponieren der System-Matrizen

Die System-Matrizen des Gesamt–Systems komponieren wir nun durch Hinzuaddieren der Anteile je Element.

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

$\displaystyle {\underline {\underline {K}}}={\frac {E\,I}{\ell _{i}^{3}}}\cdot \left({\begin{array}{ccccc}{\color {red}{+12}}&{\color {red}{+6}}&{\color {red}{-12}}&{\color {red}{+6}}&0&0\\{\color {red}{+6\ell _{i}}}&{\color {red}{+4\ell _{i}^{2}}}&{\color {red}{-6\ell _{i}}}&{\color {red}{+2\ell _{i}^{2}}}&0&0\\{\color {red}{-12}}&{\color {red}{-6}}&{\color {red}{+12}}{\color {green}{+12}}&{\color {red}{-6}}{\color {green}{+6}}&{\color {green}{-12}}&{\color {green}{+6}}\\{\color {red}{+6\ell _{i}}}&{\color {red}{2\ell _{i}^{2}}}&{\color {red}{-6\ell _{i}}}{\color {green}{+6\ell _{i}}}&{\color {red}{+2\ell _{i}^{2}}}{\color {green}{+2\ell _{i}^{2}}}&{\color {green}{-6\ell _{i}}}&{\color {green}{+2\ell _{i}^{2}}}\\0&0&{\color {green}{-12}}&{\color {green}{-6}}&{\color {green}{+12}}&{\color {green}{-6}}\\0&0&{\color {green}{+6\ell _{i}}}&{\color {green}{2\ell _{i}^{2}}}&{\color {green}{-6\ell _{i}}}&{\color {green}{+2\ell _{i}^{2}}}\end{array}}\right)$

Die Beiträge der zwei Elemente sind hier in rot bzw. grün eingefärbt.

Dieses "Hinzuaddieren" passiert in Matlab hier:

                 :
             e = 0;% counter for all elements of all sections
             for section = 1: length(element.noe)
                            :
                 for i = 1:element.noe(section)
                     e = e+1;
                            :
                     j = 2*(e-1)+1;
                            :
                     K(j:j+3,j:j+3) = K(j:j+3,j:j+3)+b*  d(:,:,e)*element.K*d(:,:,e);
                            :
                 end    
                 :

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

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

Einarbeiten der Randbedingungen

Randbedingungen einarbeiten.

Die Randbedingungen arbeiten wir ein, indem wir zeilenweise (für δW und δΦ) und spaltenweise (für W und Φ) streichen.

Hier sind es

• die erste Zeile/Spalte für W₀ (blau) und
• die zweite Zeile / Spalte für Φ₀ (grün).

NONE - see below

Solving

Die Lösung des Anfangs- und Randwertproblems ist in diesem Lösungsansatz mit Matlab^® um die Klasse "FEM_Container" herum aufgebaut - in ihr sind alle Parameter und Lösungsprozesse beschrieben.

Hier heißt die Instanz des Modells

• mathModel

In diesem Container sind alle Parameter und Zustandsgrößen des Modells bereits dimensionslos gemacht - und zwar mit den oben genannten Bezugsgrößen, die in der Excel-Eingabedatei im Blatt "0) Bezugsgrößen" definiert sind.

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:

${\underline {Q}}_{t}(t)={\underline {Q}}_{p}(t)+{\underline {Q}}_{h}(t)$

Particular Solution

Particulare Lösung.

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

${\underline {\underline {K}}}\cdot {\underline {Q}}_{p}={\underline {G}}$

und die ist - mit der Normierung durch w_s  -

 >> Qp
 
 Qp =
 
          0
          0
     0.1055
     0.0042
     0.3542
     0.0064
     0.6680
     0.0072
     1.0000
     0.0073

Diese Lösung tragen wir - elementweise - auf.

Homogene Lösung

Zur Lösung der homogenen Bewegungsgleichung

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

setzen wir an

$\displaystyle {\underline {Q}}_{h}={\underline {\hat {Q}}}_{h}\cdot e^{\lambda \cdot t}$

und erhalten das Eigenwertproblem

$\underbrace {\left(\lambda ^{2}{\underline {\underline {M}}}+{\underline {\underline {K}}}\right)} _{=:\displaystyle {\underline {\underline {D}}}}\cdot {\underline {\hat {Q}}}_{h}={\underline {0}}$

mit

den Eigenwerten	$\lambda$ und
den Egenvektoren	$\displaystyle {\underline {\hat {Q}}}_{h}$.

Für die Berechnung der Eigenwerte λ müssen wir die Abkürzung

$\lambda ^{2}=\Lambda$

einführen. Damit arbeitet die Matlab^®-Routine "eig()", die in "mathModel.eigensystem()" implementiert ist.

Acht Eigenwerte kommen aus der Bedingung $\det({\underline {\underline {D}}})=0$, diese $\Lambda _{i}$ werden auf der Spur der Matrix D abgelegt:

 >> D
 
 D =
 
    1.0e+06 *
 
    -2.9006         0         0         0         0         0         0         0
          0   -1.0774         0         0         0         0         0         0
          0         0   -0.4287         0         0         0         0         0
          0         0         0   -0.1662         0         0         0         0
          0         0         0         0   -0.0480         0         0         0
          0         0         0         0         0   -0.0123         0         0
          0         0         0         0         0         0   -0.0016         0
          0         0         0         0         0         0         0   -0.0000

Die zugehörigen dimensionslosen Periodendauern sind

${\begin{array}{ll}T_{1}=&0.0037\\T_{2}=&0.0061\\T_{3}=&0.0096\\T_{4}=&0.0154\\T_{5}=&0.0287\\T_{6}=&0.0566\\T_{7}=&0.1594\\T_{8}=&1.0000\end{array}}$.

Und offensichtlich fällt die längste Periodendauer T₈ mit der analytisch berechneten untersten Schwingunsperiode T^* zusammen.

Die zugehörigen Eigenvektoren stehen in V:

 >> V
 
 V =
 
    -0.0579    0.2951   -0.0687    0.3043    0.4561    0.5309   -0.3197   -0.0780
    -0.0129    0.0882   -0.0805   -0.0578   -0.0115    0.0064   -0.0096   -0.0032
    -0.0884   -0.1020    0.4873   -0.0417   -0.4707    0.0160   -0.5469   -0.2721
    -0.0364    0.1230    0.0134    0.0659   -0.0007   -0.0222    0.0019   -0.0051
    -0.1652   -0.4100   -0.2420   -0.1979    0.4058   -0.4268   -0.1035   -0.5270
    -0.0763    0.0330    0.0704   -0.0601    0.0098    0.0106    0.0163   -0.0059
    -0.9540   -0.8313   -0.8232    0.9212   -0.6357    0.7308    0.7663   -0.8013
    -0.2099   -0.1390   -0.0999    0.0809   -0.0389    0.0315    0.0200   -0.0060

Die Eigenvektoren in V spannen nun den Nullraum (nullspace) der Matrix auf, es ist

${\underline {\underline {V}}}=\left(\;{\underline {\hat {Q}}}_{h,1}\;,\;{\underline {\hat {Q}}}_{h,2}\;,\;{\underline {\hat {Q}}}_{h,3}\;,\;{\underline {\hat {Q}}}_{h,4}\;,\;{\underline {\hat {Q}}}_{h,5}\;,\;{\underline {\hat {Q}}}_{h,6}\;,\;{\underline {\hat {Q}}}_{h,7}\;,\;{\underline {\hat {Q}}}_{h,8}\;\right)$

Die homogene Lösung der Bewegungsgleichung lautet damit

$\displaystyle {\underline {Q}}_{h}(t)=\sum _{i=1}^{8}C_{i}\cdot {\underline {\hat {Q}}}_{h,i}\cdot e^{\displaystyle j\cdot \omega _{0,i}\cdot \tau }$

mit den acht Integrationskonstanten C_i . Durch die komplexen Eigenwerte sind die Integrationskonstanten nun (eigentlich) auch komplexwertig. Darum kommen in diesem Fall herum, weil die Anfangsgeschwindigkeit des Balkens beim Loslassen Null ist - wir also nur cos-Anteile berücksichtigen müssen. Und die gehören wiederum zum Realteil der Exponential-Funktion.

Die C_i sind nun die Konstanten, die wir brauchen, um die Lösung an Anfangsbedingungen anzupassen.

Vorher schauen wir uns die Lösung jeweils zu einer Eigenfrequenz ω_i an. Diese Funktionen heißen Modalformen ϕ(x) und deren Schwingungen können wir plotten:

Mode	Modalform ϕj(x)	Mode	Modalform ϕj(x)
#8 $\displaystyle \omega _{0,8}=1.000\cdot {\frac {2\pi }{T_{ref}}}$	Mode #1	#7 $\displaystyle \omega _{0,7}=6.2742\cdot {\frac {2\pi }{T_{ref}}}$	Mode #2
#6 $\displaystyle \omega _{0,6}=17.6833\cdot {\frac {2\pi }{T_{ref}}}$	Mode #3	#5 $\displaystyle \omega _{0,5}=34.8854\cdot {\frac {2\pi }{T_{ref}}}$	Mode #4
#4 $\displaystyle \omega _{0,4}=64.8852\cdot {\frac {2\pi }{T_{ref}}}$	Mode #5	#3 $\displaystyle \omega _{0,3}=104.2059\cdot {\frac {2\pi }{T_{ref}}}$	Mode #6

Modalformen

Alle acht Moden ϕ_j können wir auch zum Zeitpunkt τ=0 übereinander darstellen:

✔ Zur Euler-Bernoulli-Hypothese:

Praktisch haben die höheren Moden kaum Relevanz - meist klingen sie durch Dämpfung schnell ab. Man sieht allerdings bereits an der Schwingungs-Form, dass hier die Länge zwischen zwei Knoten nicht mehr sehr klein ist im Vergleich zur Höhe. Und das führt dazu, dass Schubverformungen eine Rolle spielen. Wir müssen zum Timoshenko-Balken wechseln ....

Matlab^®

Main Script

Main - part I (preprocessor)

% DESCRIPTION % _____________________________________ % computes solutions for straight rods % reads parameters from % FEAG.xlsx % writes results to % feag-*.gif % feag-*.png % process % uses own functions % - readSystemParameters(fileName); % - getUnitConversion % - % uses classes % - FEM_Container % - FEM_Element_Modell % - FEM_Section % - FEM_BoundaryCondition % uses MATLAB functions % - linsolve % - eig cd('C:\Users\abs384\OneDrive\Confluence Sources\FEAG'); addpath(['../Matlab/Functions']); addpath(['../Matlab/Classes']); %%**************************************** % PREPROCESSOR * % **************************************** fileName = 'FEAG'; [ bezug, system, boundaries, element, err ] = readSystemParameters(fileName); %% mathModel = FEM_Container(bezug,system,boundaries,element); %%  :  :  :

Main - part II (partikulare Lösung)

:  :  : %% ***************************************************************************/ % p a r t i c u l a r s o l u t i o n */ % ***************************************************************************/ NOV = length(mathModel.free)+length(mathModel.disabled); % initiate particular solution Qp Qp = zeros(NOV,1); Q = mathModel.particularSolution(); Qp(mathModel.free,:)=Q; % get coordinates of nodes L = mathModel.l(length(mathModel.l)); % plot particular solution for ele= 1: mathModel.NOE xi = linspace(mathModel.l(ele)/L,mathModel.l(ele+1)/L,mathModel.samples); plot(xi,Qp(2*ele-1:2*ele+2,1)'*(mathModel.d(:,:,ele)*phi), 'LineWidth',4); hold on; end plot(mathModel.l/mathModel.l(length(mathModel.l)),Qp(1:2:length(Qp)),'ro', 'LineWidth',2); xlabel('\xi\rightarrow'); ylabel('w/w_s\rightarrow'); hold off;

Main - part III (homogene Lösung)

:  :  : %% ***************************************************************************/ % * h o m o g e n e o u s s o l u t i o n */ % ***************************************************************************/ [V,D] =mathModel.eigensystem(); % eigenfrequencies of system omega = sqrt(diag(-D)); Qh = zeros(NOV,length(mathModel.free)); Qh(mathModel.free,:) = V(:,:); % plot modes for mode =length(mathModel.free):-1:1 for ele= 1: mathModel.NOE xi = linspace(mathModel.l(ele)/L,mathModel.l(ele+1)/L,mathModel.samples); plot(xi,Qh(2*ele-1:2*ele+2,mode)'*(mathModel.d(:,:,ele)*phi), 'LineWidth',2); hold on; end plot(mathModel.l/mathModel.l(length(mathModel.l)),Qh(1:2:length(Qp),mode),'ro', 'LineWidth',2); xlabel('\xi\rightarrow'); ylabel('w/w_r\rightarrow'); end hold off;

Adapt to Initial Conditions

Die acht reell-wertigen Integrationskonstanten C_i bestimmen wir aus den Anfangsbedingungen für den Balken, nämlich

${\underline {Q}}_{t}(0)={\underline {0}}{\text{ und }}{\underline {\dot {Q}}}_{t}(0)={\underline {0}}$ ,

also

${\begin{array}{llcll}W_{1}(0)&=0&\;\;\;&{\dot {W}}_{1}(0)&=0\\\Phi _{1}(0)&=0&&{\dot {\Phi }}_{1}(0)&=0\\W_{2}(0)&=0&&{\dot {W}}_{2}(0)&=0\\\Phi _{2}(0)&=0&&{\dot {\Phi }}_{2}(0)&=0\\W_{3}(0)&=0&&{\dot {W}}_{4}(0)&=0\\\Phi _{3}(0)&=0&&{\dot {\Phi }}_{3}(0)&=0\\W_{4}(0)&=0&&{\dot {W}}_{4}(0)&=0\\\Phi _{4}(0)&=0&&{\dot {\Phi }}_{4}(0)&=0\\\end{array}}$.

Wenn wir die C_i komplexwertig denken, so stehen hier 16 Anfangsbedingungen für acht komplexe C_i .

Denken wir die C_i gleich reelwertig, so können wir die acht Anfangsbedingungen in der rechten Spalte weglassen - sie sind dann implizit mit erfüllt - für unser spezielles Anfangswertproblem.

Die Anfangsbedingung für die Gesamtlösung lautet nun:

${\underline {Q}}_{p}+{\underline {Q}}_{h}|_{\tau =0}={\underline {0}}$ ,

also

${\underline {\underline {V}}}\cdot \underbrace {\left({\begin{array}{c}C_{1}\\C_{2}\\C_{3}\\C_{4}\\C_{5}\\C_{6}\\C_{7}\\C_{8}\end{array}}\right)} _{\displaystyle =:{\underline {C}}}=-{\underline {Q}}_{p}$

Wir finden

${\underline {C}}=\left({\begin{array}{l}+6.87897955904839E^{-07}\\-4.64941679695326E^{-06}\\+1.66295172107590E^{-05}\\+4.37871009978632E^{-05}\\+0.000288375976158946\\+0.00145539720136673\\-0.0186576428599944\\-1.26472054045511\end{array}}\right)$ .

und damit

$\displaystyle {\underline {Q}}_{t}(t)={\underline {Q}}_{p}+\sum _{i=1}^{8}C_{i}\cdot {\underline {\hat {Q}}}_{h,i}\cos {\omega _{0,i}\cdot \tau }$ .

Offensichtlich ist der Beitrag der Mode 8 (also der langsamsten Mode) bei weitem am größten Anteil an Q_h. Das können wir auch in der Animation der Lösung - also das Loslassen des Balken aus der Ruhe heraus - sehen.

Matlab^®

Main Script

Main - part IV (adapt to initial conditions)

:  :  : %% ***************************************************************************/ % t o t a l ( c o m p o s e d ) s o l u t i o n */ % ***************************************************************************/ C = linsolve(Qh(mathModel.free,:),Qp(mathModel.free,1)); Qi = zeros(NOV,length(mathModel.free)); Qi(mathModel.free,:) = Qh(mathModel.free,:)*diag(C); for ele= 1: mathModel.NOE xi = linspace(mathModel.l(ele)/L,mathModel.l(ele+1)/L,mathModel.samples); plot(xi,Qp(2*ele-1:2*ele+2,1)'*(mathModel.d(:,:,ele)*phi), 'Color', [0.8 0.8 0.8], 'LineWidth',8); hold on; end for mode =1 : length(mathModel.free) for ele= 1: mathModel.NOE xi = linspace(mathModel.l(ele)/L,mathModel.l(ele+1)/L,mathModel.samples); plot(xi,Qi(2*ele-1:2*ele+2,mode)'*(mathModel.d(:,:,ele)*phi), 'LineWidth',2); hold on; end end hold off;

Post-Processing

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

Lösung im Zeitbereich - hier für I=4.

Animation in Time Domain - hier für I=4.

Matlab^®

Main Script

Main - part IV (adapt to initial conditions)

:  :  : %% animated shape .... tau = linspace(0,3,3*21); h = figure; axis tight manual % this ensures that getframe() returns a consistent size fileName = strcat('feag-initial-value-solution-animation.gif'); for step = tau hold off; for ele= 1: mathModel.NOE xi = linspace(mathModel.l(ele)/L,mathModel.l(ele+1)/L,mathModel.samples); disp = Qp(2*ele-1:2*ele+2,1)'*(mathModel.d(:,:,ele)*phi); for mode =length(mathModel.free):-1:1 disp = disp - Qi(2*ele-1:2*ele+2,mode)'*(mathModel.d(:,:,ele)*phi)*cos(2*pi*omega(mode)/omega(length(omega))*step); end plot(xi,disp, 'LineWidth',2); hold on; axis([0 1 -0.1 2.1]); set(gca,'Ydir','reverse') xlabel('\xi\rightarrow'); ylabel('w/w_r\rightarrow'); end title('initial value solution'); drawnow  % Capture the plot as an image frame = getframe(h); im = frame2im(frame); [imind,cm] = rgb2ind(im,256);  % Write to the GIF File if step == 0 imwrite(imind,cm,fileName,'gif', 'Loopcount',inf,'DelayTime',0); else imwrite(imind,cm,fileName,'gif','WriteMode','append','DelayTime',0); end end hold off; %% time domain .... tau = linspace(0,3,3*71); fileName = 'feag-initial-value-solution-W.gif'; pick = 3:2:mathModel.NOQ; plot(tau, Qp(pick) - Qi(pick,:)*cos(omega(1:8)*tau)); set(gca,'Ydir','reverse') xlabel('\tau\rightarrow'); ylabel('w/w_s\rightarrow'); hold off; %% E - N - D

Links

• ...

Literature

• ...