[PYTHON] Der Beginn der Finite-Elemente-Methode (eindimensionale Elementsteifigkeitsmatrix)

Über diesen Artikel

Dieser Artikel wurde in Java verfasst und bezieht sich auf Erste Finite-Elemente-Methode von Java, veröffentlicht von Corona. Es ist ein Programm, das mit Python erstellt wurde. Außerdem ist es diesmal Sache des Programms, die Elementsteifigkeitsmatrix zu erstellen.

Überblick über die Finite-Elemente-Methode

Die Finite-Elemente-Methode wird zur Simulation der Spannung und Dehnung der Festkörperdynamik verwendet. In der Festkörpermechanik wird ein Objekt in viele kleine Bereiche (Elemente) unterteilt, und der Zustand jedes Elements wird berechnet, um den Verformungszustand der gesamten Substanz zu bestimmen. Der Zustand eines Objekts wird durch drei Differentialgleichungen ausgedrückt (parallele Spannungsgleichung, relationaler Dehnungsverschiebungsausdruck, relationaler Spannungs-Dehnungs-Ausdruck (Konstruktionsausdruck)). Bei der Finite-Elemente-Methode werden diese Differentialgleichungen jedoch als äquivalente Integralausdrücke ausgedrückt ( Durch Ersetzen durch das Prinzip der virtuellen Arbeit werden die Verschiebung und Dehnung ungefähr erhalten. Die integrale Darstellung ersetzt die Differentialgleichungen durch die starren Gleichungen, die simultane lineare Gleichungen sind. Darüber hinaus besteht der Berechnungsablauf der Finite-Elemente-Methode darin, das Objekt in jedes Element zu unterteilen, eine Steifigkeitsgleichung (Elementsteifigkeitsmatrix) für jedes Element zu erstellen und dann die Steifigkeitsgleichungen aller Elemente zu integrieren, um die Gesamtsteifigkeitsgleichung (insgesamt) zu erstellen. Starrheitsmatrix) wird zusammengesetzt und eine ungefähre Lösung berechnet.

Grundlegende Differentialgleichung

Der Einfachheit halber ist ein Endpunkt eines Objekts mit einer Querschnittsfläche von $ A \ \ rm {[L ^ {2}]} $ und einer Länge von $ l \ \ rm {[L]} $ festgelegt und der andere Endpunkt Betrachten wir den Fall, in dem eine eindimensionale und gleichmäßige Zugbelastung $ F \ \ rm {[F]} $ wirkt. (Hier sei $ \ rm {[L]} $ die Dimension der Länge und $ \ rm [F] $ die Dimension der Kraft. Mit einer solchen Notation ist beispielsweise Druck eine Flächeneinheit. Da es sich um eine Trefferkraft handelt, kann sie als Kraftdimension geteilt durch das Quadrat der Längendimension geschrieben werden, dh $ \ rm {[FL ^ {-2}]} $. Die Parameterdimension ist klar definiert. Dadurch können Sie leichter erkennen, was dieser Parameter darstellt und wie er sich auf andere Parameter bezieht.) Außerdem wird diese Substanz zu einem elastischen Körper verarbeitet (eine Substanz, die sich bei Anwendung einer Kraft verformt und beim Stoppen der Anwendung in ihre ursprüngliche Form zurückkehrt), und ihr Young-Verhältnis beträgt $ E \ \ rm {[FL ^ {-2}]} Angenommen, es ist $. Das Youngsche Verhältnis ist der proportionale Koeffizient, der auf die Dehnung $ \ varepsilon $ [Dimension] des Hookschen Gesetzes angewendet wird, dass "Spannung proportional zur Dehnung ist", die in einem elastischen Körper gilt. (Ähnlich wie die Beziehung zwischen Last (oder Spannung) und Verschiebung, wenn eine Kraft auf eine Feder ausgeübt wird, aber beachten Sie, dass die Abmessungen unterschiedlich sind.)

Unter dieser Annahme kann das Kräftegleichgewicht in diesem Objekt wie folgt geschrieben werden, vorausgesetzt, die im Objekt erzeugte Spannung beträgt $ \ sigma \ \ rm {[FL ^ {-2}]} $.

$ F = \ sigma A $ (Kraftausgleich, Spannungsausgleichsgleichung) (1 ')

Unter der Annahme, dass der Stress $ \ sigma $ ist, kann das Hooksche Gesetz wie folgt geschrieben werden, indem der Young-Satz $ E $ und der Stamm $ \ varepsilon $ verwendet werden.

$ \ Sigma = E \ \ varepsilon $ (Materialeigenschaften, Spannungs-Dehnungs-Gleichung) (2 ')

Die Dehnung $ \ varepsilon $ kann wie folgt ausgedrückt werden, wobei die Dehnung des Objekts unter Zugbelastung $ u \ rm {[L]} $ beträgt.

$ \ Varepsilon = u / L $ (Formänderung, Dehnungs-Verschiebungs-Beziehung) (3 ')

Diese Gleichungen (1 ') - (3') sind Sonderfälle grundlegender Differentialgleichungen. In Anbetracht der allgemeineren Situation werden die folgenden drei grundlegenden Differentialgleichungen abgeleitet.

$ d (\ sigma A) / dx = 0 $ (Kraftausgleich, Spannungsausgleichsgleichung) (1) $ \ Sigma = E \ \ varepsilon $ (Materialeigenschaften, Spannungs-Dehnungs-Gleichung) (2) $ \ Varepsilon = du / dx $ (Formänderung, relationaler Ausdruck von Dehnung und Verschiebung) (3)

Bei der analytischen Berechnung der in einem elastischen Körper erzeugten Dehnung müssen die Differentialgleichungen der Gleichungen (1) - (3) gelöst werden.

Prinzip der virtuellen Arbeit

Das Prinzip der virtuellen Arbeit wird wie folgt erklärt. Bitte beachten Sie die Referenzen, da die ausführliche Erklärung lang sein wird.

Wenn das Objekt deformiert und in einem ausgeglichenen Zustand ist, entspricht die von der externen Kraft geleistete virtuelle Arbeit der von der internen Kraft geleisteten virtuellen Arbeit.

Das Prinzip der virtuellen Arbeit kann durch eine mathematische Formel wie folgt ausgedrückt werden.

Über willkürliche virtuelle Verschiebung und virtuelle Belastung bei ausgeglichener Kraft         \int_V \sigma \delta\varepsilon dV = F\delta u Ist festgelegt.

Wobei $ \ delta \ varepsilon $ die virtuelle Belastung und $ \ delta u $ die virtuelle Verschiebung ist.

Begrenzte Elementmethode bei einer Dimension (bis zur Herstellung der Elementsteifigkeitsmatrix)

Elementteilung, Prinzip der virtuellen Arbeit

Jetzt ist die Finite-Elemente-Methode fertig. Bei der Finite-Elemente-Methode wird das Objekt zunächst in kleine Bereiche, "Elemente", unterteilt. Wenden Sie dann das Prinzip der virtuellen Arbeit auf jedes Element an und erstellen Sie für jedes Element eine Elementsteifigkeitsgleichung (Elementsteifigkeitsmatrix). Betrachten wir das Prinzip der virtuellen Arbeit für jedes Element, wenn ein eindimensionales Objekt geteilt wird. Da ein Element zwei Knoten (Verbindungen zwischen Elementen) hat, ist einer Knoten 1 (Koordinate $ x_1 $) und der andere Knoten 2 (Koordinate $ x_2 $), und die auf jeden wirkende innere Kraft ist $ f_1. Die Punktverschiebungen seien $ u_1 $ bzw. $ u_2 $ für $ und $ f_2 $. Wenn die Spannung hier $ \ sigma $ ist und die virtuelle Verschiebung jedes Knotens $ \ delta u_1 $, $ \ delta u_2 $ ist, ist das Prinzip der virtuellen Arbeit dieses Elements

$\int_{x_1}^{x_2} \sigma \delta\varepsilon A dx = f_1 \cdot \delta u_1 + f_2 \cdot \delta u_2 $  (4)

Kann geschrieben werden.

Lineare Annäherung der Verschiebung

Nehmen wir nun an, dass sich in jedem Element die Verschiebung $ u $ mit einer linearen Funktion der Koordinaten $ x $ ändert.

u(x) = ax+b  (5)

Diese Verschiebung muss also gleich $ u_1 $, $ u_2 $ an den Knoten 1 und 2 sein

u(x_1) = ax_1+b =u_1  (6)         u(x_1) = ax_2+b =u_2  (7)

Zusammenfassend für $ a $ und $ b $,

a=\frac{u_2 - u_1}{x_2 - x_1} = \frac{u_2 - u_1}{l}  (8)

b=\frac{x_2 u_1 - x_1 u_2}{x_2 - x_1}=\frac{x_2 u_1 - x_1 u_2}{l}  (9)

Es wird sein. Wobei $ l = x_2-x_1 $ die Länge des Elements ist. Mit (8) und (9) ist (5)

u(x) = ax+b = \frac{u_2 - u_1}{l}x+\frac{x_2 u_1 - x_1 u_2}{l}  (10)

Durch Umschreiben von Gleichung (10)

u(x) = \frac{x_2 - x}{l}u_1+\frac{x - x_1}{l} u_2 =N_1(x)u_1 + N_2(x)u_2  (11)

Daher war es möglich, die Gleichung mit den Knotenverschiebungen $ u_1 $ und $ u_2 $ als unbekannte Zahlen umzuschreiben. Außerdem werden $ N_1 (x) $ und $ N_2 (x) $ als Formfunktionen bezeichnet. Auch kann Gleichung (11) in einer Matrix ausgedrückt werden,

$u(x) = {\bf N}^{{\rm T}} {\bf U} $  (12)

Kann geschrieben werden. Wo $ {\ bf N} $ und $ {\ bf U} $ sind


  {\bf N} = \left(
    \begin{array}{c}
      N_1(x) \\
      N_2(x) 
    \end{array}
  \right),

  {\bf U} = \left(
    \begin{array}{c}
      u_1 \\
      u_2
    \end{array}
  \right)  (13)

Es ist ein Vektor.

Ungefähre Anzeige von Belastung und Stress

Unter Verwendung der Verschiebungsgleichung (11) und der Dehnungs-Verschiebungs-Beziehungsgleichung (3) ist die Dehnung $ \ varepsilon $

\varepsilon = \frac{d}{dx}(N_1(x)u_1 + N_2(x)u_2) = \frac{dN_1(x)}{dx}u_1 + \frac{dN_2(x)}{dx}u_2 (14)

Kann geschrieben werden. Dies ist auch eine Matrixnotation,

$\varepsilon = {\bf B}^{{\rm T}} {\bf U} $  (15)

Kann umgeschrieben werden als. Hier wird $ {\ bf B} $ als $ {\ bf B} $ -Matrix oder Dehnungsverschiebungsmatrix bezeichnet und kann wie folgt geschrieben werden.

  {\bf B} = \left(
    \begin{array}{c}
      \frac{dN_1(x)}{dx} \\
      \frac{dN_2(x)}{dx} 
    \end{array}
  \right)
  =
   \left(
    \begin{array}{c}
      \frac{-1}{l} \\
      \frac{1}{l} 
    \end{array}
  \right)
 (16)

Die Spannung $ \ sigma $ wird aus den Spannungs-Dehnungs-Beziehungsgleichungen (2) und (15) berechnet.

$\sigma = E \varepsilon = E{\bf B}^{{\rm T}} {\bf U} $  (17)

Kann ausgedrückt werden als.

Virtuelle Verschiebung und Belastung

Aus den Gleichungen (12) und (15) werden die virtuelle Verschiebung $ \ delta u $ und die virtuelle Dehnung $ \ delta \ varepsilon $ in Matrixnotation ausgedrückt.

$\delta u = {\bf N}^{{\rm T}} \delta{\bf U} $  (18)         $\delta \varepsilon = {\bf B}^{{\rm T}} \delta{\bf U} $  (19)

Kann geschrieben werden als.

Elementsteifigkeitsgleichung (Elementsteifigkeitsmatrix)

Betrachten wir nun Gleichung (4) des Prinzips der virtuellen Arbeit.

$\int_{x_1}^{x_2} \sigma \delta\varepsilon A dx = f_1 \cdot \delta u_1 + f_2 \cdot \delta u_2 $  (4)

Schreiben Sie zunächst das Innere des Integrals auf der linken Seite von Gleichung (4) mit Matrixnotation um.

\sigma \delta \varepsilon = E{\bf B}^{{\rm T}} {\bf U}{\bf B}^{{\rm T}} \delta{\bf U} = E {\bf U}^{{\rm T}}{\bf B}{\bf B}^{{\rm T}} \delta{\bf U}  (20)

Wenn dies auf die linke Seite von Gleichung (4) angewendet wird,

$\int_{x_1}^{x_2} \sigma \delta\varepsilon A dx = \int_{x_1}^{x_2} E {\bf U}^{{\rm T}}{\bf B}{\bf B}^{{\rm T}} \delta{\bf U}A dx $  (21)

In Anbetracht der Tatsache, dass die Querschnittsfläche $ A $ und die junge Rate $ E $ in den Elementen in Gleichung (21) konstant sind,

$\int_{x_1}^{x_2} \sigma \delta\varepsilon A dx = E {\bf U}^{{\rm T}}{\bf B}{\bf B}^{{\rm T}} \delta{\bf U}A \int_{x_1}^{x_2} dx = E V {\bf U}^{{\rm T}}{\bf B}{\bf B}^{{\rm T}} \delta{\bf U} $  (22)

Dabei ist $ V $ $ A \ int_ {x_1} ^ {x_2} dx = El = V $, was das Volumen des Elements darstellt. Wenn Sie $ {\ bf K} = E V {\ bf B} {\ bf B} ^ {{\ rm T}} $ setzen, können Sie es wie folgt umschreiben.

$\int_{x_1}^{x_2} \sigma \delta\varepsilon A dx = {\bf U}^{{\rm T}} {\bf K} \delta{\bf U} $  (23)

$ {\ bf K} $ heißt die Elementsteifigkeitsmatrix.

Wenn als nächstes die rechte Seite von Gleichung (4) in einer Matrix ausgedrückt wird,

f_1 \cdot \delta u_1 + f_2 \cdot \delta u_2 = {\bf F}^{{\rm T}} \delta{\bf U}  (24)

Wo $ {\ bf F} $ ist:

  {\bf F} = \left(
    \begin{array}{c}
      f_1 \\
      f_2 
    \end{array}
  \right)
 (25)

Auf diese Weise kann Gleichung (4) des Prinzips der virtuellen Arbeit unter Verwendung der Gleichungen (23) und (24) wie folgt umgeschrieben werden.

{\bf U}^{{\rm T}} {\bf K} \delta{\bf U} = {\bf F}^{{\rm T}} \delta{\bf U}  (26)

Wenn die obige Formel durch $ \ delta {\ bf U} $ zusammengefasst wird,

({\bf U}^{{\rm T}} {\bf K} - {\bf F}^{{\rm T}} )\delta{\bf U} = {\bf 0}  (27)

Da hier das Prinzip der virtuellen Arbeit für jede virtuelle Verschiebung $ \ delta {\ bf U} $ gilt,

{\bf U}^{{\rm T}} {\bf K} = {\bf F}^{{\rm T}}  (28)

Wenn Sie diese Umkehrung nehmen,

$ {\bf K}{\bf U} = {\bf F}$    (29)

Es wird sein. Gleichung (29) ist eine Gleichung, die als Elementsteifigkeitsgleichung bezeichnet wird. Wie eingangs erwähnt, werden nach der Formulierung dieser Elementsteifigkeitsgleichung für jedes Element die Steifigkeitsgleichungen aller Elemente integriert, um die Gesamtsteifigkeitsgleichung (Gesamtsteifigkeitsmatrix) zu erstellen und die ungefähre Lösung zu berechnen. Es wird sein.

Berechnungsprogramm der Elementsteifigkeitsgleichung (Elementsteifigkeitsmatrix) eines eindimensionalen Elements

Nun, ich habe ein Programm erstellt, um eine eindimensionale Elementsteifigkeitsgleichung mit Python zu berechnen. Die verwendete Bibliothek ist numpy und der Code ist unten. Aus den gegebenen Bedingungen können Sie die $ {\ bf B} $ -Matrix und die Elementsteifigkeitsmatrix $ {\ bf K} $ berechnen, Spannung für ein Element.


import numpy as np

class Element:
    #Anzahl der Kontakte
    NumberOfNodeInElement = 2
    
    def __init__(self):
        #Elementlänge
        self.length = 0
        #Querschnittsfläche des Elements
        self.area = 0
        #Elastizitätsmodul
        self.young = 0
        #Kontaktnummern, aus denen das Element besteht
        self.ix = np.zeros(self.NumberOfNodeInElement)
        #Koordinatenwerte der Kontakte, aus denen das Element besteht
        self.xe = np.zeros(self.NumberOfNodeInElement)
        #Kontaktverschiebung des Elements
        self.ue = np.zeros(self.NumberOfNodeInElement)
        #Bmatrix
        self.Bmatrix = np.zeros(self.NumberOfNodeInElement)
        #Elementsteifigkeitsmatrix
        self.Kmatrix = np.zeros((self.NumberOfNodeInElement,self.NumberOfNodeInElement))
    
    #Bmatrix-Berechnung
    def makeBmatrix(self):
        self.length = self.xe[1]-self.xe[0]
        self.Bmatrix = np.array([[-1/self.length, 1/self.length]])

    #Berechnung der Elementsteifigkeitsmatrix
    def makeElementStiffness(self):
        self.volume = self.area*self.length
        self.Kmatrix = self.young*self.volume*np.dot(self.Bmatrix.T,self.Bmatrix)

    #Berechnung von Stress
    def stressCalculation(self):
        return self.young*np.vdot(self.Bmatrix,self.ue)

Beispielsweise wird die Berechnung für ein Element durchgeführt, dessen Koordinaten von Knoten 1 und Knoten 2 $ x_1 = 50 $ (mm), $ x_2 = 150 $ (mm), eine Querschnittsfläche von 100 $ mm ^ 2 $ und eine Young-Rate von 200 GPa sind. Wann

Ele = Element()
Ele.xe = np.array([50,150])
Ele.area = 100
Ele.young = 200000
Ele.ue = np.array([0.01,0.025])

Ele.makeBmatrix()
Ele.makeElementStiffness()
scal = Ele.stressCalculation()

Das Berechnungsergebnis ist in diesem Fall wie folgt.

#Bmatrix
[[-0.01  0.01]]
#Bmatrix
[[ 200000. -200000.]
 [-200000.  200000.]]
#Stress
30.000000000000004

Zusammenfassung

Dieses Mal erklärte ich bis zur eindimensionalen Elementsteifigkeitsmatrix der Finite-Elemente-Methode und stellte das Herstellungsprogramm für die Elementsteifigkeitsmatrix vor. Das nächste Mal werden wir die Gesamtsteifigkeitsgleichung (Gesamtsteifigkeitsmatrix) zusammenstellen und berechnen.

Nachschlagewerk

Shigeru Nagagi (2010) "Die erste Finite-Elemente-Methode von Java" Corona

Recommended Posts

Der Beginn der Finite-Elemente-Methode (eindimensionale Elementsteifigkeitsmatrix)
Begrenzte Elementmethode Anfang 2 (Erstellung einer eindimensionalen Gesamtsteifigkeitsmatrix)
Lösen einer eindimensionalen Wellengleichung mit der Differenzmethode (Python)
Verwenden Sie IvyFEM (Finite-Elemente-Methodenbibliothek für .NET) aus Python