'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-----
	# Endre gjerne selv
    first_index = 0

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

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

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

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

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

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

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

    # -----Løser ligningssystemet------
	# Lag funksjonen selv
    rot = ...
	# Hint, se side for løsing av lineære systemer i Python
    
    #------Finner endemoment for hvert element-----
    # Lag funksjonen 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(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' læresetning
    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 kan være et motiverende og effektivt verktøy.

Filen under kan også lastes ned ved å trykke på denne linken: structure_visualization.py

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(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 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

Page tree

Prosjektoppgave i TMR4167 Marin teknikk - Konstruksjoner