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.

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. Allerdings können wir Matlab-Inhalte nicht gut auf dieser Seite unterbringen - deshalb gibt es dafür die Seite FEAG-Matlab, die der gleichen Struktur folgt.

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 ω₀ 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 ω₀ 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 α₁ 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)

.

C O D E

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

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).

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

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 -

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

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

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.

Post-Processing

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

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

Links

...

Literature

...

Header

Matlab

Classes

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

Class: 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

Class: 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

Class: 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

Class: FEM-Section

Gelöste Aufgaben/FEAG

Inhaltsverzeichnis

Aufgabenstellung

Lösung mit Matlab

Header

Bezugsgrößen wählen

System-Parameter des FEM-Modells

Equilibrium Conditions

Komponieren der System-Matrizen

Einarbeiten der Randbedingungen

Solving

Particular Solution

Homogene Lösung

Adapt to Initial Conditions

Post-Processing

Header

Matlab

Classes

Eingabe-Parameter aus Excel

Navigationsmenü

Gelöste Aufgaben/FEAG

Aufgabenstellung

Lösung mit Matlab

Header

Bezugsgrößen wählen

System-Parameter des FEM-Modells

Equilibrium Conditions

Komponieren der System-Matrizen

Einarbeiten der Randbedingungen

Solving

Particular Solution

Homogene Lösung

Adapt to Initial Conditions

Post-Processing

Header

Matlab

Classes

Eingabe-Parameter aus Excel

Navigationsmenü

Suche