Gelöste Aufgaben/Kw96

Aufgabenstellung

Ein Stab ABC ist durch eine lineare veränderliche Streckenlast q mit den Eckwerten q_A in A und q_B in B sowie dem Moment M_B in B belastet. Der Stab (E-Modul: E) besteht aus zwei Sektionen mit den Längen l₁ bzw. l₂ sowie den Flächenmomenten I₁ bzw. I₂. Der Stab ist in A durch ein gelenkiges Festlager, in C durch eine Schiebehülse gelagert, in B sind die beiden Sektionen fest miteinander verbunden. Die Feder in A ist eine Drehfester mit Steifigkeit K_A, die Federn in B und C sind Translationsfedern mit den Steifigkeiten k_B, k_C.

Lageplan

Gesucht ist die FEM Lösung für den Euler-Bernoulli-Balken unter Verwendung von zwei Finiten Elementen.

Systemparameter

Ermitteln Sie für ein Euler-Bernoulli-Modell die analytischen Verläufe der Schnittgrößen und Verschiebungen im Balken für die angegebenen Parameter:

Die analytische Lösung finden Sie in Aufgabe Kw98.

Lösung mit Maxima

Header

Wir arbeiten mit den Standard-System-Matrizen nach Abschnitt "FEM-Formulierung für den Euler-Bernoulli-Balken".

📑 toggle code

/*******************************************************/
/* MAXIMA script                                       */
/* version: wxMaxima 18.10.1                           */
/* author: Andreas Baumgart                            */
/* last updated: 2019-02-10                            */
/* ref: TM-C, Labor 4                                  */
/* description: finds the FE solution for              */
/*              lab problem #4                         */
/*******************************************************/

Declarations

System-Parameter sind:

$\begin{array}{l}{q}_A=3000\frac{N}{m},\\ {\ell }_1=\frac{7\cdot m}{10},\\ {{E}{I}}_1=33600N{m}^2,\\ {\ell }_2=\frac{21}{40}m,\\ {{E}{I}}_2=16800N{m}^2,\\ K_A=96000Nm,\\ k_C=\frac{256}{229}\cdot k_B,\\ k_B=\frac{256}{229}Nm,\\ q_B=12000\frac{N}{m},\\ M_B=1470Nm\end{array}$.

📑 toggle code

/* declare variational variables - see 6.3 Identifiers */
declare("δA", alphabetic);
declare("δΠ", alphabetic);
declare("δW", alphabetic);
declare("δΦ", alphabetic);
declare("δw", alphabetic);
declare( "ℓ", alphabetic);
 
/* declarations */
assume(ℓ[i]>0);
 
/* system parameter */
units  : [mm = m/1000, cm = m/100];
params : [q[A]=3*N/mm, ℓ[1]=700*mm, EI[1] = 2.1*10^11*N/m^2 * 3*cm*(4*cm)^3/12];
simple : [ℓ[2] = 3/4*ℓ[1], EI[2] = EI[1]/2,
          K[A] = 2*EI[1]/ℓ[1], K[C] = 512/229*K[B]/2, K[B] = 2*EI[1]/ℓ[1]^3,
          q[B] = 4*q[A], M[B] = q[A]*ℓ[1]^2];
 
params : append(params,makelist(
               lhs(simple[i])=subst(params,rhs(simple[i])),i,1,length(simple)));
params : subst(units,params);

Formfunctions

Die Ansatzfunktion für die Trial-Functions ist ein Polynom 3. Grades:

$\mathrm{w}\left(\xi \right):=c_3\cdot {\xi }^3+c_2\cdot {\xi }^2+c_1\cdot \xi +c_0$

An den Rändern müssen die Auslenkung und Kippung mit den Knoten-Variablen übereinstimmen:

$\begin{array}{l}{c}_0=W_{i-1},\\ \frac{c_1}{l_i}={\mathrm{\Phi }}_{i-1},\\ c_3+c_2+c_1+c_0=W_i,\\ \frac{c_1+2\cdot c_2+3\cdot c_3}{l_i}={\mathrm{\Phi }}_i\end{array}$

Trial-Functions

Damit ist die Ansatzfunktion des Finiten Elements mit den vier Knotenvariablen

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

📑 toggle code

/**** define form functions ***/
/* coordinates */
coords : [[ W[i-1], Φ[i-1], W[i], Φ[i]],
          [δW[i-1],δΦ[i-1],δW[i],δΦ[i]]];
 
/* generic polynomials */
define(w(xi),sum(c[i]*xi^(i),i,0,3));
/* … and boundary conditions */
bc : [w(0) = W[i-1],
      subst([xi=0],diff(w(xi),xi)/ℓ[i])=Φ[i-1],
      w(1) = W[i],
      subst([xi=1],diff(w(xi),xi)/ℓ[i])=Φ[i]];
 
/* solve for c's to comply with geometric boundary conditions      */
coeffs : solve(bc, makelist(c[i],i,0,3))[1];  
print("w(ξ) = ",facsum(expand(subst(coeffs,w(xi))),coords[1]));
trialfcts: makelist(phi[i]=coeff(expand(subst(coeffs,w(xi))),coords[1][i]),i,1,4);
 
plot2d(subst([ℓ[i]=5],subst(trialfcts,makelist(phi[i],i,1,4))),
                    [xi,0,1],[xlabel,"ξ/1 →"],
                    [legend,"φ1","5*φ2/ℓ","φ3","5*φ4/ℓ"]);

Equilibrium Conditions

So sind die Element-Steifigkeitsmatrix

${\underset{\_}{\underset{\_}{K}}}_i=\frac{E{I}_i}{{\ell }_i^3}\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)$

die Koordinaten des FE-Modells - hier für das Element "1":

$\underset{\_}{Q}=\left(\begin{array}{c}{\mathrm{\Phi }}_0\\ W_1\\ {\mathrm{\Phi }}_1\\ W_2\end{array}\right)$.

Wir komponieren daraus die System-Steifigkeitsmatrix - durch Aufaddieren der Beiträge der beiden Elemente und Einarbeiten der Randbedingugnen - zu

${\underset{\_}{\underset{\_}{K}}}_0=\left(\begin{array}{cccc}{K}_A+\frac{4{{E}{I}}_1}{{\ell }_1}& -\frac{6{{E}{I}}_1}{{\ell }_1^2}& \frac{2{{E}{I}}_1}{{\ell }_1}& 0\\ -\frac{6{{E}{I}}_1}{{\ell }_1^2}& k_B+\frac{12{{E}{I}}_2}{{\ell }_2^3}+\frac{12{{E}{I}}_1}{{\ell }_1^3}& \frac{6{{E}{I}}_2}{{\ell }_2^2}-\frac{6{{E}{I}}_1}{{\ell }_1^2}& -\frac{12{{E}{I}}_2}{{\ell }_2^3}\\ \frac{2{{E}{I}}_1}{{\ell }_1}& \frac{6{{E}{I}}_2}{{\ell }_2^2}-\frac{6{{E}{I}}_1}{{\ell }_1^2}& \frac{4{{E}{I}}_2}{{\ell }_2}+\frac{4{{E}{I}}_1}{{\ell }_1}& -\frac{6{{E}{I}}_2}{{\ell }_2^2}\\ 0& -\frac{12{{E}{I}}_2}{{\ell }_2^3}& -\frac{6{{E}{I}}_2}{{\ell }_2^2}& k_C+\frac{12{{E}{I}}_2}{{\ell }_2^3}\end{array}\right)$

Wie das geht, steht in Abschnitt Finite Elemente Methode.

📑 toggle code

/********************************************/
/* EUILIBRIUM CONDITIONS                    */
/* generic stiffness matrix */
K[i] : EI[i]/ℓ[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[i] : subst(trialfcts,K[i]);
K[i] : ev(K[i],nouns);
K[i] : funmake('matrix,K[i]);
 
/* compose system matrix */
NoN : 3; /* Number of Nodes*/
K[0] : zeromatrix(2*NoN,2*NoN);
 
for m:1 thru 4 do
   for n:1 thru 4 do
       (K[0][  m,  n] : K[0][  m,  n] + subst([i=1],K[i][m,n]),
        K[0][2+m,2+n] : K[0][2+m,2+n] + subst([i=2],K[i][m,n]));
/* add springs */
K[0][2,2] : K[0][2,2] + K[A]; /* Φ[0] */
K[0][3,3] : K[0][3,3] + K[B]; /* W[1] */
K[0][5,5] : K[0][5,5] + K[C]; /* W[2] */
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[0] : submatrix(1,submatrix(6,K[0],6),1);
Q : submatrix(1,submatrix(6,Q));
 
/* compose righ-hand-side */

P : transpose(funmake('matrix,
  [subst([i=1],append(makelist(integrate(
            ℓ[i]*subst(trialfcts,(q[A]*(1-xi)+q[B]*xi)*phi[j]),
            xi,0,1),j,1,4),[0,0]))]
      ));
P[4,1] : P[4,1]+M[B];

/* eliminate BCs */
P : submatrix(1,6,P);
 
print('K[i]," = ",EI[i]/ℓ[i]^3, ratsimp(K[i]/(EI[i]/ℓ[i]^3)))$
print('K[0]," = ", K[0])$
print("Q = ", Q)$
print("P = ", P)$

Solving

Die Knotenvariablen sind damit

$\begin{array}{ll}{W}_0=0& \\ {\mathrm{\Phi }}_0=0.00624& \\ W_1=0.00657m& \\ {\mathrm{\Phi }}_1=0.0123& \\ W_2=0.00846m& \\ {\mathrm{\Phi }}_2=0& \end{array}$.

📑 toggle code

/********************************************/
/* SOLVE                                    */
fpprintprec: 3;
sol[1] : linsolve_by_lu(subst(params,K[0]),subst(params,P))[1];
sol[1] : append([W[0]=0],makelist(Q[i][1] = sol[1][i][1],i,1,4), [Φ[2]=0]);

Post-Processing

Biegelinie w(x)

Die Biegelinie des Balkens sieht damit so aus:

📑 toggle code

/********************************************/
/* POSTPROCESS                              */
/* user parametric plot- hence t as independant variable */
f1 : subst([xi=t],subst(params,expand(subst(sol[1],subst([i=1],subst(coeffs,w(xi)))))));
f2 : subst([xi=t],subst(params,expand(subst(sol[1],subst([i=2],subst(coeffs,w(xi)))))));
scale : subst(params,(ℓ[2]/ℓ[1]));
 
plot2d([[parametric,   t,       f1/m, [t, 0, 1]],
            [parametric, 1+t*scale, f2/m, [t, 0, 1]]],
            [gnuplot_preamble, "set yrange [] reverse"],
            [legend, "sec. I", "sec. II"],
                             [xlabel, "x/ℓ[1] →"],
                             [ylabel, "← w/m"])$

---

Links

• Analytische Lösung des Problems: Kw98

Literature

• ...