Gelöste Aufgaben/StaF: Unterschied zwischen den Versionen
Keine Bearbeitungszusammenfassung |
Keine Bearbeitungszusammenfassung |
||
(2 dazwischenliegende Versionen desselben Benutzers werden nicht angezeigt) | |||
Zeile 31: | Zeile 31: | ||
|code= | |code= | ||
<syntaxhighlight lang="lisp" line start=1> | <syntaxhighlight lang="lisp" line start=1> | ||
1 | /*******************************************************/ | ||
/* 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]; | |||
</syntaxhighlight> | </syntaxhighlight> | ||
}} | }} | ||
Zeile 128: | Zeile 173: | ||
|code= | |code= | ||
<syntaxhighlight lang="lisp" line start=1> | <syntaxhighlight lang="lisp" line start=1> | ||
1+1 | /*-----------------------------------------------------*/ | ||
/* 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 ]); | |||
</syntaxhighlight> | </syntaxhighlight> | ||
}} | }} | ||
Zeile 135: | Zeile 196: | ||
{{MyCodeBlock|title=Transformation der Koordinaten in das globale System | {{MyCodeBlock|title=Transformation der Koordinaten in das globale System | ||
|text= | |text= | ||
In den Ausdrücken der virtuellen Formänderungsenergie stehen die Koordinaten des lokalen Koordinatensystems von ''k''. Die müssen wir, wie in [ | In den Ausdrücken der virtuellen Formänderungsenergie stehen die Koordinaten des lokalen Koordinatensystems von ''k''. Die müssen wir, wie in Aufgabe [[Gel%C3%B6ste_Aufgaben/StaB|StaB]] mit der Euler-Drehmatrix ineinander überführen. | ||
Dafür haben wir | Dafür haben wir | ||
Zeile 239: | Zeile 300: | ||
|code= | |code= | ||
<syntaxhighlight lang="lisp" line start=1> | <syntaxhighlight lang="lisp" line start=1> | ||
1+1 | /* 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$ | |||
</syntaxhighlight> | </syntaxhighlight> | ||
}} | }} | ||
Zeile 285: | Zeile 360: | ||
|code= | |code= | ||
<syntaxhighlight lang="lisp" line start=1> | <syntaxhighlight lang="lisp" line start=1> | ||
1 | /*-----------------------------------*/ | ||
/* 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)$ | |||
</syntaxhighlight> | </syntaxhighlight> | ||
}} | }} | ||
Zeile 292: | Zeile 372: | ||
{{MyCodeBlock|title=Solving | {{MyCodeBlock|title=Solving | ||
|text= | |text= | ||
Das Ergebnis ist - für die gleichen Parameter wie in [ | Das Ergebnis ist - für die gleichen Parameter wie in Aufgabe [[Gel%C3%B6ste_Aufgaben/StaB|StaB]]: | ||
::<math> | ::<math> | ||
\begin{array}{l} | \begin{array}{l} | ||
Zeile 309: | Zeile 389: | ||
\end{array} | \end{array} | ||
</math> | </math> | ||
|code= | |||
<syntaxhighlight lang="lisp" line start=1> | |||
/*-----------------------------------*/ | |||
/* 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)))))$ | |||
</syntaxhighlight> | |||
}} | |||
<!--------------------------------------------------------------------------------> | |||
{{MyCodeBlock|title=Post-Processing | |||
|text= | |||
Mit diesen Ergebnissen können wir nun das Stabwerk in seiner verformten Konfiguration zeichnen. | Mit diesen Ergebnissen können wir nun das Stabwerk in seiner verformten Konfiguration zeichnen. | ||
[[Datei:StaF-21.png|350px|right|mini|Stabwerk in verformter Konfiguration.]] | [[Datei:StaF-21.png|350px|right|mini|Stabwerk in verformter Konfiguration.]] | ||
Für die Darstellung müssen wir die Koodinaten der Verschiebung | 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 [[Gelöste_Aufgaben/StaB|StaB]], sind auch die Verläufe der Schnittlasten in den Stäben identisch - wir brauchen also die Ergebnisse nicht noch einmal aufzutragen. | Da die Koordinaten der Verschiebung im globalen Koordinatensystem die gleichen sind wie in Aufgabe [[Gelöste_Aufgaben/StaB|StaB]], sind auch die Verläufe der Schnittlasten in den Stäben identisch - wir brauchen also die Ergebnisse nicht noch einmal aufzutragen. | ||
|code= | |code= | ||
<syntaxhighlight lang="lisp" line start=1> | <syntaxhighlight lang="lisp" line start=1> | ||
1+1 | /*-----------------------------------*/ | ||
/*-----------------------------------------------------*/ | |||
/* 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"]); | |||
</syntaxhighlight> | </syntaxhighlight> | ||
}} | }} | ||
<hr/> | <hr/> | ||
'''Links''' | '''Links''' | ||
* | * * [[Gelöste_Aufgaben/StaB|Gelöste_Aufgaben / StaB]] | ||
'''Literature''' | '''Literature''' | ||
* ... | * ... |
Aktuelle Version vom 27. November 2024, 08:06 Uhr
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.
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
benötigen wir die virtuelle Formänderungsenergie und die virtuelle Arbeit der äußeren Kraft der äußeren Kräfte und Momente.
Mit den Konventionen für die Knoten-Verschiebungen aus Aufgabe StaB ist
- .
Für gilt
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
- und .
Damit haben wir
Für den Stab k definieren wir
- sowie
und finden damit
mit der Element-Steifigkeitsmatrix des Elements k im k-Koordinatensystem
/*-----------------------------------------------------*/
/* 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
mit der Transformationsmatrix
Es ist praktisch, an dieser Stelle die Abkürzung
für die Koordinaten eines Knoten im Referenzsystem einzuführen. Also ist
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
Mit diesen ist
Die resultierenden Element-Steifigkeitsmatrizen sind im folgenden aufgeschreiben:
Element-Steigigkeitsmatrizen mit globalen Koordinaten |
---|
Element #1 |
|
Element #2 |
|
Element #3 |
|
Element #4 |
|
Wir sammeln nun alle Koordianten der Knoten im 0-System in
und schreiben die Gleichgewichtsbedingungen in der Form
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:
/* 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 sowie der passenden Spalten für ein.
Das resultierende Gleichungssystem ist dies:
/*-----------------------------------*/
/* 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:
/*-----------------------------------*/
/* 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.
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
- ...