Knickung

Copyright: Rolf Wirz

1. Knickung eines nicht eingespannten druckbelasteten Stabes mit einer hypothetischen Biegelinie

a    Definitionen, Vorbereitungen:

Balken- oder Trägerlänge xL = 5 m; 
Balken, Rohrquerschnitt: Aussenradius R = 0.01 m; Innenradius r = 0.005 m; 
axiales Flächenträgheitsmoment Iy  = Pi/4 (R^4 - r^4) (Masse übernommen); 
Elastizitätsmodul eE = 210000 *1/(1/1000^2) N/ (m^2); 
Kraft F = 10^6 N;

Remove["Global`*"]

a.    Approximierte die Differentialgleichung, exakte Lösung der Approximation der Differentialgleichung, Shooting-Methode

M[y_]:= F y;
y''[x] == -M[y[x]]/(eE Iy)

solv = DSolve[{y''[x] == -M[y[x]]/(eE Iy)},y,x]/.E^(0. x)->1

solv = DSolve[{y''[x] == -M[y[x]]/(eE Iy),y[0]==0},y,x]/.E^(0. x)->1

xL = 5; R = 0.1; r = 0.01; Iy = Pi/4 (R^4 - r^4); 
eE = 210000 *1/(1/1000^2); F = 5 10^6;  eE

M[y_]:= F y;
y''[x] == -M[y[x]]/(eE Iy)

solv = DSolve[{y''[x] == -M[y[x]]/(eE Iy)},y,x]/.E^(0. x)->1

solv = DSolve[{y''[x] == -M[y[x]]/(eE Iy),y[0]==0},y,x]/.E^(0. x)->1

solv = DSolve[{y''[x] == -M[y[x]]/(eE Iy),y[0]==0,y[xL]==0}, y,x] /.E^(0. x)->1

DSolve[{y''[x] == -k y[x]},y,x]

Keine Lösung für dieses xL ausser die Null-Lösung, welche nicht interessiert.
Eine Lösung kann man nur finden, wenn xL so geartet ist, dass (Pi/xL)^2= k=F/(eE Iy)

DSolve[{y''[x] == -(Pi/xL0)^2 y[x],y[0]==0},y,x]

sol = DSolve[{y''[x] == -(Pi/xL0)^2 y[x],y[xL0]==0},y,x] // Flatten

y''[x]== -M[y[x]]/(eE Iy) /.sol

Solve[(((π^2) )/(xL0^2)) == 
        F1/(eE1 Iy1), {xL0}]

C[2] lässt sich aus dieser Gleichung nicht bestimmen!

Solve[(((π^2) )/(xL0^2)) == 
        F1/(eE1 Iy1), {F1}]

Wenn F1(xL0) den oben errechneten Wert hat, so hat die Differentialgleichung eine Lösung, welche einer sinusartigen Biegelinie entspricht. Weil aber die Auslenkung y[x] infolge fehlendem C[2] nicht bekannt ist, kann auch die maximale Biegespannung wegen dem folglich fehlenden Biegemoment auf diesem Weg nicht ermittelt werden. Es ist daher vorläufig nicht entscheidbar, ob der Stab bricht.
Man beachte, dass die Bestimmung von C[2] aus diesem Ranswertproblem aus logischen Gründen unmöglich ist. Das Problem lässt sich abstrakt wie folgt formulieren:
y''[x] = -k*y[x],  y[x] = y[xL0] = 0. Wenn man alle Gleichungen hier mit c multipliziert, so ergibt sich:  c*y''[x] = -c*k*y[x],  c*y[x] = c*y[xL0] = 0, also mit z[x] = c*y[x]:
z''[x] = -k*z[x],  z[x] = z[xL0] = 0. Das Randwertproblem ist daher mit y[x] auch für c*y[x] erfüllt, unabhängig von c.
Glücklicherweise hilft hier noch ein anderes Argument: In der eben gewonnenen Formel  F  = ((eE * Iy * (π^2))/(xL^2)) ist F diejenige Kraft, welche bei der Trägerlänge xL für die vorhandene Biegelinie sorgt. In einem solchen durch F verursachten Spannungszustand herrschen in den Punkten des Trägers die Gleichgewichtsbedingungen: Die Summe der Kräfte und die Summe der Momente sind gleich null. Ist die bei der Stablänge xL gegebene Kraft nun grösser als das F in der letzten Formel, so muss xL bei gleich bleigendem eE und Iy kleiner werden. Wenn das nicht möglich ist, so fällt der Gleichgewichtszustand zusammen. Der Stab muss brechen!
F  = ((eE * Iy * (π^2))/(xL^2)) heisst übrigens "Euler-Knickkraft".
In der nun unlösbaren Situation mit dem unbekannten C[2]  versuchen wir es noch mit der exakten Differentialgleichung und einer numerischen Lösung:

b. Exakte Differentialgleichung, numerische Lösung (Shooting-Methode). 
    (Exakt zu schwierig allgemein zu lösen)

(* Remove["Global`*"] *)
Remove[x,y]

M[y_]:= F y;

Remove[x, y, solution1];
ys0 = 
        0.1;
solution1=NDSolve[{Derivative[2][y][x]+F y[
                          x]/(eE Iy) (1 + (y'[x])^2)^(3/
                            2)==0,y[0] == 0,y'[0]==ys0},y,{x, 0, 1.2 xL}];
k1 = 
        Plot[y[x] /. solution1, {x, 0, 
            1.2 xL}];
ys0 = 
        0.2;
solution2=NDSolve[{Derivative[2][y][x]+F y[
                          x]/(eE Iy) (1 + (y'[x])^2)^(3/
                            2)==0,y[0] == 0,y'[0]==ys0},y,{x, 0, 1.2 xL}];
k2 = 
        Plot[y[x] /. solution2, {x, 0, 
            1.2 xL}];
ys0 = 
        0.3;
solution3=NDSolve[{Derivative[2][y][x]+F y[
                          x]/(eE Iy) (1 + (y'[x])^2)^(3/
                            2)==0,y[0] == 0,y'[0]==ys0},y,{x, 0, 1.2 xL}];
k3 = 
        Plot[y[x] /. solution3, {x, 0, 
            1.2 xL}];
ys0 = 
        0.4;
solution4=NDSolve[{Derivative[2][y][x]+F y[
                          x]/(eE Iy) (1 + (y'[x])^2)^(3/
                            2)==0,y[0] == 0,y'[0]==ys0},y,{x, 0, 1.2 xL}];
k4 = 
        Plot[y[x] /. solution4, {x, 0, 
            1.2 xL}];
ys0 = 
        1;
solution5=NDSolve[{Derivative[2][y][x]+F y[
                          x]/(eE Iy) (1 + (y'[x])^2)^(3/
                            2)==0,y[0] == 0,y'[0]==ys0},y,{x, 0, 1.2 xL}];
k5 = 
        Plot[y[x] /. solution5, {x, 0, 
            1.2 xL}];
ys0 = 
        1.1;
solution6=NDSolve[{Derivative[2][y][x]+F y[
                          x]/(eE Iy) (1 + (y'[x])^2)^(3/
                            2)==0,y[0] == 0,y'[0]==ys0},y,{x, 0, 1.2 xL}];
k6 = 
        Plot[y[x] /. solution6, {x, 0, 1.2 xL}];

Show[k1,k2,k3,k4,k5,k6];

ys0 = 1.05742;
solution7=NDSolve[{Derivative[2][y][x]+F y[
                          x]/(eE Iy) (1 + (y'[x])^2)^(3/
                            2)==0,y[0] == 0,y'[0]==ys0},y,{x, 0, 1.2 xL}];
k7 = 
        Plot[y[x] /. solution7, {x, 0, 
            1.2 xL}];
yE = (y[x] /. 
            solution7) /. x -> xL

ys0 = 1.05741;
solution7=NDSolve[{Derivative[2][y][x]+F y[
                          x]/(eE Iy) (1 + (y'[x])^2)^(3/
                            2)==0,y[0] == 0,y'[0]==ys0},y,{x, 0, 1.2 xL}];
k7 = 
        Plot[y[x] /. solution7, {x, 0, 
            1.2 xL}];
yE = (y[x] /. 
            solution7) /. x -> xL

FAZIT:  Das Randwertproblem ist äquivalent zu einem Anfangswertproblem mit  y[0] = 0 und  y'[0] irgendwo im Intervall 1.05741 < y'[0] <  1.05742

2. Programm für das automatische Ausführen der Shooting-Methode

Remove["Global`*"]

(* Initialisierung *)
xL = 5; R = 0.1; r = 0.01; Iy = Pi/4 (R^4 - r^4); 
eE = 210000 *1/(1/1000^2); F = 5 10^6;  
M[y_]:= F y;
ys0=1; i=1;
solution=
NDSolve[{y''[x] + F y[x]/(eE Iy) (1+(y'[x])^2)^(3/2)==0,
y[0]==0, y'[0] == ys0 },y,{x,0,1.2 xL}];
Endwert=(y[x]/.solution)/.x->xL;
If[Endwert[[1]] > 0, steuerung = 1, If[Endwert[[1]] < 0, steuerung = -1, steuerung = 0]];
 h=steuerung;
pl[i]=Plot[y[x]/.solution,{x,0,1.2 xL},DisplayFunction->Identity];
Print[{i,h,steuerung,Endwert " = y[xL]",ys0}];

(* Loop *)
While[(h-steuerung) == 0 ,
Print["Vor ",{i,h,steuerung,Endwert " = y[xL]",ys0//N}];
ys0=1+i/1000; 
i++; 
h=steuerung;
solution=
NDSolve[{y''[x] + F y[x]/(eE Iy) (1+(y'[x])^2)^(3/2)==0,
y[0]==0, y'[0] == ys0 },y,{x,0,1.2 xL}];
Endwert=(y[x]/.solution)/.x->xL;
If[Endwert[[1]] > 0, steuerung = 1, If[Endwert[[1]] < 0, steuerung = -1, steuerung = 0]];
pl[i]=Plot[y[x]/.solution,{x,0,1.2 xL},DisplayFunction->Identity];
Print["Nach ",{i,h,steuerung,Endwert " = y[xL]",ys0//N},If[(h-steuerung) == 0,""," Umschlagpunkt"]];
];
Show[Table[pl[n],{n,1,i}],DisplayFunction->$DisplayFunction];
Plot[y[x]/.solution,{x,0,1.2 xL}];

3. An einem Ende fest eingespannter, am anderen Ende frei geführter Stab

Damit an der Einspannstelle der Stab senkrecht gehalten wird, muss in diesem Fall auch eine horizontale Querkraft FQ vorhanden sein. Aus Gleichgewichtsgründen finden wir diese Situation an anderen Stabende gespiegelt wieder. Daher entsteht jetzt im Abstand x vom Stabende noch ein Moment mit dem Betrage FQ * x, welches entgegengesetzt zu F * y gereichtet ist. Damit erhalten wir bei den obigen Bedingungen die folgenden Gleichungen:

Remove["Global`*"]

M[y_,x_]:= F y -FQ x;
y''[x] == -M[y[x]]/(eE Iy)

solv = DSolve[{y''[x] == -M[y[x],x]/(eE Iy)},y,x]/.E^(0. x)->1

Wegen der Randbedingung y[0] = 0 muss hier C[1] = 0 sein. Aus den verbleibenden Teilen erhält man dann wieder Bedingungen für die Abhängigkeiten zwischen FQ, F und C[2] bei gegebenen eE und Iy. Ebenso erhält man weitere Bedingungen mittels des bekannten Randwerts der Ableitung von y[x].

{
FQ xL /F == - C[2] Sin[Sqrt[F] xL /(Sqrt[eE Iy])], 
FQ /F == - C[2]Sqrt[F]/Sqrt[eE Iy] Cos[Sqrt[F] xL/(Sqrt[eE Iy])] == 0
}

xL == Sqrt[eE Iy]/Sqrt[F] Tan[Sqrt[F] xL /(Sqrt[eE Iy])]

xL = 5; R = 0.1; r = 0.01; Iy = Pi/4 (R^4 - r^4); 
eE = 210000 *1/(1/1000^2);

xL Sqrt[F]== Sqrt[eE Iy] Tan[Sqrt[F] xL /(Sqrt[eE Iy])]

FindRoot[xL Sqrt[F]== Sqrt[eE Iy] Tan[Sqrt[F] xL /(Sqrt[eE Iy])],{F,1}]

FindRoot[xL Sqrt[F]== Sqrt[eE Iy] Tan[Sqrt[F] xL /(Sqrt[eE Iy])],{F,0.2}]

FindRoot[xL Sqrt[F]== Sqrt[eE Iy] Tan[Sqrt[F] xL /(Sqrt[eE Iy])],{F,0.5}]

FindRoot[xL Sqrt[F]== Sqrt[eE Iy] Tan[Sqrt[F] xL /(Sqrt[eE Iy])],{F,10^6}]

FindRoot[xL Sqrt[F]== Sqrt[eE Iy] Tan[Sqrt[F] xL /(Sqrt[eE Iy])],{F,5*10^6}]

FindRoot[xL Sqrt[F]== Sqrt[eE Iy] Tan[Sqrt[F] xL /(Sqrt[eE Iy])],{F,10^7}]

Plot[{xL Sqrt[F], 
          Sqrt[eE Iy] Tan[Sqrt[F] xL/(Sqrt[eE Iy])]}, {F,  0,  
          3*10^7}];

F lässt sich bei günstigen Koeffizientenverhältnissen aus dieser Gleichung numerisch bestimmen. Bei klein gewählten F entstehen hier Werte, welche wegen der numerischen Ungenauigkeit keinen praktischen Wert besitzen. F ist ungefähr 1.33 * 10^7. 
Hier wollen wir jedoch das Problem einmal numerisch weiter verfolgen:

xL = 5; R = 0.1; r = 0.01; Iy = Pi/4 (R^4 - r^4); 
eE = 210000 *1/(1/1000^2); F = 5 10^6;

solv = DSolve[{y''[x] == -M[y[x],x]/(eE Iy)},y,x]/.E^(0. x)->1

solv = DSolve[{y''[x] == -M[y[x],x]/(eE Iy),y[0]==0,y[xL]==0}, y,x] /.E^(0. x)->1

solv = DSolve[{y''[x] == -M[y[x],x]/(eE Iy),y[0]==0,y'[0]==0,y[xL]==0}, y,x] /.E^(0. x)->1

Die Gleichung 2. Ordnung hat bei den drei Rand- und Anfangsbedingungen  
y[0] = 0,y'[0] = 0,y[xL] = 0 keine allgemeine Lösung. 
Wenn man nun die gegebene Differentialgleichung nochmals differenziert, so erhöht sich die Ordnung. Dafür geht andererseits Information verloren. (Die Ableitung einer Konstanten ist 0.)

Remove["Global`*"]

M[y_,x_]:= F y -FQ x;

(* Einmal differenzieren: *)

y'''[x] == D[-M[y[x],x]/(eE Iy),x]

solv = DSolve[{y'''[x] == D[-M[y[x],x]/(eE Iy),x]},y,x]/.E^(0. x)->1

solv = (DSolve[{y'''[x] == D[-M[y[x],x]/(eE Iy),x],y[0]==0,y'[0]==0,y[xL]==0}, y,x] /.E^(0. x)->1 ) // Simplify

Nochmals differenzieren:

solv = DSolve[{y''''[x] == D[-M[y[x],x]/(eE Iy),{x,2}]},y,x]/.E^(0. x)->1

solv = (DSolve[{y''''[x] == D[-M[y[x],x]/(eE Iy),{x,2}],y[0]==0,y'[0]==0,y[xL]==0}, y,x] /.E^(0. x)->1 )

(y[x]/.solv) /.x->xxxxxx

Die Lösung ist bei fixem A, B, C, D von der Form:

y[x] == C[4] (A + B x + C Sin[x Sqrt[F/(eE Iy)]] + D Cos[x Sqrt[F/(eE Iy)]])

Wegen y[0] == 0 ist D = -A.

xL = 5; R = 0.1; r = 0.01; Iy = Pi/4 (R^4 - r^4); 
eE = 210000 *1/(1/1000^2); F = 5 10^6;

solv = (DSolve[{y''''[x] == D[-M[y[x],x]/(eE Iy),{x,2}],y[0]==0,y'[0]==0,y[xL]==0}, y,x] /.E^(0. x)->1 )

((y[x]/.solv) /.x->0) // Chop

(y[x]/.solv) /.x->xL

Problem: Erfüllt diese Lösung auch die Gleichungen 3. und 2. Ordnung?

solv = (DSolve[{y''''[x] == D[-M[y[x],x]/(eE Iy),{x,2}],y[0]==0,y'[0]==0,y[xL]==0}, y,x] /.E^(0. x)->1 )

(* 3. Ordnung: *)
0 == ((y'''[x] - D[-M[y[x],x]/(eE Iy),x])/.solv)[[1]] // Simplify  // Chop

Hier ist also C[4] durch FQ bestimmt.

(* 2. Ordnung: *)
0 == ((y''[x] + M[y[x],x]/(eE Iy))/.solv)[[1]] // Simplify  // Chop //ExpandAll

Hier lässt sich x nicht ausklammern und wegkürzen. (z.B. für x = 0 wird C[4] = 0.)  Man sieht sofort, dass die letzte Gleichung nur für C[4] = FQ = 0 bestehen kann. Nun lassen wir F frei:

Remove["Global`*"]

xL = 5; R = 0.1; r = 0.01; Iy = Pi/4 (R^4 - r^4); 
eE = 210000 *1/(1/1000^2);

M[y_,x_]:= F y -FQ x;

solv = (DSolve[{y''''[x] == D[-M[y[x],x]/(eE Iy),{x,2}],y[0]==0,y'[0]==0,y[xL]==0}, y,x] /.E^(0. x)->1 )

(* 2. Ordnung: *)
((y''[x] + M[y[x],x]/(eE Iy)/.solv)[[1]] // Simplify // Together  // Chop //ExpandAll // Together) == 0

Hier hat man einen Ausdruck der Form
FQ x (1 - Cos[A Sqrt[F]]) - F x C[4](1 - Cos[A Sqrt[F]]) == B Sqrt[F] C[4] Sin[A Sqrt[F]] - C F C[4] für alle x.
Daraus folgt FQ = F C[4] und B Sin[A Sqrt[F]] = C Sqrt[F]. In Anbetracht der Konstanten sollte man F approximativ berechnen können, C[4] oder FQ bleiben hier jedoch wiederum frei. Das bedeutet, wie oben, dass man die exakte Differentialgleichung behandeln muss.