Aufgabenstellung
Ein Stab ABC ist durch eine lineare veränderliche Streckenlast q mit den Eckwerten qA in A und qB in B sowie dem Moment MB in B belastet. Der Stab (E-Modul: E ) besteht aus zwei Sektionen mit den Längen l1 bzw. l2 sowie den Flächenmomenten I1 bzw. I2 . 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 KA , die Federn in B und C sind Translationsfedern mit den Steifigkeiten kB , kC .
Lageplan
Gesucht ist die analytische Lösung für den Euler-Bernoulli-Balken.
Systemparameter
Ermitteln Sie für ein Euler-Bernoulli-Modell die analytischen Verläufe der Schnittgrößen und Verschiebungen im Balken für diese Parameter:
Lösung mit Maxima
Die Aufgabe ist ein klassisches Randwertproblem:
zwei Gebiete, in denen ein Euler-Bernoulli-Balken in AB und BC durch eine Streckenlast q belastet ist (in Bereich II ist die Streckenlast allerdings Null) und somit durch die DifferentialbeziehungE I i w i I V ( x i ) = q ( x i ) , i = { 1 , 2 } berschrieben wird.
Rand- und Übergangsbedingungen in den Punkten A, B, C
Wir verwenden xi und ξi als Koordinaten je Bereich, in der Übersicht sieht das Randwertproblem so aus:
Rand A Bereich I Übergang B Bereich II Rand C
Diese Aufgabe mit der Methode der Finiten Elemente in KW96 gelöst.
📑 toggle code
/*******************************************************/
/* MAXIMA script */
/* version: wxMaxima 15.08.2 */
/* author: Andreas Baumgart */
/* last updated: 2017-09-06 */
/* ref: TM-C, Labor 1 */
/* description */
/* */
/*******************************************************/
Declarations
Wir definieren die Parameter
q A = 3 ⋅ k N m , l 1 = 7 ⋅ m 1 0 , E I 1 = 3 3 6 0 0 N m 2 , l 2 = 2 1 m 4 0 , E I 2 = 1 6 8 0 0 N m 2 , K A = 9 6 k N m , k C = 2 2 k N m , k B = 9 8 k N m , q B = 1 2 N m m , M B = 1 4 7 0 N m .
und die Formfunktionen für die Streckenlast
ϕ 0 ( ξ ) : = 1 − ξ ϕ 1 ( ξ ) : = ξ .
📑 toggle code
/* system parameter */
units : [mm = m/1000, cm = m/100] ;
params : [q[A]=3*N/mm, l[1]=700*mm, EI[1] = 2.1*10^11*N/m^2 * 3*cm* ( 4*cm ) ^3/12] ;
simple : [l[2] = 3/4*l[1], EI[2] = EI[1]/2,
K[A]=2*EI[1]/l[1], k[C] = 512/229*EI[1]/l[1]^3, k[B] = EI[1]/l[1]^3,
q[B] = 4*q[A], M[B] = q[A]*l[1]^2] ;
params : append ( params,makelist ( lhs ( simple[i] ) =subst ( params,rhs ( simple[i] )) , i,1,length ( simple ))) ;
params : subst ( units,params ) ;
/* form - functions */
phi[0] ( xi ) := 1 - xi ;
phi[1] ( xi ) := xi ;
Formfunctions
In Bereich I und II gilt dieselbe Bewegungs-Differentialgleichung
E I i w i I V ( x i ) = q ( x i ) , i = { 1 , 2 } mit q ( x i ) = q 0 ⋅ ϕ 0 ( ξ ) + q 1 ⋅ ϕ 1 ( ξ ) ,
die wir durch Integration lösen und dann bereichsweise anpassen.
So gilt für Bereich II: q0 = 0 und q1 = 0 .
Die allgemeine Lösung ist mit
ϕ i ( x ) = d w ( x ) d x
... für Bereich I:
w 1 ( x ) : = 1 2 0 ⋅ l 1 ⋅ C 1 , 0 + 1 2 0 ⋅ l 1 ⋅ C 1 , 1 ⋅ x + 6 0 ⋅ l 1 ⋅ C 1 , 2 ⋅ x 2 + 2 0 ⋅ l 1 ⋅ C 1 , 3 ⋅ x 3 + 5 ⋅ l 1 ⋅ x 4 ⋅ q A + x 5 ⋅ ( q B − q A ) 1 2 0 ⋅ l 1 ⋅ E I 1 ϕ 1 ( x ) : = 1 2 0 ⋅ l 1 ⋅ C 1 , 1 + 1 2 0 ⋅ l 1 ⋅ C 1 , 2 ⋅ x + 6 0 ⋅ l 1 ⋅ C 1 , 3 ⋅ x 2 + 2 0 ⋅ l 1 ⋅ x 3 ⋅ q A + 5 ⋅ x 4 ⋅ ( q B − q A ) 1 2 0 ⋅ l 1 ⋅ E I 1 M 1 ( x ) : = − 1 2 0 ⋅ l 1 ⋅ C 1 , 2 + 1 2 0 ⋅ l 1 ⋅ C 1 , 3 ⋅ x + 6 0 ⋅ l 1 ⋅ x 2 ⋅ q A + 2 0 ⋅ x 3 ⋅ ( q B − q A ) 1 2 0 ⋅ l 1 Q 1 ( x ) : = − 1 2 0 ⋅ l 1 ⋅ C 1 , 3 + 1 2 0 ⋅ l 1 ⋅ x ⋅ q A + 6 0 ⋅ x 2 ⋅ ( q B − q A ) 1 2 0 ⋅ l 1
... für Bereich II:
w 2 ( x ) : = 1 2 0 ⋅ l 2 ⋅ C 2 , 0 + 1 2 0 ⋅ l 2 ⋅ C 2 , 1 ⋅ x + 6 0 ⋅ l 2 ⋅ C 2 , 2 ⋅ x 2 + 2 0 ⋅ l 2 ⋅ C 2 , 3 ⋅ x 3 1 2 0 ⋅ l 2 ⋅ E I 2 ϕ 2 ( x ) : = 1 2 0 ⋅ l 2 ⋅ C 2 , 1 + 1 2 0 ⋅ l 2 ⋅ C 2 , 2 ⋅ x + 6 0 ⋅ l 2 ⋅ C 2 , 3 ⋅ x 2 1 2 0 ⋅ l 2 ⋅ E I 2 M 2 ( x ) : = − 1 2 0 ⋅ l 2 ⋅ C 2 , 2 + 1 2 0 ⋅ l 2 ⋅ C 2 , 3 ⋅ x 1 2 0 ⋅ l 2 Q 2 ( x ) : = − C 2 , 3 .
📑 toggle code
/* solve .... */
dgl : EI[i]*diff ( w ( x ) , x,4 ) = q[0]*phi[0] ( x/l[i] ) + q[1]*phi[1] ( x/l[i] ) ;
/* generic solution */
displ : solve ( integrate ( integrate ( integrate ( integrate ( dgl,x ) , x ) , x ) , x ) , w ( x )) ;
sections: [[i=1, %c4=C[1,0], %c3=C[1,1], %c2=C[1,2], %c1=C[1,3], q[0]=q[A], q[1]=q[B]],
[i=2, %c4=C[2,0], %c3=C[2,1], %c2=C[2,2], %c1=C[2,3], q[0]= 0 , q[1]= 0 ]] ;
/* section I */
define ( w[1] ( x ) , subst ( sections[1],subst ( displ,w ( x )))) ;
define ( Phi[1] ( x ) , diff ( w[1] ( x ) , x )) ;
define ( M[1] ( x ) , -EI[1]*diff ( w[1] ( x ) , x,2 )) ;
define ( Q[1] ( x ) , -EI[1]*diff ( w[1] ( x ) , x,3 )) ;
/* section II */
define ( w[2] ( x ) , subst ( sections[2],subst ( displ,w ( x )))) ;
define ( Phi[2] ( x ) , diff ( w[2] ( x ) , x )) ;
define ( M[2] ( x ) , -EI[2]*diff ( w[2] ( x ) , x,2 )) ;
define ( Q[2] ( x ) , -EI[2]*diff ( w[2] ( x ) , x,3 )) ;
Boundary Conditions
Für die 2*4 = 8 Integrationskonstanten
[ C 1 , 0 , C 1 , 1 , C 1 , 2 , C 1 , 3 , C 2 , 0 , C 2 , 1 , C 2 , 2 , C 2 , 3 ]
suchen wir jetzt die passenden Gleichungen aus Rand- und Übergangsbedingungen.
Zur besseren Übersicht nennen wir die Schnitt-Momente und -Kräfte nach den jeweiligen Knotenpunkten A, B, C und fügen als Index ein + / - hinzu, um die Seite (+: rechts vom Knoten, -: links vom Knoten) zu kennzeichnen.
Aus Rand "A"
Geometrische Randbedingungen
w 1 ( 0 ) = 0
Kraft- und Momenten-Randbedingungen
K A ⋅ ϕ A + M A , + = 0 mit M A , + = − E I 1 ⋅ w ″ ( x ) | x = 0
Aus Übergang "B"
Geometrische Randbedingungen
w 1 ( ℓ 1 ) = w 2 ( ℓ 1 )
ϕ 1 ( ℓ 1 ) = ϕ 2 ( ℓ 1 )
Kraft- und Momenten-Randbedingungen
− M B , − − M B + M B , + = 0
− Q B , − − k B ⋅ w B + − Q B , + = 0
Aus Rand "C"
Geometrische Randbedingungen
ϕ 2 ( ℓ 2 ) = 0
Kraft- und Momenten-Randbedingungen
− Q C , − − k C ⋅ w C = 0
Und das liefert das Gleichungssystem aus 8 Gleichungen
( C 1 , 0 E I 1 = 0 C 1 , 1 ⋅ K A E I 1 − C 1 , 2 = 0 ℓ 1 4 ⋅ q B 1 2 0 ⋅ E I 1 + ℓ 1 4 ⋅ q A 3 0 ⋅ E I 1 + ℓ 1 3 ⋅ C 1 , 3 6 ⋅ E I 1 + ℓ 1 2 ⋅ C 1 , 2 2 ⋅ E I 1 + ℓ 1 ⋅ C 1 , 1 E I 1 + C 1 , 0 E I 1 = C 2 , 0 E I 2 ℓ 1 3 ⋅ q B 2 4 ⋅ E I 1 + ℓ 1 3 ⋅ q A 8 ⋅ E I 1 + ℓ 1 2 ⋅ C 1 , 3 2 ⋅ E I 1 + ℓ 1 ⋅ C 1 , 2 E I 1 + C 1 , 1 E I 1 = C 2 , 1 E I 2 ℓ 1 ⋅ q B 2 − C 2 , 0 ⋅ k B E I 2 + ℓ 1 ⋅ q A 2 − C 2 , 3 + C 1 , 3 = 0 − M B + ℓ 1 2 ⋅ q B 6 + ℓ 1 2 ⋅ q A 3 − C 2 , 2 + ℓ 1 ⋅ C 1 , 3 + C 1 , 2 = 0 ℓ 2 2 ⋅ C 2 , 3 2 ⋅ E I 2 + ℓ 2 ⋅ C 2 , 2 E I 2 + C 2 , 1 E I 2 = 0 − ℓ 2 3 ⋅ C 2 , 3 ⋅ k C 6 ⋅ E I 2 − ℓ 2 2 ⋅ C 2 , 2 ⋅ k C 2 ⋅ E I 2 − ℓ 2 ⋅ C 2 , 1 ⋅ k C E I 2 − C 2 , 0 ⋅ k C E I 2 + C 2 , 3 = 0 ) .
für die Integrationskonstanten.
📑 toggle code
/* integration constants */
ICs : [C[1,0],C[1,1],C[1,2],C[1,3],C[2,0],C[2,1],C[2,2],C[2,3]] ;
/* boundary conditions */
node[A]: [ w[1] ( 0 ) = 0 ,
K[A]*Phi[1] ( 0 ) +M[1] ( 0 ) = 0] ;
node[B]: [ w[1] ( l[1] ) = w[2] ( 0 ) ,
Phi[1] ( l[1] ) = Phi[2] ( 0 ) ,
-Q[1] ( l[1] ) -k[B]*w[2] ( 0 ) +Q[2] ( 0 ) = 0 ,
-M[1] ( l[1] ) -M[B]+M[2] ( 0 ) = 0] ;
node[C]: [ Phi[2] ( l[2] ) = 0 ,
-Q[2] ( l[2] ) - k[C]*w[2] ( l[2] ) = 0] ;
BCs : expand ( append ( node[A],node[B],node[C] )) ;
Prepare for Solver
Das Gleichungssystem wollen wir als
A _ _ ⋅ x _ = b _
schreiben, also
( 1 E I 1 0 0 0 0 0 0 0 0 K A E I 1 − 1 0 0 0 0 0 1 E I 1 ℓ 1 E I 1 ℓ 1 2 2 ⋅ E I 1 ℓ 1 3 6 ⋅ E I 1 − 1 E I 2 0 0 0 0 1 E I 1 ℓ 1 E I 1 ℓ 1 2 2 ⋅ E I 1 0 − 1 E I 2 0 0 0 0 0 1 − k B E I 2 0 0 − 1 0 0 1 ℓ 1 0 0 − 1 0 0 0 0 0 0 1 E I 2 ℓ 2 E I 2 ℓ 2 2 2 ⋅ E I 2 0 0 0 0 − k C E I 2 − ℓ 2 ⋅ k C E I 2 − ℓ 2 2 ⋅ k C 2 ⋅ E I 2 − ℓ 2 3 ⋅ k C − 6 ⋅ E I 2 6 ⋅ E I 2 ) ⋅ x _ = ( 0 0 − ℓ 1 4 ⋅ q B 1 2 0 ⋅ E I 1 − ℓ 1 4 ⋅ q A 3 0 ⋅ E I 1 − ℓ 1 3 ⋅ q B 2 4 ⋅ E I 1 − ℓ 1 3 ⋅ q A 8 ⋅ E I 1 − ℓ 1 ⋅ q B 2 − ℓ 1 ⋅ q A 2 M B − ℓ 1 2 ⋅ q B 6 − ℓ 1 2 ⋅ q A 3 0 0 )
Die Matrix-Elemente sind für die Koeffizientenmatrix
a 1 , 1 = 1 / E I 1 a 2 , 2 = K A / E I 1 a 2 , 3 = − 1 a 3 , 1 = 1 / E I 1 a 3 , 2 = ℓ 1 / E I 1 a 3 , 3 = ℓ 1 2 / ( 2 ⋅ E I 1 ) a 3 , 4 = ℓ 1 3 / ( 6 ⋅ E I 1 ) a 3 , 5 = − 1 / E I 2 a 4 , 2 = 1 / E I 1 a 4 , 3 = ℓ 1 / E I 1 a 4 , 4 = ℓ 1 2 / ( 2 ⋅ E I 1 ) a 4 , 6 = − 1 / E I 2 a 5 , 4 = 1 a 5 , 5 = − k B / E I 2 a 5 , 8 = − 1 a 6 , 3 = 1 a 6 , 4 = ℓ 1 a 6 , 7 = − 1 a 7 , 6 = 1 / E I 2 a 7 , 7 = ℓ 2 / E I 2 a 7 , 8 = ℓ 2 2 / ( 2 ⋅ E I 2 ) a 8 , 5 = − k C / E I 2 a 8 , 6 = − ( ℓ 2 ⋅ k C ) / E I 2 a 8 , 7 = − ( ℓ 2 2 ⋅ k C ) / ( 2 ⋅ E I 2 ) a 8 , 8 = − ( ℓ 2 3 ⋅ k C − 6 ⋅ E I 2 ) / ( 6 ⋅ E I 2 )
und für die rechte Seite
b 1 = 0 b 2 = 0 b 3 = ( − ( ℓ 1 4 ⋅ q B ) / ( 1 2 0 ⋅ E I 1 ) ) − ( ℓ 1 4 ⋅ q A ) / ( 3 0 ⋅ E I 1 ) b 4 = ( − ( ℓ 1 3 ⋅ q B ) / ( 2 4 ⋅ E I 1 ) ) − ( ℓ 1 3 ⋅ q A ) / ( 8 ⋅ E I 1 ) b 5 = ( − ( ℓ 1 ⋅ q B ) / 2 ) − ( ℓ 1 ⋅ q A ) / 2 b 6 = M B − ( ℓ 1 2 ⋅ q B ) / 6 − ( ℓ 1 2 ⋅ q A ) / 3 b 7 = 0 b 8 = 0 .
📑 toggle code
/* augmented coeff matrix */
ACM: augcoefmatrix ( BCs,ICs ) ;
AA : submatrix ( ACM,9 ) ;
bb : - col ( ACM,9 ) ;
for i: 1 thru 8 do
print ( simplode ( [ "b[" , i, "] = " , string ( bb[i][1] ) ] )) $
for i: 1 thru 8 do
for j: 1 thru 8 do
if not AA[i][j] = 0 then
print ( simplode ( [ "A[" , i, "," , j, "] = " , string ( AA[i][j] ) ] )) $
Solving
Das Lösen des Gleichungssystems liefert
C 1 , 0 = 0 , C 1 , 1 = 2 4 6 . 1 N m 2 , C 1 , 2 = 7 0 3 . 2 N m , C 1 , 3 = − 2 4 0 4 . 3 N , C 2 , 0 = 1 2 7 . 6 N m 3 , C 2 , 1 = 2 2 4 . 7 N m 2 , C 2 , 2 = − 9 7 9 . 8 N m , C 2 , 3 = 2 1 0 1 . 8 N .
📑 toggle code
/* solving */
D : ratsimp ( determinant ( AA )) $
[ P, L, U] : ratsimp ( get_lu_factors ( lu_factor ( AA ))) $
cc : ratsimp ( linsolve_by_lu ( AA,bb ) [1] ) $
sol : makelist ( ICs[i] = cc[i][1],i,1,8 ) $
Post-Processing
Und die Ergebnisse können wir uns anschauen ...
... für w(x):
Biegelinie w(x)
... für Φ(x) :
Kippung w'(x)
... für M(x):
Biegemoment M(x)
... für Q(x):
Querkraft Q(x)
... für die Lager-Reaktionskräfte:
A z = 2 4 0 4 . 3 N , M A = 7 0 3 . 2 N m , B z = 7 4 3 . 9 N , C z = 2 1 0 1 . 8 N , M C = − 1 2 3 . 7 N m
📑 toggle code
/* bearing forces and moments */
reactForces: [A[z]=Q[1] ( 0 ) ,
M[A] = K[A]*Phi[1] ( 0 ) ,
B[z] = k[B]*w[2] ( 0 ) ,
C[z] = k[C]*w[2] ( l[2] ) ,
M[C] = M[2] ( l[2] ) ] ;
expand ( subst ( params,subst ( sol, reactForces ))) ;
/* plot displacements */
fcts: [[ w [1] ( x ) , w [2] ( x-l[1] ) ],
[Phi[1] ( x ) , Phi[2] ( x-l[1] ) ],
[ M [1] ( x ) , M [2] ( x-l[1] ) ],
[ Q [1] ( x ) , Q [2] ( x-l[1] ) ]] ;
facts: [EI[1]/ ( l[1]^4*q[A] ) , EI[1]/ ( l[1]^3*q[A] ) , 1/ ( l[1]^2*q[A] ) , 1/ ( l[1]^1*q[A] ) ] ;
subst ( M[B]/l[1]^2,q[A],facts ) ;
textlabels : [ "w(x)/(M[B]*l^2/EI[1])→" , "w'(x)/(M[B]*l/EI[1])→" , "M(x)/M[B]→" , "Q(x)/(M[B]/l[1]→" ] ;
for i: 1 thru 4 do (
f : expand ( subst ( simple,subst ( xi*l[1],x,facts[i]*[subst ( sol, fcts[i][1] ) ,
subst ( sol, fcts[i][2] ) ] ))) ,
f1 : f[1], f2 : f[2],
toplot : [if xi<=1 then f1 else 0 ,
if xi < 1 then 0 else f2],
plot2d ( toplot,[xi,0,1+subst ( simple,l[2]/l[1] ) ], [legend, "sec. I" , "sec. II" ],
[gnuplot_preamble, "set yrange [] reverse" ] ,
[xlabel, "x/l[1] ->" ],
[ylabel, textlabels[i]] )) $
Links
Literature