You are viewing an old version of this page. View the current version.

Compare with Current View Page History

« Previous Version 56 Next »

Prosjektoppgaven i TMR4167 innebærer å analysere en konstruksjon med matrisemetoden. Hensikten med prosjektet er å få trening i å utføre beregningsoppgaver med Python, samt få innsikt i oppbyggingen av dataprogram for analyse av rammekonstruksjoner.

Koden vi har presentert under krever at følgende Python-pakker er installert: NumPy, SciPy og Matplotlib. Mer informasjon om installering av pakker er tilgjengelig på siden Python.

Page content

'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, lengder og all visualisering allerede er implementert. Flere detaljer om disse funksjonene kommer i neste avsnitt.

Rammeanalyse
from structure_visualization import *

# -----Rammeanalyse-----
def main():
    # -----Initialiserer figurer-----
    fig_init, ax_init, fig_def, ax_def = setup_plots()
    
    # -----Til visualiseringen, velg første indeks brukt i nummerering av noder og element-----
    first_index = 1

    # -----Leser input-data-----
    npunkt, punkt, nelem, elem, nlast, last = lesinput()

	# -----Plott initalramme-----
    plot_structure(fig_init, ax_init, punkt, elem, 1, first_index)

    # -----Regner ut lengder til elementene------
    elementlengder = lengder(punkt, elem, nelem)

    # -----Fastinnspenningsmomentene------
    # Lag funksjon selv
    fim = moment(npunkt, punkt, nelem, elem, nlast, last, elementlengder)

    # -----Setter opp lastvektor-----
    # Lag funksjon selv
    b = lastvektor(fim, npunkt, punkt, nelem, elem)

    # ------Setter opp systemstivhetsmatrisen-----
    # Lag funksjon selv
    K = stivhet(nelem, elem, elementlengder, npunkt)

    # ------Innfører randbetingelser------
    # Lag funksjon selv
    Kn, Bn = bc(npunkt, punkt, K, b)

    # -----Løser ligningssystemet------
	# Lag funksjon selv
    rot = ...
	# Hint, se side for løsing av lineære systemer i Python
    
    #------Finner endemoment for hvert element-----
    # Lag funksjon selv
    endemoment = endeM(npunkt, punkt, nelem, elem, elementlengder, rot, fim)

    #-----Skriver ut hva rotasjonen ble i de forskjellige nodene-----
    print("Rotasjoner i de ulike punktene:")
    print(rot)

    #-----Skriver ut hva momentene ble for de forskjellige elementene-----
    print("Elementvis endemoment:")
    print(endemoment)

	#-----Plott deformert ramme-----
    skalering = 100;     # Du kan endre denne konstanten for å skalere de synlige deformasjonene til rammen
    plot_structure_def(fig_def, ax_def, punkt, elem, 1, first_index, skalering*rot)
    plt.show()

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

Elementlengder
def lengder(knutepunkt, element, nelem):

    elementlengder = np.zeros((nelem, 1))
    # Beregner elementlengder med Pythagoras' laeresetning
    for i in range (0, nelem):
        # 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]
        dy = knutepunkt[element[i, 0], 1] - knutepunkt[element[i, 1], 1]
        elementlengder[i] = np.sqrt(dx*dx + dy*dy)

    return elementlengder

Lesinput

Les input
def lesinput():

    # Åpner inputfilen
    fid = 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
    # Nodenummer tilsvarer radnummer i "Node-variabel"
    # x-koordinat lagres i kolonne 1, y-koordinat i kolonne 2
    # Grensebetingelse lagres i kolonne 3; 1 = fast innspent og 0 = fri rotasjon
    punkt = 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 konnektivitet: Sammenheng mellom elementender og knutepunktnummer samt EI for elementene
    # Elementnummer tilsvarer radnummer i "elem"-variabel
    # Knutepunktnummer for lokal ende 1 lagres i kolonne 1
    # Knutepunktnummer for lokal ende 2 lagres i kolonne 2
    # Det anbefales at knutepunktnummerering starter på 0, slik at 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)

    # 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.

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 er en mer motiverende måte å debugge på.


structure_visualization.py
# 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(fig, 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)
    fixed_dof = np.array(np.where(punkt[:, 2] == 1)[0], copy=1, dtype=int) + 1
    el_nod = np.array(elem[:, 0:2], copy=1, dtype=int) + 1
    nod_dof = np.arange(1, nodes.shape[0] + 1, 1, dtype=int)

    # 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(fig, 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)
    fixed_dof = np.array(np.where(punkt[:, 2] == 1)[0], copy=1, dtype=int) + 1
    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)

Inputfil - Eksempel

Følgende er et eksempel på hvordan input.txt kan se ut. I dette eksemplet er det ikke lagt til noen laster på konstruksjonen. 

OBS! Husk å fjerne kommentarene fra egen inputfil.


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
  • No labels