Gelöste Aufgaben/TC13

Aufgabenstellung

Das System ist eine Variante von Aufgabe TC12. Hier ist eine Näherungslösung mit der Methode der Finiten Elemente gefragt. Das Sytem besteht aus zwei Sektionen mit den Längen ℓ₁ bzw. ℓ₂ sowie den Flächenmomenten I₁ bzw. I₂. Sektion AB ist durch eine konstante Streckenlast q₀ belastet, in B wirkt das Moment M_B0. Der Stab ist in A durch ein gelenkiges Festlager gelagert. In C ist das Stabende fest mit dem Umfang einer Rolle vom Radius r verbunden, die in D frei drehbar gelagert ist. In B sind die beiden Sektionen fest miteinander verbunden. Die Feder in B ist eine Translationsfeder mit der Steifigkeit  k_B.

Lageplan

Gesucht ist die FEM-Lösung mit zwei Elementen für den Euler-Bernoulli-Balken.

Die Systemparameter sind die gleichen wie in TC12, das dort gesuchte Flächenmoment I₁ übernehmen wir zu

$I_1=54000{\text{ mm}}^4$.

Lösung mit Maxima

Header

Wir suchen eine Lösung mit den Standard-Trial-Functions (Hermitesche-Polynome) für einen Euler-Bernoulli-Balken. Aufpassen müssen wir am Rand "C". Hier sind Verschiebung und Verdrehung durch die Rolle gekoppelt. Die Element-Steifigkeitsmatrix müssen wir passend umschreiben.

/*******************************************************/
/* MAXIMA script                                       */
/* version: wxMaxima 15.08.2                           */
/* author: Andreas Baumgart                            */
/* last updated: 2018-04-15                            */
/* ref: TM-C, Labor 4                                  */
/* description: finds the FE solution for              */
/*              lab problem #4                         */
/*******************************************************/

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

Declarations

Wir übernehmen die Systemparameter aus TC12 zu

$\begin{array}{l}{\ell }_1=1000\cdot {m}{m},\\ {\ell }_2=\frac{{\ell }_1}{2},\\ r=\frac{{\ell }_1}{4},\\ E=\frac{210000.0\cdot N}{{{m}{m}}^2},\\ I_1=54000\cdot {{m}{m}}^4,\\ I_2=2\cdot I_1,\\ k_B=\frac{3\cdot I_1\cdot E}{{\ell }_1^3},\\ q_0=\frac{10\cdot N}{{m}{m}},\\ M_B=q_0\cdot {\ell }_1^2\end{array}$,

und wählen als Referenzgröße w_ref für die Auslenkung des Balkens die maximale Verschiebung eines Kragbalkens unter konstanter Streckenlast:

$w_{ref}=\frac{q_o{\ell }_1^4}{8E{I}_1}$.

/*** declarations ***/
assume(l[i]>0);

/* system parameter */
units  : [mm = m/1000, cm = m/100];
params : [l[1]=1000*mm, l[2] = l[1]/2, r=l[1]/4,
          E    = 2.1*10^5*N/mm^2,
          I[1] = 54000*mm^4, I[2]=2*I[1],
          k[B] = 3*E*I[1]/l[1]^3,
          q[0]=10*N/mm, M[B]=q[0]*l[1]^2];

params : subst(units,makelist(lhs(params[i])=subst(params,rhs(params[i])),i,1,length(params)));

/* kinematics of boundary condition at "C" */
kinematics : [W[2]=-r*Φ[2], δW[2] = -r*δΦ[2]];

/* max. displacement of cantilevered beam under load q[0] */
dimless : [w[ref] = q[0]*l[1]^4/(8*E*I[1])];
print(subst(params,dimless))$

Formfunctions

Die Trial-Functions je Element "i" zur Komposition der Form-Funktion kopieren wir aus Finite Elemente Methode zu

$\begin{array}{l}{\varphi }_1=2\cdot {\xi }^3-3\cdot {\xi }^2+1,\\ {\varphi }_2={\ell }_i\cdot ({\xi }^3-2\cdot {\xi }^2+\xi ),\\ {\varphi }_3=3\cdot {\xi }^2-2\cdot {\xi }^3,\\ {\varphi }_4={\ell }_i\cdot ({\xi }^3-{\xi }^2)\end{array}$,

die Koordinaten der Verschiebung sind

$Q_i=\left(\begin{array}{l}{W}_{i-1}\\ {\mathrm{\Phi }}_{i-1}\\ W_i\\ {\mathrm{\Phi }}_i\end{array}\right)$.

Damit ist

$\begin{array}{lcll}{w}_i(x_i)=& & W_{i-1}& \cdot (2\cdot {\xi }^3-3\cdot {\xi }^2+1)\\ & +& {\mathrm{\Phi }}_{i-1}& \cdot {\ell }_i\cdot ({\xi }^3-2\cdot {\xi }^2+\xi )\\ & +& W_i& \cdot (3\cdot {\xi }^2-2\cdot {\xi }^3)\\ & +& {\mathrm{\Phi }}_i& \cdot {\ell }_i\cdot ({\xi }^3-{\xi }^2)\end{array}$.

/**** define form functions ***/
/* coordinates */
coords : [[ W[i-1], Φ[i-1], W[i], Φ[i]],
          [δW[i-1],δΦ[i-1],δW[i],δΦ[i]]];

/* tial-functions */
phi : [ 2*xi^3-3*xi^2+1,
       (xi^3-2*xi^2+xi)*l[i],
        3*xi^2-2*xi^3,
       (xi^3-xi^2)*l[i]];

/* form-function */
form : w(xi) = sum(coords[1][j]*phi[j],j,1,4);

Equilibrium Conditions

Die Gleichgewichtsbedingungen kommen aus dem Prinzip der virtuellen Verrückungen zu

$\begin{array}{ll}\delta W& \stackrel{!}{=}0\\ & =\delta W_a-\delta \mathrm{\Pi }\end{array}$.

Wir sortieren die Elemente der virtuellen Arbeit nach den virtuellen Koordinaten und den Koordinaten des System und konstruieren daraus das gewöhnliche Gleichungssystem

$\underset{\_}{\underset{\_}{K}}\cdot \underset{\_}{Q}=\underset{\_}{P}$.

Die Anteile an der gesamten virtuellen Arbeit kommen aus der virtuellen Formänderungsenergie δΠ und der virtuellen Arbeit der äußeren Lasten δW^a.

Dabei ist für unsere zwei Elemente δΠ = δΠ₁ + δΠ₂' mit

$\delta {\mathrm{\Pi }}_i=\delta {\underset{\_}{Q}}_i^T\cdot {\underset{\_}{\underset{\_}{k}}}_0\cdot \underset{\_}{Q_i}$

und der Element-Steifigkeitsmatrix

${\underset{\_}{\underset{\_}{k}}}_0=\frac{{E}{I}}{{\ell }_i^3}\cdot \left(\begin{array}{cccc}12& 6\cdot {\ell }_i& -12& 6\cdot {\ell }_i\\ 6\cdot {\ell }_i& 4\cdot {\ell }_i^2& -6\cdot {\ell }_i& 2\cdot {\ell }_i^2\\ -12& -6\cdot {\ell }_i& 12& -6\cdot {\ell }_i\\ 6\cdot {\ell }_i& 2\cdot {\ell }_i^2& -6\cdot {\ell }_i& 4\cdot {\ell }_i^2\end{array}\right)$.

Für Element 1 können wir diese Matrix direkt übernehmen, das heißt k₁ = k₀ für i=1.

Element 2 birgt eine Tücke: Hier sind Verschiebung W₂ und Verdrehung Φ₂ gekoppelt über die Kinematik der Rolle:

$W_2=-r\cdot {\mathrm{\Phi }}_2$

Damit ist auch die Variation in diesen Koordinaten gekoppelt, also

$\delta W_2=-r\cdot \delta {\mathrm{\Phi }}_2$.

Mit dieser Beziehung eliminieren wir demnach W₂ aus der Formulierung der Gleichgewichtsbedingungen. Das Einsetzen der kinematischen Beziehungen in δΠ liefert nun als Element-Steifigkeitsmatrix

${\underset{\_}{\underset{\_}{k}}}_2=\frac{E{I}_2}{l_2^3}\cdot \left(\begin{array}{cccc}12& 6\cdot {\ell }_2& 0& 12\cdot r+6\cdot {\ell }_2\\ 6\cdot {\ell }_2& 4\cdot {\ell }_2^2& 0& 6\cdot {\ell }_2\cdot r+2\cdot {\ell }_2^2\\ 0& 0& 0& 0\\ 12\cdot r+6\cdot {\ell }_2& 6\cdot {\ell }_2\cdot r+2\cdot {\ell }_2^2& 0& 12\cdot r^2+12\cdot {\ell }_2\cdot r+4\cdot {\ell }_2^2\end{array}\right)$,

bei der die Elemente der dritten Zeile (zu δW_i) und der dritten Spalte (zu W_i) verschwinden.

Struktur der Steifigkeitsmatrix.

Wir setzen die Gesamt-Steifigkeitsmatrix K aus den beiden Element-Steifigkeitsmatrizen so zusammen:

Und das ist die anteilige Gesamt-Steifigkeitsmatrix der beiden Balken-Abschnitte:

$\underset{\_}{\underset{\_}{K}}=\left(\begin{array}{cccccc}\frac{12\cdot I_1\cdot E}{{\ell }_1^3}& \frac{6\cdot I_1\cdot E}{{\ell }_1^2}& -\frac{12\cdot I_1\cdot E}{{\ell }_1^3}& \frac{6\cdot I_1\cdot E}{{\ell }_1^2}& 0& 0\\ \frac{6\cdot I_1\cdot E}{{\ell }_1^2}& \frac{4\cdot I_1\cdot E}{{\ell }_1}& -\frac{6\cdot I_1\cdot E}{{\ell }_1^2}& \frac{2\cdot I_1\cdot E}{{\ell }_1}& 0& 0\\ -\frac{12\cdot I_1\cdot E}{{\ell }_1^3}& -\frac{6\cdot I_1\cdot E}{{\ell }_1^2}& \frac{12\cdot I_2\cdot E}{{\ell }_2^3}+\frac{12\cdot I_1\cdot E}{{\ell }_1^3}& \frac{6\cdot I_2\cdot E}{{\ell }_2^2}-\frac{6\cdot I_1\cdot E}{{\ell }_1^2}& 0& \frac{12\cdot I_2\cdot r\cdot E}{{\ell }_2^3}+\frac{6\cdot I_2\cdot E}{{\ell }_2^2}\\ \frac{6\cdot I_1\cdot E}{{\ell }_1^2}& \frac{2\cdot I_1\cdot E}{{\ell }_1}& \frac{6\cdot I_2\cdot E}{{\ell }_2^2}-\frac{6\cdot I_1\cdot E}{{\ell }_1^2}& \frac{4\cdot I_2\cdot E}{{\ell }_2}+\frac{4\cdot I_1\cdot E}{{\ell }_1}& 0& \frac{6\cdot I_2\cdot r\cdot E}{{\ell }_2^2}+\frac{2\cdot I_2\cdot E}{{\ell }_2}\\ 0& 0& 0& 0& 0& 0\\ 0& 0& \frac{12\cdot I_2\cdot r\cdot E}{{\ell }_2^3}+\frac{6\cdot I_2\cdot E}{{\ell }_2^2}& \frac{6\cdot I_2\cdot r\cdot E}{{\ell }_2^2}+\frac{2\cdot I_2\cdot E}{{\ell }_2}& 0& \frac{12\cdot I_2\cdot r^2\cdot E}{{\ell }_2^3}+\frac{12\cdot I_2\cdot r\cdot E}{{\ell }_2^2}+\frac{4\cdot I_2\cdot E}{{\ell }_2}\end{array}\right)$

Auch die Feder ist ein elastisches Element und trägt damit eine virtuelle Formänderungsenergie, zu δΠ hinzu, nämlich

$\delta {\mathrm{\Pi }}_3=k_B{W}_1\cdot \delta W_1$

Diesen Beitrag müssen wir in Zeile 3, Spalte 3 hinzuaddieren, also mit der Anweisung

$K_{3,3}=K_{3,3}+k_B$.

Auf der rechten Seite des Gleichungssystem bleiben die Beiträge der virtuellen Arbeiten äußerer Lasten stehen. Dies sind

$\delta W^a=\underset{{\ell }_1}{\int }{q}_0\cdot \delta w(x)dx+M_B\cdot \delta {\mathrm{\Phi }}_1$,

wir erhalten

$\underset{\_}{P}=\left(\begin{array}{c}\frac{q_0\cdot {\ell }_1}{2}\\ \frac{q_0\cdot {\ell }_1^2}{12}\\ \frac{q_0\cdot {\ell }_1}{2}\\ M_B-\frac{q_0\cdot {\ell }_1^2}{12}\\ 0\\ 0\end{array}\right)$.

Jetzt fehlen nur noch die Randbedingungen. Wir setzen

$W_0=0$,

die kinematische Beziehung

$W_2=-r\cdot {\mathrm{\Phi }}_2$

haben wir schon eingearbeit - bleibt also nur noch, die entsprechenden Spalten und Zeilen in den Matrizen K, Q und P zu eliminieren. Es bleibt das Gleichungssystem

$\left(\begin{array}{cccc}\frac{4\cdot I_1\cdot E}{{\ell }_1}& -\frac{6\cdot I_1\cdot E}{{\ell }_1^2}& \frac{2\cdot I_1\cdot E}{{\ell }_1}& 0\\ -\frac{6\cdot I_1\cdot E}{{\ell }_1^2}& \frac{12\cdot I_2\cdot E}{{\ell }_2^3}+\frac{12\cdot I_1\cdot E}{{\ell }_1^3}+k_B& \frac{6\cdot I_2\cdot E}{{\ell }_2^2}-\frac{6\cdot I_1\cdot E}{{\ell }_1^2}& \frac{12\cdot I_2\cdot r\cdot E}{{\ell }_2^3}+\frac{6\cdot I_2\cdot E}{{\ell }_2^2}\\ \frac{2\cdot I_1\cdot E}{{\ell }_1}& \frac{6\cdot I_2\cdot E}{{\ell }_2^2}-\frac{6\cdot I_1\cdot E}{{\ell }_1^2}& \frac{4\cdot I_2\cdot E}{{\ell }_2}+\frac{4\cdot I_1\cdot E}{{\ell }_1}& \frac{6\cdot I_2\cdot r\cdot E}{{\ell }_2^2}+\frac{2\cdot I_2\cdot E}{{\ell }_2}\\ 0& \frac{12\cdot I_2\cdot r\cdot E}{{\ell }_2^3}+\frac{6\cdot I_2\cdot E}{{\ell }_2^2}& \frac{6\cdot I_2\cdot r\cdot E}{{\ell }_2^2}+\frac{2\cdot I_2\cdot E}{{\ell }_2}& \frac{12\cdot I_2\cdot r^2\cdot E}{{\ell }_2^3}+\frac{12\cdot I_2\cdot r\cdot E}{{\ell }_2^2}+\frac{4\cdot I_2\cdot E}{{\ell }_2}\end{array}\right)\cdot \left(\begin{array}{c}{\mathrm{\Phi }}_0\\ W_1\\ {\mathrm{\Phi }}_1\\ {\mathrm{\Phi }}_2\end{array}\right)=\left(\begin{array}{c}\frac{q_0\cdot {\ell }_1^2}{12}\\ \frac{q_0\cdot {\ell }_1}{2}\\ M_B-\frac{q_0\cdot {\ell }_1^2}{12}\\ 0\end{array}\right)$.

/**** equilibrium conditions ***/
/* generic stiffness matrix */
k[0] : E*I[i]/l[i]^3*makelist(makelist(
              integrate(
                      diff(phi[i],xi,2)*diff(phi[j],xi,2),
                                              xi,0,1),
                                               j,1,4),i,1,4);
k[0] : funmake('matrix,k[0]);

/* element stiffness matrices */
/* element 1 */
k[1] : subst([i=1],k[0]);
/* element 2 */
δΠ[2] : subst([i=2],coords[2].k[0].transpose(coords[1]));
δΠ[2] : expand(subst(kinematics, δΠ[2]));
k[2]  : funmake('matrix,makelist(makelist(coeff(coeff(δΠ[2],subst([i=2],coords[2])[m]),subst([i=2],coords[1])[n]),m,1,4),n,1,4));




/* compose system matrix */
NoN : 3; /* Number of Nodes*/
K : zeromatrix(2*NoN,2*NoN);

for m:1 thru 4 do
   for n:1 thru 4 do
      (K[  m,  n] : K[  m,  n] + k[1][m,n],
       K[2+m,2+n] : K[2+m,2+n] + k[2][m,n])$

/* add spring */
K[3,3] : K[3,3] + k[B]; /* W[1] */

/* compose righ-hand-side */
P : zeromatrix(6,1);
for m:1 thru 4 do
    P[m,1] : l[1]*integrate(q[0]*subst([i=1],phi[m]), xi,0,1);
P[4,1] : P[4,1] + M[B];


/* coordinates of displacement */
Q : matrix([W[0]],[Φ[0]],[W[1]],[Φ[1]],[W[2]],[Φ[2]]);


/* incorporate geometric boundary conditions */
/* eliminate rows / columns for W[0], Φ[2] (positions 1, 6) */
K : submatrix(1,submatrix(5,K,5),1);
Q : submatrix(1,submatrix(5,Q));
P : submatrix(1,submatrix(5,P));

print("k[1] = ", k[1])$
print("k[2] = ", k[2])$
print("K = ", K)$
print("Q = ", Q)$
print("P = ", P)$

print(K," * ",Q," = ",P)$

Solving

Auflösen nach Q liefert

$\begin{array}{ll}{\mathrm{\Phi }}_0& =-0.0157,\\ W_1& =-4.86\cdot {10}^{-18}\text{m},\\ {\mathrm{\Phi }}_1& =+0.0682,\\ {\mathrm{\Phi }}_2& =-0.0262\end{array}$,

zusammen mit den Randbedingungen also

$\begin{array}{ll}{W}_0& =0,\\ {\mathrm{\Phi }}_0& =-0.0157,\\ W_1& =-4.86\cdot {10}^{-18}\text{m},\\ {\mathrm{\Phi }}_1& =+0.0682,\\ W_2& =+0.00656\text{ m},\\ {\mathrm{\Phi }}_2& =-0.0262\end{array}$.

/**** solve K Q = P ***/

sol[1] : linsolve_by_lu(subst(params,K),subst(params,P))[1]$
sol[1] : makelist(Q[i][1] = sol[1][i][1],i,1,4);
sol[1] : append(subst(params,[W[0]=0, W[2]=-r*subst(sol[1], Φ[2])]), sol[1]);

Post-Processing

Biegelinie w(x)

Die Auslenkung w(x) der Querschnitte tragen wir jetzt auf. Für die dimensionslose Darstellung wählen wir aus der analytischen Lösung für den Kragbalken unter konstanter Streckenlast

$\begin{array}{ll}{w}_{{r}{e}{f}}& =\frac{q_0\cdot {\ell }_1^4}{8\cdot E{I}_1}\\ & =0.11\text{ m}\end{array}$

und tragen w(x)/w_ref auf:

Die maximale Auslenkung des Systems ist dann

$\begin{array}{ll}{w}_{max}& \approx 0.1\cdot 0.11\text{ m}\\ & \approx 11\text{ mm}\end{array}$,

und das entspricht der Vorgabe aus Aufgabe TC12.

/**** post-process  ***/
f1 : subst([xi=t],subst(params,subst(dimless,expand(subst(sol[1],subst([i=1],subst(form,w(xi)/w[ref])))))));
f2 : subst([xi=t],subst(params,subst(dimless,expand(subst(sol[1],subst([i=2],subst(form,w(xi)/w[ref])))))));
scale : subst(params,(l[2]/l[1]));

plot2d([[parametric,   t,       f1, [t, 0, 1]],
        [parametric, 1+t*scale, f2, [t, 0, 1]]],
       [gnuplot_preamble, "set yrange [] reverse"],
       [legend, "sec. I", "sec. II"],
       [xlabel, "x/l[1] →"], [ylabel, "w/w[ref] →"])$

---

Links

• ...

Literature

• ...