Skip to content
Snippets Groups Projects
Commit 14691b02 authored by fsteinme's avatar fsteinme
Browse files

Delete Übung 4 - Lösung.ipynb

parent 68b89c2b
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Übung 4 - Lösung.ipynb deleted 100644 → 0
+ 0
− 319
View file @ 68b89c2b
%% Cell type:markdown id:2457c52a tags:
# Übung 4 - FEM Implementierung Fachwerk
%% Cell type:markdown id:aaf32bb0 tags:
## Imports
%% Cell type:code id:2d991054 tags:
``` python
import numpy as np
import matplotlib.pyplot as plt
import json
import sys, os
from numpy.linalg import solve
```
%% Cell type:markdown id:e5fc02cc tags:
## Problemdefinition: Mesh, Materialien, Randbedingungen und Lasten
Hinweis: Hier muss nichts verändert werden
%% Cell type:code id:875ce8a1 tags:
``` python
# Spezifizieren der Dateinamen -> Geometriedaten werden in andere Dateien ausgegliedert
nodes_file = "nodes.csv"
elems_file = "elems.csv"
disp_file = "disp.csv"
mat_file = "mat.csv"
force_file = "force.csv"
ndf = 2 # Festlegen der Freiheitsgrade (number degrees of freedom)
# Laden der Materialparameter -> par ist ein Python Dictionary
# jf -> json file
with open(os.path.join("data_tower", "par_prob.json"), 'r') as jf:
par = json.load(jf)
# Berechnen der Materialparameter der Stäbe -> [EA, sqrt(2)*EA]
# Notwendig, da die Stäbe unterschiedliche Materialparameter besitzen
# Beispiel mat_par[0] -> Materialparameter des Stabes mit Materialnummer 0
mat_par = []
for key in par.keys():
mat_par.append(par[key]["E"] * par[key]["A"])
# Laden der Materialzuordnung (Index -> Materialparameter ID)
# Beispiel mat[2] -> Materialnummer des Stabs mit ID 2
with open(os.path.join("data_tower", mat_file), "r") as mf:
mat = [int(l.strip()) for l in mf]
# Laden der Knotenkoordinaten
# Aufbau: nodes = [[x0,y0], [x1,y1], ...]
# Zugriff: nodes[Achse][KnotenID]; 0-> x-Achse, 1-> y-Achse
# Beispiel: nodes[0,1] -> x1; nodes[1,3] -> y3
nodes = np.genfromtxt(os.path.join("data_tower", nodes_file), delimiter=',')
# Laden der Elementdaten (Knotenverbindungen)
# Aufbau: elems = [[KnotenID0, KnotenID1], [KnotenID0, KnotenID1], ...]
# Zugriff: elems[ElementID] -> Knoten ID beider Knoten, die das Element verbindet
elems = []
with open(os.path.join("data_tower", elems_file), "r") as ef:
for line in ef:
ids = [int(i) for i in line.strip().split(",")]
elems.append(ids)
# Laden der externen Kräfte
# Aufbau: ID des Knotens, Freiheitsgrad (0:x, 1:y), Betrag der Kraft
force_nodes = []
with open(os.path.join("data_tower", force_file)) as ff:
for line in ff:
sp = line.strip().split(',')
force_nodes.append([int(sp[0]), int(sp[1]), float(sp[2])])
# Laden der Verschiebungsrandbedingungen (nur gesperrte Verschiebung)
# Aufbau: ID des Knotens, Freiheitsgrad (0:x, 1:y)
disp_nodes = []
with open(os.path.join("data_tower", disp_file), 'r') as bcf:
for line in bcf:
sp = line.strip().split(",")
disp_nodes.append([int(s) for s in sp])
```
%% Cell type:markdown id:385bcd6b tags:
## Aufbau der Finite Element Rechnung
%% Cell type:code id:b65cdc82 tags:
``` python
n = len(nodes) # Anzahl der Knoten bestimmen
#TODO: Initialisieren der globalen Stiffness Matrix und des globalen Lastvektors
# in Abhängigkeit der Anzahl der (gesamt) Knoten und der (gesamt) Freiheitsgrade
K_glob = np.zeros((n*ndf, n*ndf))
F_glob = np.zeros(n*ndf)
#Iteration über die Elemente
for elem_id, elem in enumerate(elems):
# Aufbau der Elementsteifigkeitsmatrix
#TODO: Berechnung der Elementlänge
dx = nodes[elem[0]][0] - nodes[elem[1]][0]
dy = nodes[elem[0]][1] - nodes[elem[1]][1]
length = np.sqrt(dx**2 + dy**2)
cos = dx/length
sin = dy/length
#TODO: Aufbau der Elementsteifigkeitsmatrix
k_elem = np.array([[cos**2, cos*sin, -cos**2, -cos*sin], [cos*sin, sin**2, -cos*sin, -sin**2],
[-cos**2, -cos*sin, cos**2, cos*sin], [-cos*sin, -sin**2, cos*sin, sin**2]])
#TODO: Material berücksichtigen
k_elem = k_elem * mat_par[mat[elem_id]] / length
#TODO: Aufbau des Elementresiduums (ohne Volumenkräfte 0)
f_elem = np.zeros(4)
#Element Assembly
for I in range(4):
for J in range(4):
#TODO: Koinzidenzschema implementieren
K_glob[ndf*elem[I//ndf] + I %
ndf, ndf*elem[J//ndf] + J % ndf] += k_elem[I, J]
F_glob[ndf*elem[I//ndf] + I % ndf] += f_elem[I]
```
%% Cell type:markdown id:6f90c5b6 tags:
## Einarbeiten der Randbedingungen
%% Cell type:code id:0c03c294 tags:
``` python
# TODO: Kraftrandbedingungen auf der rechten Seite aufbringen
for fp in force_nodes:
F_glob[ndf*fp[0] + fp[1]] = fp[2]
# IDs der Zeilen/Spalten berechnen, bei denen die Verschiebungsrandbedingungen
# (fester Rand) aufgebracht werden
rem_rows = [ndf*bcp[0] + bcp[1] for bcp in disp_nodes]
# TODO: Löschen der Zeilen/Spalten im LGS - glob_red ist reduziertes LGS
K_glob_red = np.delete(K_glob, rem_rows, 0)
K_glob_red = np.delete(K_glob_red, rem_rows, 1)
F_glob_red = np.delete(F_glob, rem_rows, 0)
#ind ist einfach eine Liste an Indizes von 0 - größe des LGS
ind = np.array([i for i in range(len(nodes)*ndf)])
# hier werden alle gelöschenten Zeilen/Spalten des reduzierten LGS aus ind
# entfernt -> es bleiben so die ursprünglichen Indizes übrig, die noch weiter
# verwendet werden
ind_d = np.delete(ind, rem_rows, 0)
#Initialisierung der kompletten Lösung (inkl. Randbedingungen) als 0-Vektor
result = np.zeros(len(nodes)*ndf)
```
%% Cell type:markdown id:cdcf2fbf tags:
## Lösen des Linearen Gleichungssystems
%% Cell type:code id:f4095332 tags:
``` python
w = solve(K_glob_red, F_glob_red)
#Hier werden die Lösungen der "freien" Knoten in die Gesamtlösung inkl.
#Randbedingungen eingef+gt
result[ind_d] = w
```
%% Cell type:markdown id:22106203 tags:
## post processing (Grafische Ausgabe)
%% Cell type:code id:0b368f1f tags:
``` python
# Berechnung der Reaktionskräfte
F_full = K_glob @ result
# Ausgabe der Reatkionskräfte, Indizes können aus Aufgabenstellung abgelesen werden
print("Reaktion 4: {}, Reaktion 5:{}".format(F_full[4], F_full[5]))
print("Reaktion 6: {}, Reaktion 7:{}".format(F_full[6], F_full[7]))
# Plot des unverformten Fachwerks
fig, ax = plt.subplots(figsize=(10, 10))
# Plot und Beschriftung der Knoten
ax.plot(nodes[:, 0], nodes[:, 1], '.', c='black')
for i in range(nodes.shape[0]):
ax.annotate(str(i), (nodes[i, 0], nodes[i, 1]))
# Plot der Elemente (Stäbe)
for i, e in enumerate(elems):
x = [nodes[e[pi], 0] for pi in [0, 1]]
y = [nodes[e[pi], 1] for pi in [0, 1]]
ax.plot(x, y, '-', c='black')
# Plot des Verformten Fachwerks
scale = 10
# Das Ergebnis (die Knotenverschiebungen) werden auf die Ursprungskoordinaten
# der Knoten addiert und zur besseren Visualisierung mit einem Faktor skaliert
nodes_dp = nodes + result.reshape(result.shape[0]//2,2) * scale
ax.plot(nodes_dp[:, 0], nodes_dp[:, 1], '.', c='red')
for i, e in enumerate(elems):
x = [nodes_dp[e[pi], 0] for pi in [0, 1]]
y = [nodes_dp[e[pi], 1] for pi in [0, 1]]
ax.plot(x, y, '-', c='red')
plt.title("deformed mesh plot - sacled with factor {}".format(scale))
plt.xlabel("x")
plt.ylabel("y")
plt.show()
```
%% Output
Reaktion 4: 0.0, Reaktion 5:-4.163336342344337e-17
Reaktion 6: -5.551115123125783e-17, Reaktion 7:-2.7755575615628914e-17
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment