Gelöste Aufgaben/StaF

Aus numpedia
Zur Navigation springen Zur Suche springen


Aufgabenstellung

Wir untersuchen die Belastung eines ebenen Stabwerks. Die Stäbe haben wie skizziert die Länge ℓ bzw. ℓ/2. Die Struktur wird mit der Kraft F belastet.


Stabwerk mit in den Knoten fest verbundenen Stäben.

Gesucht ist ein Vergleich zwischen der klassischen Stabwerkstheorie und einer Herangehensweise, bei der wir eine feste Verbindung der Stäbe in den Knoten ansetzten. Grundlage des Modells ist die FEM-Lösung der Felddifferentialgleichung im Vergleich zur Lösung in Problemstellung „StaB“.

Wir stellen das Modell des Stabwerks mit dem Prinzip der virtuellen Verrückungen auf und vergleichen, wie sich diese von der Herangehensweise aus „Aufgabe StaB“ mit der analytischen Lösung unterscheidet.

Lösung mit Maxima

Wir nutzen das Computer-Algebra-System Maxima zur Lösung. Das macht hier Sinn, weil wir die Herangehensweise mit der aus Stab vergleichen wollen – für die wir ebenfalls Maxima eingesetzt haben.

Declarations

Wir übernehmen alle Vereinbarungen und Parameter aus der Problemformulierung „Stab“.


/*******************************************************/
/* MAXIMA script                                       */
/* version: wxMaxima 21.05.2                           */
/* author: Andreas Baumgart                            */
/* last updated: 2024-09-20                            */
/* ref: NMM, Labor 2, dimensionsbehaftete              */
/*                           Vorgehensweise            */
/* description: finds the FEM solution for             */
/*              lab problem #2                         */
/*******************************************************/

assume(A>0,H>h,h>0, a>0);

/*-----------------------------------------------------*/
/* Euler-Matrix                                        */
DR(α) := matrix([ cos(α),sin(α), 0],
                [-sin(α),cos(α), 0],
		[   0   ,   0  , 1]);
/* inverse of Euler-Matrix                                        */
DI(α) := transpose(DR(α));
/* compose transforation matrix for rod (two ends ....)*/
DE(α) := matrix([ cos(α),-sin(α), 0,   0   ,  0    , 0],
                [ sin(α), cos(α), 0,   0   ,  0    , 0],
                [   0   ,    0  , 1,   0   ,  0    , 0],
                [   0   ,    0  , 0, cos(α),-sin(α), 0],
                [   0   ,    0  , 0, sin(α), cos(α), 0],
                [   0   ,    0  , 0,   0   ,  0    , 1]);

/* for each rod, define
            * angle α, and
	    * node-IDs
   at its start and end                                */
index : [[ 0  ,[1,3]],  /* rod #1*/
         [α[2],[2,3]],  /* rod #2*/
	 [α[3],[3,4]],  /* rod #3*/
	 [ 0  ,[2,4]]   /* rod #4*/];
/* node - reference locations                          */
nodes: [[0,0],[0,sqrt(3)/2],[1/2,0],[0,sqrt(3)/2]]$

/*-----------------------------------------------------*/
/* system parameters                                   */
params: [ℓ[1]=ℓ[0]/2, ℓ[2]=ℓ[0], ℓ[3]=ℓ[0], ℓ[4]=ℓ[0],
         α[2]=%pi/3,α[3]= -%pi/3,
	 EA = E*A, EI = E*I, I=η*A^2/12];

moreParams: [A = a^2, ℓ[0]= 100*a];




Gleichgewichtsbedingungen

Für die Gleichgewichtsbedingung nach dem Prinzip der virtuellen Verrückungen

δW=!0=δΠδWa

benötigen wir die virtuelle Formänderungsenergie δΠ und die virtuelle Arbeit der äußeren Kraft δWa der äußeren Kräfte und Momente.

Mit den Konventionen für die Knoten-Verschiebungen aus Aufgabe StaB ist

δWa=+δW4,0F.

Für δΠ gilt

δΠ=i=04δΠi

mit den virtuellen Formänderungsarbeiten der vier Stäbe.

Dabei haben wir Anteile der Arbeit aus der Biegung und der Längs-Dehnung des Stabes.

Für den Stab k mit den Knoten I und J haben wir als Koodinaten der Knoten

UI,k,WI,k,ΦI,k und UJ,k,WJ,k,ΦJ,k.

Damit haben wir

δΠk=(δWI,k,δΦI,k,δΦJ,k,δWJ,k)EIi3(126i126i6i4i26i2i2126i126i6i2i26i4i2)(WI,kΦI,kWJ,kΦJ,k)+(δUI,k,δUJ,k)EAi(1111)(UI,kUJ,k)

Für den Stab k definieren wir

Q_k,k=(UI,kWI,kΦI,kUJ,kWJ,kΦJ,k) sowie δQ_k,k=(δUI,kδWI,kδΦI,kδUJ,kδWJ,kδΦJ,k)

und finden damit

δΠk=δQ_k,kTK__k,kQ_k,k

mit der Element-Steifigkeitsmatrix des Elements k im k-Koordinatensystem

K__k,k=EIi3(0000000126i0126i06i4i206i2i20000000126i0126i06i2i206i4i2)+EAi(100100000000000000100100000000000000)

/*-----------------------------------------------------*/
/* nodal coordinates                                   */
nodalCoord : flatten(makelist([U[k,0],W[k,0],Φ[k,0]],k,1,4));
/*****************************************************/
/* Element-Striffness Matrix in local coordinates    */
kI : (EI/ℓ[i]^3)*matrix([ 0,   0  ,     0  , 0,    0  ,     0  ],
                        [ 0,  12  ,  6*ℓ[i], 0,  -12  ,  6*ℓ[i]],			
                        [ 0,6*ℓ[i],4*ℓ[i]^2, 0,-6*ℓ[i],2*ℓ[i]^2],
			[ 0,   0  ,     0  , 0,    0  ,     0  ],
                        [ 0, -12  , -6*ℓ[i], 0,   12  , -6*ℓ[i]],
                        [ 0,6*ℓ[i],2*ℓ[i]^2, 0,-6*ℓ[i],4*ℓ[i]^2])+
     (EA/ ℓ[i] )*matrix([ 1,   0  ,	0  ,-1,    0  ,     0  ],
			[ 0,   0  ,  	0  , 0,    0  ,     0  ],
			[ 0,   0  ,	0  , 0,    0  ,     0  ],
			[-1,   0  ,	0  , 1,    0  ,     0  ],
			[ 0,   0  ,	0  , 0,    0  ,     0  ],
			[ 0,   0  ,	0  , 0,    0  ,     0  ]);




Transformation der Koordinaten in das globale System

In den Ausdrücken der virtuellen Formänderungsenergie stehen die Koordinaten des lokalen Koordinatensystems von k. Die müssen wir, wie in Aufgabe StaB mit der Euler-Drehmatrix ineinander überführen.

Dafür haben wir

(UI,0WI,0ΦI,0)=D__R(αk)(UI,kWI,kΦI,k)

mit der Transformationsmatrix

D__R(αk)=(cos(αk)sin(αk)0sin(αk)cos(αk)0001).

Es ist praktisch, an dieser Stelle die Abkürzung

q_I,0=(UI,0WI,0ΦI,0)

für die Koordinaten eines Knoten im Referenzsystem einzuführen. Also ist

(UI,kWI,kΦI,k)=D__RT(αk)q_I,0.

Damit wir für die Elementsteifigkeitsmatrix - mit beiden Anfangs- und Endknoten des Elements - vom "0"-System ins "k"-System transformieren, brauchen wir die neue Transformations-Matrix

T__(αk):=(D__RT(αk)0__0__D__RT(αk))

Mit diesen ist

δΠ=k=14δQ_k,0TT__T(αk)K__k,kT__(αk):=K__k,0Q_k,0

Die resultierenden Element-Steifigkeitsmatrizen sind im folgenden aufgeschreiben:

Element-Steigigkeitsmatrizen mit globalen Koordinaten
Element #1

K__1,0=(8A2Eη0302A2Eη02(8A2Eη03)02A2Eη0202AE000(2AE0)02A2Eη0202A2Eη30(2A2Eη02)0A2Eη30(8A2Eη03)0(2A2Eη02)8A2Eη030(2A2Eη02)0(2AE0)002AE002A2Eη020A2Eη30(2A2Eη02)02A2Eη30)

Element #2

K__2,0=(A2Eη+302AE4033A2Eη302AE403A2Eη402(A2Eη+302AE403)(3A2Eη302AE403)A2Eη4023A2Eη302AE4033A2Eη+02AE4033A2Eη402(3A2Eη302AE403)(3A2Eη+02AE403)3A2Eη402A2Eη4023A2Eη402A2Eη30(A2Eη402)(3A2Eη402)A2Eη60(A2Eη+302AE403)(3A2Eη302AE403)(A2Eη402)A2Eη+302AE4033A2Eη302AE403(A2Eη402)(3A2Eη302AE403)(3A2Eη+02AE403)(3A2Eη402)3A2Eη302AE4033A2Eη+02AE403(3A2Eη402)A2Eη4023A2Eη402A2Eη60(A2Eη402)(3A2Eη402)A2Eη30)

Element #3

K__3,0=(A2Eη+302AE403(3A2Eη302AE403)A2Eη402(A2Eη+302AE403)3A2Eη302AE403A2Eη402(3A2Eη302AE403)3A2Eη+02AE403(3A2Eη402)3A2Eη302AE403(3A2Eη+02AE403)(3A2Eη402)A2Eη402(3A2Eη402)A2Eη30(A2Eη402)3A2Eη402A2Eη60(A2Eη+302AE403)3A2Eη302AE403(A2Eη402)A2Eη+302AE403(3A2Eη302AE403)(A2Eη402)3A2Eη302AE403(3A2Eη+02AE403)3A2Eη402(3A2Eη302AE403)3A2Eη+02AE4033A2Eη402A2Eη402(3A2Eη402)A2Eη60(A2Eη402)3A2Eη402A2Eη30)

Element #4

K__4,0=(A2Eη030A2Eη202(A2Eη03)0A2Eη2020AE000(AE0)0A2Eη2020A2Eη30(A2Eη202)0A2Eη60(A2Eη03)0(A2Eη202)A2Eη030(A2Eη202)0(AE0)00AE00A2Eη2020A2Eη60(A2Eη202)0A2Eη30)

Wir sammeln nun alle Koordianten der Knoten im 0-System in

Q_0T=(q_1,0T,q_2,0T,q_3,0T,q_4,0T)

und schreiben die Gleichgewichtsbedingungen in der Form

δW=δQ_0T(K__0Q_0P_)=!0

an. Dabei kommen die Beiträge zur Gesamt-Steifigkeitsmatrix (hier noch in der Fassung ohne Berücksichtigung der Randbedingungen) aus den vier Beiträgen der virtuellen Formänderungsenergie - die wir hier farblich gekennzeichnet haben:

Einsortieren der Anteile der virtuellen Formänderungsenergie nach den Knoten.
Einsortieren der Anteile der virtuellen Formänderungsenergie nach den Knoten.

/* apply transformation to global coordinates         */
elemStiffMat: ratsimp(subst(params,makelist(transpose(DE(index[k][1])).subst(i=k,kI).DE(index[k][1]),k,1,4)))$
for i:1 thru 4 do
    print(k[i],"=",elemStiffMat[i])$
    
K: zeromatrix(length(nodalCoord),length(nodalCoord))$
for elem: 1 thru 4 do
   (pivot:append(3*(index[elem][2][1]-1)*[1,1,1],
                 3*(index[elem][2][2]-1)*[1,1,1])+[1,2,3,1,2,3],
    for row: 1 thru 6 do
        for col: 1 thru 6 do
       	    K[pivot[row],pivot[col]] : K[pivot[row],pivot[col]]+elemStiffMat[elem][row,col])$
        
B: zeromatrix(length(nodalCoord),1)$
B[11,1]:F$




Einarbeitung der Randbedingungen

Die Randbedingungen arbeiten wir hier durch das Streichen der passenden Zeilen für δU1,0,δW1,0,δU2,0,δW2,0 sowie der passenden Spalten für U1,0,W1,0,U2,0,W2,0 ein.

Das resultierende Gleichungssystem ist dies:

(2A2Eη3000(2A2Eη02)A2Eη3000002A2Eη30(3A2Eη402)(A2Eη402)A2Eη600(A2Eη202)A2Eη600(3A2Eη402)3A2Eη+02AE203+2AE00(3A2Eη202)(3A2Eη+02AE403)3A2Eη302AE403(3A2Eη402)(2A2Eη02)(A2Eη402)0A2Eη+302AE203+8A2Eη03(2A2Eη02)3A2Eη302AE403(A2Eη+302AE403)A2Eη402A2Eη30A2Eη60(3A2Eη202)(2A2Eη02)4A2Eη303A2Eη402(A2Eη402)A2Eη6000(3A2Eη+02AE403)3A2Eη302AE4033A2Eη4023A2Eη+02AE403+AE0(3A2Eη302AE403)3A2Eη4020(A2Eη202)3A2Eη302AE403(A2Eη+302AE403)(A2Eη402)(3A2Eη302AE403)A2Eη+302AE403+A2Eη03(3A2Eη402)0A2Eη60(3A2Eη402)A2Eη402A2Eη603A2Eη402(3A2Eη402)2A2Eη30)(Φ1,0Φ2,0W3,0U3,0Φ3,0W4,0U4,0Φ4,0)=(00000F00)

/*-----------------------------------*/
/* boundary conditions */
K: submatrix(1,2,4,5,K,1,2,4,5)$
X: submatrix(1,2,4,5,transpose(nodalCoord))$
B: submatrix(1,2,4,5,B)$
print(K,"*",X,"=",B)$




Solving

Das Ergebnis ist - für die gleichen Parameter wie in Aufgabe StaB:

U1,0=0W1,0=0Φ1,0=3.97FEa2U2,0=0W2,0=0Φ2,0=2.45FEa2U3,0=57.7FEaW3,0=1.67102FEaΦ3,0=2.06FEa2U4,0=57.7FEaW4,0=3.67102FEaΦ4,0=3.12FEa2

/*-----------------------------------*/
/* solving                           */
sol: ratsimp(linsolve_by_lu(K,B)[1]);
sol: append(makelist([U[1,0],W[1,0],U[2,0],W[2,0]][i]=0,i,1,4),makelist(X[i][1]=sol[i][1],i,1,length(X)));
sol: makelist(nodalCoord[i]=subst(sol,nodalCoord[i]),i,1,length(nodalCoord));

fpprintprec: 3;
print(transpose(float(subst(η=1,subst(moreParams,sol)))))$




Post-Processing

Mit diesen Ergebnissen können wir nun das Stabwerk in seiner verformten Konfiguration zeichnen.

Stabwerk in verformter Konfiguration.

Für die Darstellung müssen wir die Koodinaten der Verschiebung in die lokalen Koordinatensysteme zurücktransformieren - und das geht wieder über die Euler-Drehmatrizen (s.o.).

Da die Koordinaten der Verschiebung im globalen Koordinatensystem die gleichen sind wie in Aufgabe StaB, sind auch die Verläufe der Schnittlasten in den Stäben identisch - wir brauchen also die Ergebnisse nicht noch einmal aufzutragen.


/*-----------------------------------*/
/*-----------------------------------------------------*/
/* post-processing: show deformed structure            */

/* employ these Trial-Functions - shared from FEM      */
φ : [      (ξ-1)^2*(2*ξ+1),
      ℓ[i]* ξ     *(  ξ-1)^2,
     -      ξ^2   *(2*ξ-3), 
      ℓ[i]* ξ^2   *(  ξ-1)];

/* construct u[i](x[i]),w[i](x[i])                     */
/* step 1: re-construct local nodal-coordinates (DI(α))*/
/* step 2: use trial function to derive u[i],w[i]      */
/* step 3: transform into global coordinates (DR(α))   */
sol: ratsimp(subst([η=1],subst(moreParams,sol)))$

fct:[]$
for e:1 thru 4 do
   (I: index[e][2][1],
    J: index[e][2][2],
    locals: [subst(params,matrix([U[I,e]],[W[I,e]],[Φ[I,e]])=DI(index[e][1]).matrix([U[I,0]],[W[I,0]],[Φ[I,0]])),
             subst(params,matrix([U[J,e]],[W[J,e]],[Φ[J,e]])=DI(index[e][1]).matrix([U[J,0]],[W[J,0]],[Φ[J,0]]))],
    trafo: flatten(makelist(makelist(lhs(locals[j])[i][1]=rhs(locals[j])[i][1],i,1,3),j,1,2)),
    /* ξ for x-Axis .....                              */
    localFct:matrix([ξ*ℓ[i]/ℓ[0] + U[I,e]*(1-ξ)+ U[J,e]   ],
                    [              W[I,e]*φ[1] + Φ[I,e]*φ[2]+W[J,e]*φ[3]+Φ[J,e]*φ[4]],[1]),
    localFct: subst(trafo,localFct),	    
    localFct: subst(params,subst([i=e],DR(index[e][1]).localFct)),
    fct: append(fct,[flatten(makelist(localFct[i],i,1,2))]))$

/* choose scale of deflection                          */
scale: [E=1000*F/a];
fct: float(expand(subst([scale],subst(moreParams,subst(params,subst(sol,fct))))));

/* and plot .....                                      */
toPlot: makelist([parametric, nodes[i][1] + fct[i][1],
                              nodes[i][2] + fct[i][2], [ξ,0,1]],i,1,4);
plot2d(toPlot, [legend,"Stab 1","Stab 2","Stab 3","Stab 4"]);





Links


Literature

  • ...