'Main'-funksjon
I main har vi laget en hovedramme for hvordan programmet kan bygges opp. Her er også forslag til flere av funksjonskallene inkludert slik at det blir lettere å komme igang med programmeringen. Legg merke til at funksjonene lesinput og , lengder allerede og all visualisering allerede er implementert. Flere detaljer om disse funksjonene kommer i neste avsnitt.
| Code Block | ||||||
|---|---|---|---|---|---|---|
| ||||||
importfrom numpystructure_visualization asimport np* # -----Rammeanalyse----- def main(): # Rammeanalyse # -----LeserInitialiserer input-datafigurer----- npunktfig_init, punktax_init, nelemfig_def, elem, nlast, last = lesinputax_def = setup_plots() # -----Regner ut lengder til elementeneTil visualiseringen, velg første indeks brukt i nummerering av noder og element------ # Endre gjerne selv elementlengder = lengder(punkt, elem, nelem) first_index = 0 # -----FastinnspenningsmomenteneLeser input-data----- # Lag funksjon selv fim = moment(npunkt, punktnpunkt, punkt, nelem, elem, nlast, last, = elementlengderlesinput() # -----SetterPlott opp lastvektorinitalramme----- # Lag funksjon selvplot_structure(ax_init, punkt, elem, 1, first_index) b# -----Regner ut lengder til elementene------ elementlengder = lastvektorlengder(fimpunkt, npunkt, punktelem, nelem, elem) # -----Fastinnspenningsmomentene-Setter opp systemstivhetsmatrisen----- # Lag funksjonfunksjonen selv Kfim = stivhet(moment(npunkt, punkt, nelem, elem, elementlengdernlast, npunktlast, elementlengder) # ------Innfører randbetingelser-Setter opp lastvektor----- # Lag funksjonfunksjonen selv Kn, Bnb = bclastvektor(fim, npunkt, punkt, Knelem, belem) # -----Løser ligningssystemet--Setter opp systemstivhetsmatrisen----- # Lag funksjonfunksjonen selv rotK = ... # Hint, se side for løsing av lineære systemer i Python stivhet(nelem, elem, elementlengder, npunkt) # ------Finner endemoment for hvert elementInnfører randbetingelser------ # Lag funksjonfunksjonen selv Kn, endemomentBn = endeMbc(npunkt, punkt, nelemK, elem, elementlengder, rot, fimb) # -----Skriver ut hva rotasjonen ble i de forskjellige nodene----- print("Rotasjoner i de ulike punktene:") print(rot) #Løser ligningssystemet------ # Lag funksjonen selv rot = ... # Hint, se side for løsing av lineære systemer i Python #------Skriver ut hva momentene bleFinner endemoment for dehvert forskjellige elementeneelement----- print("Elementvis endemoment:") # Lag funksjonen selv endemoment print(endemoment) |
| Info |
|---|
Dersom du ønsker tips til å forbedre kodingen din, besøk siden Tips and tricks for coding. Andre nyttige sider for dette prosjektet kan være LaTeX, Inkscape TexText og Python. |
Utdelte funksjoner
| Code Block | ||||||
|---|---|---|---|---|---|---|
| ||||||
import numpy as np def lengder(knutepunkt, element, nelem):= endeM(npunkt, punkt, nelem, elem, elementlengder, rot, fim) elementlengder = np.zeros((nelem, 1)) # Beregner elementlengder med Pythagoras' laeresetning#-----Skriver ut hva rotasjonen ble i de forskjellige nodene----- forprint("Rotasjoner i inde range (0, nelem):ulike punktene:") print(rot) # OBS! Grunnet indekseringsyntaks i Python-arrays vil ikke denne funksjonen fungere naar vi bare har ett element.-----Skriver ut hva momentene ble for de forskjellige elementene----- print("Elementvis endemoment:") dx = knutepunkt[element[i, 0], 0] - knutepunkt[element[i, 1], 0] print(endemoment) #-----Plott deformert ramme----- skalering = 100; # Du kan dyendre = knutepunkt[element[i, 0], 1] - knutepunkt[element[i, 1], 1]denne konstanten for å skalere de synlige deformasjonene til rammen plot_structure_def(ax_def, punkt, elem, elementlengder[i] = np.sqrt(dx*dx + dy*dy) return elementlengder1, first_index, skalering*rot) plt.show() |
| Info |
|---|
Dersom du ønsker tips til å forbedre kodingen din, besøk siden Tips and tricks for coding. Andre nyttige sider for dette prosjektet kan være LaTeX, Inkscape TexText og Python. |
Utdelte funksjoner
Lengder
| Code Block | ||||||
|---|---|---|---|---|---|---|
| ||||||
import numpy as np def lesinput(def lengder(knutepunkt, element, nelem): # Åpner inputfilen fid = open("input.txt", "r") elementlengder = np.zeros((nelem, 1)) # LeserBeregner totaltelementlengder antallmed punkt npunkt = int(fid.readline())Pythagoras' læresetning for i in range (0, nelem): # 'fid.readline()' leser en linje, 'int(...)' gjør at linjen tolkes som et heltall # LESER INN XY-KOORDINATER TIL KNUTEPUNKTENE # Nodenummer tilsvarer radnummer i "Node-variabel"OBS! Grunnet indekseringsyntaks i Python-arrays vil ikke denne funksjonen fungere naar vi bare har ett element. dx = knutepunkt[element[i, 0], 0] - knutepunkt[element[i, 1], 0] # x-koordinat lagres i kolonne 1, y-koordinat i kolonne 2 dy = knutepunkt[element[i, 0], 1] - knutepunkt[element[i, 1], 1] # Grensebetingelse lagres elementlengder[i] = np.sqrt(dx*dx + dy*dy) return elementlengder |
Lesinput
| Code Block | ||||||
|---|---|---|---|---|---|---|
| ||||||
def lesinput(): # Åpner inputfilenkolonne 3; 1 = fast innspent og 0 = fri rotasjon punktfid = np.loadtxt(fid, dtype = int, max_rows = npunkt) # 'max_rows = npunkt' sorger for at vi bare leser # de 'npunkt' neste linjene i tekstfilen # Leser antall elementer nelem = int(fid.readline()) open("input.txt", "r") # Leser totalt antall punkt npunkt = int(fid.readline()) # 'fid.readline()' leser en linje, 'int(...)' gjør at linjen tolkes som et heltall # LESER INN XY-KOORDINATER TIL KNUTEPUNKTENE # LeserNodenummer konnektivitet:tilsvarer Sammenhengradnummer mellom elementender og knutepunktnummer samt EI for elementene # Elementnummer tilsvarer radnummer i "elem"-variabeli "Node-variabel" # Knutepunktnummerx-koordinat forlagres lokali endekolonne 1, lagresy-koordinat i kolonne 12 # Knutepunktnummer for lokal ende 2 Grensebetingelse lagres i kolonne 2 3; 1 = fast #innspent Detog anbefales0 at= knutepunktnummereringfri starterrotasjon på 0, slik atpunkt det samsvarerer med listeindeksering i Python # E-modul for materiale lagres i kolonne 3 # Tverrsnittstype lagres i kolonne 4; I-profil = 1 og rørprofil = 2 elem = np.loadtxt(fid, dtype = int, max_rows = nelem= np.loadtxt(fid, dtype = int, max_rows = npunkt) # 'max_rows = npunkt' sorger for at vi bare leser # de 'npunkt' neste linjene i tekstfilen # Leser antall elementer nelem = int(fid.readline()) # Leser antallkonnektivitet: lasterSammenheng sommellom virkerelementender påog rammen knutepunktnummer samt EI nlast = int(fid.readline()) for elementene # Leser lastdataElementnummer tilsvarer radnummer i "elem"-variabel # BestemKnutepunktnummer selvfor verdienelokal somende er1 nødvendiglagres åi lesekolonne inn,1 samt hva verdiene som# lesesKnutepunktnummer innfor skallokal representere ende 2 lagres i lastkolonne = np.loadtxt(fid, dtype = float, max_rows = nlast) 2 # Det anbefales at knutepunktnummerering starter på 0, slik at det samsvarerer med listeindeksering i Python # <E--modul for Forslagmateriale tillagres innlesingi avkolonne lastdata 3 # Lukker input-filen fid.close() return npunkt, punkt, nelem, elem, nlast, lastTverrsnittstype lagres i kolonne 4; I-profil = 1 og rørprofil = 2 elem = np.loadtxt(fid, dtype = int, max_rows = nelem) # Leser antall laster som virker på rammen nlast = int(fid.readline()) # Leser lastdata # Bestem selv verdiene som er nødvendig å lese inn, samt hva verdiene som leses inn skal representere last = np.loadtxt(fid, dtype = float, max_rows = nlast) # <-- Forslag til innlesing av lastdata # Lukker input-filen fid.close() return npunkt, punkt, nelem, elem, nlast, last |
Visualisering
Den tildelte visualiseringdelen av prosjektet består av en fil (structure_visualization) med flere funksjoner som er laget for plotte rammen både før og etter lastene er påført. Disse plottene er i hovedsak ikke lagt for å være veldig pene, men det anbefales heller at de blir brukt til å verifisere at input-filen har blitt tolket riktig, samt at de endelige rotasjonene er noenlunde realistiske.
| Info |
|---|
Det er ikke et krav å bruke disse funksjonene dersom det ikke ønskes, og det er heller ikke et krav å forstå hva som blir gjort inni funksjonene dersom man velger å bruke dem. De er kun laget for å være til hjelp ved å tilby en simpel visualisering av problemet (og delvis løsningen), da dette av erfaring kan være et motiverende og effektivt verktøy. |
| Info |
|---|
Filen under kan også lastes ned ved å trykke på denne linken: structure_visualization.py |
| Code Block | ||||||
|---|---|---|---|---|---|---|
| ||||||
# This visualization is based on the original Matlab code by Josef Kiendl, and has been modified to fit TMR4176 by Jon Arnt Kårstad
# NB! Denne filen krever at du har installert Python-pakkene: NumPy, SciPy og Matplotlib
# More detailed information regarding Python (matrices, visualizations etc.) may be found at 'https://www.ntnu.no/wiki/display/imtsoftware'
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import CubicSpline
def setup_plots():
fig_init, ax_init = plt.subplots()
fig_def, ax_def = plt.subplots()
ax_init.set_title('Initialramme')
ax_def.set_title('Deformert ramme')
ax_init.axes.set_aspect('equal')
ax_def.axes.set_aspect('equal')
return fig_init, ax_init, fig_def, ax_def
def plot_structure(ax, punkt, elem, numbers, index_start):
# This is a translation of the original function written by Josef Kiendl in Matlab
# It has been slightly modified in order to be used in TMR4176
# This function plots the beam structure defined by nodes and elements
# The bool (0 or 1) 'numbers' decides if node and element numbers are plotted or not
# Change input to the correct format
nodes = np.array(punkt[:, 0:2], copy = 1, dtype = int)
el_nod = np.array(elem[:, 0:2], copy=1, dtype=int) + 1
# Start plotting part
for iel in range(0, el_nod.shape[0]):
# Plot element
ax.plot([nodes[el_nod[iel, 0] - 1, 0], nodes[el_nod[iel, 1] - 1, 0]],
[nodes[el_nod[iel, 0] - 1, 1], nodes[el_nod[iel, 1] - 1, 1]], '-k', linewidth = 2)
if numbers == 1:
# Plot element numbers. These are not plotted in the midpoint to
# avoid number superposition when elements cross in the middle
ax.text(nodes[el_nod[iel, 0] - 1, 0] + ( nodes[el_nod[iel, 1] - 1, 0] - nodes[el_nod[iel, 0] - 1, 0] ) / 2.5,
nodes[el_nod[iel, 0] - 1, 1] + ( nodes[el_nod[iel, 1] - 1, 1] - nodes[el_nod[iel, 0] - 1, 1] ) / 2.5,
str(iel + index_start), color = 'blue', fontsize = 16)
if numbers == 1:
# Plot node number
for inod in range(0, nodes.shape[0]):
ax.text(nodes[inod, 0], nodes[inod, 1], str(inod + index_start), color = 'red', fontsize = 16)
def plot_structure_def(ax, punkt, elem, numbers, index_start, r):
# This is a translation of the original function written by Josef Kiendl in Matlab
# This function plots the deformed beam structure defined by nodes and elements
# The bool (0 or 1) 'numbers' decides if node and element numbers are plotted or not
# Change input to the correct format
nodes = np.array(punkt[:, 0:2], copy = 1, dtype = int)
el_nod = np.array(elem[:, 0:2], copy=1, dtype=int) + 1
nod_dof = np.arange(1, nodes.shape[0] + 1, 1, dtype=int)
if numbers == 1:
# Plot node number
for inod in range(0, nodes.shape[0]):
ax.text(nodes[inod, 0], nodes[inod, 1], str(inod + index_start), color = 'red', fontsize = 16)
for iel in range(0, el_nod.shape[0]):
delta_x = nodes[el_nod[iel, 1] - 1, 0] - nodes[el_nod[iel, 0] - 1, 0]
delta_z = nodes[el_nod[iel, 1] - 1, 1] - nodes[el_nod[iel, 0] - 1, 1]
L = np.sqrt(delta_x ** 2 + delta_z ** 2)
if delta_z >= 0:
psi = np.arccos(delta_x / L)
else:
psi = -np.arccos(delta_x / L)
phi = np.zeros((2, 1))
for inod in range(0, 2):
if nod_dof[el_nod[iel, inod] - 1] > 0:
phi[inod] = r[nod_dof[el_nod[iel, inod] - 1] - 1]
x = np.array([0, L])
z = np.array([0, 0])
xx = np.arange(0, 1.01, 0.01)*L
cs = CubicSpline(x, z, bc_type = ((1, -phi[0, 0]), (1, -phi[1, 0])))
zz = cs(xx)
# Rotate
xxzz = np.array([[np.cos(psi), -np.sin(psi)], [np.sin(psi), np.cos(psi)]]) @ np.vstack([xx, zz])
# Displace
xx2 = xxzz[0, :] + nodes[el_nod[iel, 0] - 1, 0]
zz2 = xxzz[1, :] + nodes[el_nod[iel, 0] - 1, 1]
ax.plot(xx2, zz2, '-k', linewidth = 2)
if numbers == 1:
# Plot element numbers. These are not plotted in the midpoint to
# avoid number superposition when elements cross in the middle
ax.text(xx2[round(xx2.size / 2.5)], zz2[round(xx2.size / 2.5)], str(iel + index_start), color = 'blue', fontsize = 16) |
Eksempel på bruk av visualisering som et verktøy
La oss se på et lite eksempel som viser hvordan visualisering kan hjelpe oss undervegs i implementeringen av matrise-metoden.
Etter å ha diskretisert eksempelproblemet til noder og elementer, ført det inn i en input-fil og til slutt prøvekjørt funksjonen lesinput er det gjerne ønskelig å verifisere det midlertidige resultatet. Ved å kjøre den første visualiseringsfunksjonen plot_structure kan vi enkelt se om rammen er slik vi på forhånd har definert den.
Basert på bildet over kan vi kjapt fastslå at alle noder og elementer både henger sammen og er riktig nummerert. Hadde noe sett feil ut ville vi ha måttet gjøre endringer i input-filen. Etter å ha implementert resten av matrisen-metoden har vi også tilgang på rotasjonene, og vi kan nå gjøre samme prosedyre for den deformerte rammen ved å kjøre funksjonen plot_structure_def.
Her er skalering satt til 100 (dette betyr at alle deformasjoner er skalert opp 100 ganger), og vi kan tolke dette som potensielt realistisk resultat (det er ikke et poeng at deformasjonene skal være av helt samme grad som den jeg har her, men derimot at de ikke skal være åpenbart feil). La oss derimot se på hva som skjer dersom vi endrer f.eks. E-modulen til stål til en verdi som er altfor liten, og dermed feil.
Her kan vi se at noe absolutt ikke er som det burde være. La oss sette skalering = 1 og se om dette endrer på noe.
Vi kan nå se elementene litt mer tydelig, men det er fortsatt ganske klart at noe er feil. Når vi har kommet så langt at koden beregner disse rotasjonene (og tilhørende endemoment), kan et tips være å sjekke hvilken bøyestivhet vi har brukt for bjelkene og rørene våre (E-modul, 'E', og arealtreghetsmoment, 'I'). Disse verdiene har stor påvirkning på stivhetsmatrisen, 'K', og dermed også stor påvirkning på de endelige rotasjonene.
Inputfil - Eksempel
Følgende er et eksempel på hvordan input.txt kan se se ut. I dette eksemplet er det ikke lagt til noen laster på konstruksjonen.
| Note |
|---|
OBS! Husk å fjerne kommentarene fra egen inputfil. |
| No Format |
|---|
9 # Antall knutepunkt 0 0 1 # [x, y, Innspenning] 0 10 0 0 20 0 10 20 0 10 10 0 10 0 0 20 20 0 20 10 0 20 0 0 10 # Antall element 0 1 70854000 1 # [Lokal ende 1, Lokal ende 2, Elastisitetsmodul, Profil] 1 2 70854000 1 1 4 141708000 2 2 3 70854000 2 4 3 70854000 1 5 4 70854000 1 3 6 70854000 2 4 7 141708000 1 7 6 70854000 2 8 7 70854000 1 0 # Antall laster |