Dreidimensionale Skelettstrukturanalyse mit Python

Große Lösung

Überblick

Ich habe ein dreidimensionales Programm zur Analyse der Skelettstruktur erstellt, das ich mir in Python schon immer gewünscht habe. In dem zuvor erstellten ebenen Skelettanalyseprogramm werden die Steifigkeitsmatrix und die Koordinatenumwandlungsmatrix hauptsächlich geändert, um Balkenelemente mit 1 Knoten und 6 Freiheitsgraden zu unterstützen.

~~ Der Test ist noch nicht genug, aber ~~ Ich wurde von dem Artikel unten inspiriert und habe ihn gepostet. http://qiita.com/sasaco/items/917429a497c6e9a08401

~~ Da das Programm auf Jupyter erstellt wurde, wird dieser Artikel auch gemäß der Beschreibung meines Jupyter beschrieben. Das Analyseprogramm von Python ist in drei Teile gegliedert. ~~ Ich mache mir immer Sorgen, aber der Dateiexportteil ist zu überladen. Kannst du nicht besser schreiben? .. ..

~~ Die Berechnung läuft. .. .. ~~

Die Programmierumgebung ist Machine: MacBook Pro (Retina, 13-inch, Mid 2014) OS: macOS Sierra Python version: Python 3.6.1

Theorie der dreidimensionalen Skelettstrukturanalyse

Grundlegende Angelegenheiten

Die folgenden Bedingungen sind in der Theorieentwicklung und im Programm enthalten.

In diesem Programm wird unter der Annahme, dass das Modell großräumig sein wird, eine spärliche Matrixverarbeitung eingeführt, wenn simultane lineare Gleichungen gelöst werden.

Elementsteifigkeitsgleichung

Die allgemeine Form der Elementsteifigkeitsgleichung ist dieselbe wie bei der zweidimensionalen Spannungsanalyse.

Gleichung zum Finden der Verschiebung

\begin{equation}
[\boldsymbol{k}]\{\boldsymbol{u_{nd}}\}=\{\boldsymbol{f}\}+\{\boldsymbol{f_t}\}+\{\boldsymbol{f_b}\}
\end{equation}

\begin{align}
&[\boldsymbol{k}]=\int_V[\boldsymbol{B}]^T[\boldsymbol{D_e}][\boldsymbol{B}]dV  & & \text{(Elementsteifigkeitsmatrix)}\\
&\{\boldsymbol{f}\}=\int_A [\boldsymbol{N}]^T\{\boldsymbol{S}\}dA  & & \text{(Nodaler externer Kraftvektor)}\\
&\{\boldsymbol{f_t}\}=\int_V[\boldsymbol{B}]^T[\boldsymbol{D_e}]\{\boldsymbol{\epsilon_0}\}dV  & & \text{(Anfänglicher Knotenverzerrungsvektor)}\\
&\{\boldsymbol{f_b}\}=\gamma\left(\int_V [\boldsymbol{N}]^T [\boldsymbol{N}]dV\right)\{\boldsymbol{w_{nd}}\} & & \text{(Knotenträgheitskraftvektor)}
\end{align}

Berechnungsformel für die interne Spannung des Elements

\begin{equation}
\{\boldsymbol{\sigma}\}=[\boldsymbol{D_e}]\{\boldsymbol{\epsilon}-\boldsymbol{\epsilon_0}\}
\end{equation}

\begin{align}
&\{\boldsymbol{\epsilon}\}  & &\text{: Dehnungsvektor des Elements, erhalten aus der Verschiebung, die durch Lösen der Steifigkeitsgleichung erhalten wurde}\\
&\{\boldsymbol{\epsilon_0}\} & &\text{: Anfangsbelastung wie Temperaturbelastung}\\
&[\boldsymbol{D_e}]          & &\text{: Spannungs-Dehnungs-Beziehungsmatrix}
\end{align}

Elementsteifigkeitsmatrix der dreidimensionalen Rahmenstruktur

Die Elementsteifigkeitsmatrix der dreidimensionalen Skelettstruktur ist eine 12 × 12-Quadratmatrix. Dies kann abgeleitet werden, indem das für das flache Skelett erweitert wird, außer um die Verdrehung zu berücksichtigen. Der Ableitungsprozess entfällt und nur die Ergebnisse werden hier beschrieben.

Elementsteifigkeitsgleichung (einfachster Ausdruck)

\begin{equation}
\boldsymbol{\{f\}}=\boldsymbol{[k]}\boldsymbol{\{u\}}
\end{equation}

Knotenverschiebungsvektor

\begin{equation}
\boldsymbol{\{u\}}=
\begin{Bmatrix}
u_i & v_i & w_i & \theta_{xi} & \theta_{yi} & \theta_{zi} & u_j & v_j & w_j & \theta_{xj} & \theta_{yj} & \theta_{zj}
\end{Bmatrix}^T
\end{equation}

\begin{align}
&u & &\text{: Verschiebung der X-Achse} & &\theta_x& &\text{: Rotationsbetrag der X-Achse (Verdrehwinkel)}\\
&v & &\text{: Verschiebung der Y-Achse} & &\theta_y& &\text{: Rotationsbetrag der Y-Achse}\\
&w & &\text{: Verschiebung der Z-Achse} & &\theta_z& &\text{: Z-Achsen-Rotationsbetrag}
\end{align}

Knotenkraftvektor

\begin{equation}
\boldsymbol{\{f\}}=
\begin{Bmatrix}
N_i & S_{yi} & S_{zi} & M_{xi} & M_{yi} & M_{zi} & N_j & S_{yj} & S_{zj} & M_{xj} & M_{yj} & M_{zj}
\end{Bmatrix}^T
\end{equation}

\begin{align}
&N  & &\text{: Axialkraft (x-Achsenrichtungskraft)} & &M_x& &\text{: X-Achsenmoment (Verdrehmoment)}\\
&S_y& &\text{: Scherkraft der Y-Achse}   & &M_y& &\text{: Biegemoment der Y-Achse}\\
&S_z& &\text{: Z-Achsen-Scherkraft}   & &M_z& &\text{: Biegemoment der Z-Achse}
\end{align}

Elementsteifigkeitsmatrix

\begin{equation}
\boldsymbol{[k]}=
\begin{bmatrix}
\cfrac{EA}{L} & 0 & 0 & 0 & 0 & 0 & -\cfrac{EA}{L} & 0 & 0 & 0 & 0 & 0 \\
0 & \cfrac{12EI_z}{L^3} & 0 & 0 & 0 & \cfrac{6EI_z}{L^2} & 0 & -\cfrac{12EI_z}{L^3} & 0 & 0 & 0 & \cfrac{6EI_z}{L^2} \\
0 & 0 & \cfrac{12EI_y}{L^3} & 0 & -\cfrac{6EI_y}{L^2} & 0 & 0 & 0 & -\cfrac{12EI_y}{L^3} & 0 & -\cfrac{6EI_y}{L^2} & 0 \\
0 & 0 & 0 & \cfrac{GJ}{L} & 0 & 0 & 0 & 0 & 0 & -\cfrac{GJ}{L} & 0 & 0 \\
0 & 0 & -\cfrac{6EI_y}{L^2} & 0 & \cfrac{4EI_y}{L} & 0 & 0 & 0 & \cfrac{6EI_y}{L^2} & 0 & \cfrac{2EI_y}{L} & 0 \\
0 & \cfrac{6EI_z}{L^2} & 0 & 0 & 0 & \cfrac{4EI_z}{L} & 0 & -\cfrac{6EI_z}{L^2} & 0 & 0 & 0 & \cfrac{2EI_z}{L} \\
-\cfrac{EA}{L} & 0 & 0 & 0 & 0 & 0 & \cfrac{EA}{L} & 0 & 0 & 0 & 0 & 0 \\
0 & -\cfrac{12EI_z}{L^3} & 0 & 0 & 0 & -\cfrac{6EI_z}{L^2} & 0 & \cfrac{12EI_z}{L^3} & 0 & 0 & 0 & -\cfrac{6EI_z}{L^2} \\
0 & 0 & -\cfrac{12EI_y}{L^3} & 0 & \cfrac{6EI_y}{L^2} & 0 & 0 & 0 & \cfrac{12EI_y}{L^3} & 0 & \cfrac{6EI_y}{L^2} & 0 \\
0 & 0 & 0 & -\cfrac{GJ}{L} & 0 & 0 & 0 & 0 & 0 & \cfrac{GJ}{L} & 0 & 0 \\
0 & 0 & -\cfrac{6EI_y}{L^2} & 0 & \cfrac{2EI_y}{L} & 0 & 0 & 0 & \cfrac{6EI_y}{L^2} & 0 & \cfrac{4EI_y}{L} & 0 \\
0 & \cfrac{6EI_z}{L^2} & 0 & 0 & 0 & \cfrac{2EI_z}{L} & 0 & -\cfrac{6EI_z}{L^2} & 0 & 0 & 0 & \cfrac{4EI_z}{L} \\
\end{bmatrix}
\end{equation}

\begin{align}
&L& &\text{: Elementlänge}        & &J& &\text{: Drehkonstante}\\
&E& &\text{: Elastizitätskoeffizient}      & &I_y& &\text{: Sekundärmoment des Querschnitts der Y-Achse}\\
&G& &\text{: Scherelastizitätskoeffizient}& &I_z& &\text{: Sekundärmoment des Querschnitts der Z-Achse}\\
&A& &\text{: Kreuzung}        & &
\end{align}

Knotenträgheitskraftvektor des ebenen Skelettelements

Betrachten Sie die Trägheitskraft als Objektkraft und geben Sie "Beschleunigung" als Knotenvektor an. Hier besteht die einfachste Methode darin, die auf das Element wirkende Trägheitskraft einfach auf die Knotenkraft zu verteilen. Die gesamte Trägheitskraft, die auf das Element wirkt, ist wie folgt als Wert im gesamten Koordinatensystem.

\begin{align}
F_{m(x)}=&\frac{\gamma\cdot A \cdot L}{g}\cdot (k_x \cdot g) \\
F_{m(y)}=&\frac{\gamma\cdot A \cdot L}{g}\cdot (k_y \cdot g) \\
F_{m(z)}=&\frac{\gamma\cdot A \cdot L}{g}\cdot (k_z \cdot g) \\
\end{align}

Daher kann die Knotenkraft $ \ {\ boldsymbol {f_ {b}} } $ aufgrund der Trägheitskraft wie folgt bestimmt werden.

\begin{equation*}
\{\boldsymbol{F_{b}}\}=\{\boldsymbol{f_{b}}\}=\frac{\gamma\cdot A \cdot L}{2}
\begin{Bmatrix}
k_x \\
k_y \\
k_z \\
0   \\
0   \\
0   \\
k_x \\
k_y \\
k_z \\
0   \\
0   \\
0
\end{Bmatrix}
\end{equation*}

Hier ist $ \ gamma $ das Volumenvolumen des Elements, $ A $ ist die Elementquerschnittsfläche, $ L $ ist die Elementlänge und $ g $ ist die Schwerkraftbeschleunigung.

Knotentemperaturlastvektor des ebenen Skelettelements

Es wird angenommen, dass die Temperaturänderung nur die Axialkraft des Elements beeinflusst. Unter der Annahme, dass der Temperaturanstieg eine positive Temperaturänderung ist, ist der Knotenkraftvektor im Elementkoordinatensystem aufgrund der Temperaturänderung wie folgt.

\begin{equation*}
\{\boldsymbol{f_{t}}\}=EA\cdot\alpha\cdot\Delta T
\begin{Bmatrix}
-1 \\
 0 \\
 0 \\
 0 \\
 0 \\
 0 \\
 1 \\
 0 \\
 0 \\
 0 \\
 0 \\
 0
\end{Bmatrix}
\end{equation*}

Hier ist $ EA $ die axiale Steifigkeit des Elements, $ \ alpha $ ist der lineare Ausdehnungskoeffizient des Elements und $ \ Delta T $ ist der Betrag der Temperaturänderung des Elements.

Durch Umwandlung in Koordinaten wird der Knotenkraftvektor im Gesamtkoordinatensystem bestimmt.

\begin{equation*}
\{\boldsymbol{F_t}\}=[\boldsymbol{T}]^T\{\boldsymbol{f_t}\}
\end{equation*}

Hier ist $ [\ boldsymbol {T}] ^ T $ die Translokationsmatrix der Koordinatentransformationsmatrix $ [\ boldsymbol {T}] $ (in diesem Fall gleich der inversen Matrix).

Knotenaußenkraftvektor des ebenen Rahmenelements

In Bezug auf den externen Kraftvektor des Knotens wird in dem hier vorgestellten Programm die äquivalente externe Kraft des Knotens direkt aus der Eingabedatendatei eingegeben, und es wird keine spezielle Berechnungsverarbeitung durchgeführt. Die Last ist die horizontale Kraft, vertikale Kraft und das Moment im gesamten Koordinatensystem.

Koordinatenumrechnungsmatrix

Die bisher erhaltenen Elementsteifigkeitsgleichungen basieren auf dem am Element festgelegten "Elementkoordinatensystem". In einer tatsächlichen Struktur sind die Elemente in verschiedene Richtungen ausgerichtet, so dass es notwendig ist, die Beziehung zu dem einheitlich definierten "gesamten Koordinatensystem" zu finden.

Wenn hier die Koordinatenkonvertierungsmatrix, die die gesamte Verschiebung des Koordinatensystems in die Verschiebung des Elementkoordinatensystems konvertiert, $ [\ boldsymbol {T}] $ ist, Die Starrheitsmatrix $ [\ boldsymbol {K}] $ im globalen Koordinatensystem und die Starrheitsmatrix $ [\ boldsymbol {k}] $ im Elementkoordinatensystem haben die folgende Beziehung.

\begin{equation}
[\boldsymbol{K}]=[\boldsymbol{T}]^T[\boldsymbol{k}][\boldsymbol{T}]
\end{equation}

Diese Beziehung kann aus der Beziehung zwischen der Elementkoordinatensystemverschiebung $ \ {\ boldsymbol {u} } $ und der Gesamtkoordinatensystemverschiebung $ \ {\ boldsymbol {U} } $ wie folgt verstanden werden. ..

\begin{equation}
\{\boldsymbol{u}\}=[\boldsymbol{T}]\{\boldsymbol{U}\} \qquad\qquad \{\boldsymbol{u}\}^T=\{\boldsymbol{U}\}^T [\boldsymbol{T}]^T
\end{equation}

Hier ist $ [\ boldsymbol {T}] $ eine Koordinatenumwandlungsmatrix, die die gesamte Verschiebung / Knotenkraft des Koordinatensystems in die Verschiebung / Knotenkraft des Elementkoordinatensystems umwandelt.

Spezifische Darstellung der Koordinatenumwandlungsmatrix

Richtungskosinus der Stabachse (x-Achse)

Die Knoten, aus denen das Element besteht, seien $ i $ und $ j $, und ihre Koordinaten seien $ (x_i, y_i, z_i) $, $ (x_j, y_j, z_j) $.

\begin{gather}
l=\cfrac{x_j-x_i}{\ell} \qquad m=\cfrac{y_j-y_i}{\ell} \qquad n=\cfrac{z_j-z_i}{\ell} \\
\ell=\sqrt{(x_j-x_i)^2+(y_j-y_i)^2+(z_j-z_i)^2}
\end{gather}

Koordinatentransformations-Submatrix

$ \ boldsymbol {[t]} $ ist eine kleine Matrix (3 x 3) der Koordinatentransformationsmatrix (12 x 12), die das gesamte Koordinatensystem in das Elementkoordinatensystem transformiert.

\begin{equation}
\boldsymbol{\{x\}}=\boldsymbol{[t]}\boldsymbol{\{X\}}
\end{equation}

\begin{align}
&\boldsymbol{\{x\}}=\{x,y,z\}^T  & &\text{: Lokale Koordinate}\\
&\boldsymbol{\{X\}}＝\{X,Y,Z\}^T & &\text{: Globale Koordinate}
\end{align}

\begin{equation}
\boldsymbol{[t]}=
\begin{bmatrix}
1 & 0 & 0 \\
0 & \cos\theta & \sin\theta \\
0 & -\sin\theta & \cos\theta
\end{bmatrix}
\cdot
\begin{bmatrix}
l & m & n \\
-\cfrac{m}{\sqrt{l^2+m^2}} & \cfrac{l}{\sqrt{l^2+m^2}} & 0 \\
-\cfrac{l \cdot n}{\sqrt{l^2+m^2}} & -\cfrac{m \cdot n}{\sqrt{l^2+m^2}} & \sqrt{l^2+m^2}
\end{bmatrix}
\qquad \text{($\theta$ : chord angle)}
\end{equation}

\begin{equation}
\boldsymbol{[t]}=
\begin{bmatrix}
1 & 0 & 0 \\
0 & \cos\theta & \sin\theta \\
0 & -\sin\theta & \cos\theta
\end{bmatrix}
\cdot
\begin{bmatrix}
0 & 0 & n \\
n & 0 & 0 \\
0 & 1 & 0
\end{bmatrix}
\qquad \text{($\theta$ : chord angle)}
\end{equation}

Koordinatenumrechnungsmatrix

\begin{equation}
\boldsymbol{[T]}=
\begin{bmatrix}
\boldsymbol{[t]} & \boldsymbol{0}   & \boldsymbol{0}   & \boldsymbol{0} \\\
\boldsymbol{0}   & \boldsymbol{[t]} & \boldsymbol{0}   & \boldsymbol{0} \\\
\boldsymbol{0}   & \boldsymbol{0}   & \boldsymbol{[t]} & \boldsymbol{0} \\\
\boldsymbol{0}   & \boldsymbol{0}   & \boldsymbol{0}   &\boldsymbol{[t]}
\end{bmatrix}
\end{equation}

Beziehung zwischen dem Elementkoordinatensystem (x-y-z) und dem Gesamtkoordinatensystem (X-Y-Z), wobei der Codewinkel 0 ist

Gesamtsteifigkeitsgleichung

Aus dem Obigen wird die Steifigkeitsgleichung im Gesamtkoordinatensystem wie folgt ausgedrückt.

\begin{equation}
[\boldsymbol{K}]\{\boldsymbol{U}\}=\{\boldsymbol{F}\}+\{\boldsymbol{F_b}\}+\{\boldsymbol{F_t}\}
\end{equation}

\begin{align}
&\{\boldsymbol{U}\}                          & & \text{Knotenverschiebungsvektor im globalen Koordinatensystem}\\
&\{\boldsymbol{F}\}                          & & \text{Knoten externer Kraftvektor im globalen Koordinatensystem}\\
&\{\boldsymbol{F_b}\}=\{\boldsymbol{f_b}\}   & & \text{Knotenträgheitskraftvektor im globalen Koordinatensystem}\\
&\{\boldsymbol{F_t}\}=[\boldsymbol{T}]^T\{\boldsymbol{f_t}\}  & & \text{Knotentemperaturlastvektor im globalen Koordinatensystem}\\
&[\boldsymbol{K}]=[\boldsymbol{T}]^T[\boldsymbol{k}][\boldsymbol{T}]   & & \text{Starrheitsmatrix im globalen Koordinatensystem}\\
&[\boldsymbol{k}]                            & & \text{Starrheitsmatrix im Elementkoordinatensystem}\\
&\{\boldsymbol{f_b}\}                          & & \text{Knotenträgheitskraftvektor im Elementkoordinatensystem}\\
&\{\boldsymbol{f_t}\}                          & & \text{Knotentemperaturlastvektor im Elementkoordinatensystem}\\
&[\boldsymbol{T}]                            & & \text{Koordinatenumrechnungsmatrix}
\end{align}

Durch Lösen der obigen Gleichung kann die Knotenverschiebung $ \ {\ boldsymbol {U} } $ im globalen Koordinatensystem erhalten werden.

Elementquerschnittskraft

Für die Elementquerschnittskraft wird die oben erhaltene Gesamtverschiebung des Koordinatensystems $ \ {\ boldsymbol {U} } $ durch die Koordinatenumwandlungsmatrix $ [\ boldsymbol {T}] $ in die Verschiebung des Elementkoordinatensystems umgewandelt. Berechnet durch Multiplikation der Elementsteifigkeitsmatrix $ [\ boldsymbol {k}] $. Das heißt, die Elementquerschnittskraft $ \ {\ boldsymbol {f_ {sec}} } $ wird wie folgt berechnet.

\begin{equation}
\{\boldsymbol{f_{sec}}\}=[\boldsymbol{k}]\cdot ([\boldsymbol{T}]\{\boldsymbol{U}\})
\end{equation}

Beachten Sie, dass die oben berechnete berechnete Axialkraft den Effekt der Temperaturänderung nicht berücksichtigt. Daher muss der Temperatureffekt der Axialkraft wie folgt korrigiert werden.

\begin{align}
N_1'=&N_1+EA\cdot\alpha\cdot\Delta T \\
N_2'=&N_2-EA\cdot\alpha\cdot\Delta T
\end{align}

Hier sind $ N_1 '$ und $ N_2' $ die Axialkräfte an den Knoten 1 und 2 unter Berücksichtigung der Temperaturänderung, und $ N_1 $ und $ N_2 $ sind die Axialkräfte, die aus der Verschiebung als Lösung der simultanen Gleichungen berechnet werden. $ EA $ ist die axiale Steifigkeit des Elements, $ \ alpha $ ist der lineare Ausdehnungskoeffizient des Elements und $ \ Delta T $ ist der Betrag der Temperaturänderung des Elements.

Programm zur Analyse dreidimensionaler Skelettstrukturen

Programmübersicht

Erklärung der Randbedingungen, die behandelt werden können

Element Beschreibung
Externe Knotenkraft Geben Sie die Anzahl der geladenen Knoten, die Anzahl der geladenen Knoten und den Lastwert an
Knotentemperaturänderung Geben Sie das Ausmaß der Temperaturänderung an allen Knoten an. Der Temperaturanstieg sei der positive Wert
Elementträgheitskraft Geben Sie in den Materialeigenschaftsdaten die Beschleunigung in x / y / z-Richtung als Verhältnis zur Schwerkraftbeschleunigung an
Knoten erzwungene Verschiebung Geben Sie die Anzahl der Knoten, die Knotennummer und den Verschiebungswert an, um die erzwungene Verschiebung zu erhalten

Programmkonfiguration (Programmname: py_fem_3dfrmx.py)

   1.Modul laden
   2.Funktion: Dateneingabe(def inpdata_3dfrm)
       (1)Grundlegende numerische Lesung
       (2)Array-Deklaration
       (3)Materialeigenschaftsdaten lesen
       (4)Lesen Sie die Knotennummer und die Materialnummer der Elementzusammensetzung
       (5)Lesen Sie die Knotenkoordinaten und die Änderung der Knotentemperatur
       (6)Ablesen der Randbedingungen
       (7)Knotenladedaten lesen
   3.Funktion: Eingabedaten exportieren(def prinp_3dfrm)
   4.Funktion: Berechnungsergebnis exportieren(def prout_3dfrm)
   5.Funktion: Erstellen Sie eine Elementsteifigkeitsmatrix(def sm_3dfrm)
   6.Funktion: Koordinatenkonvertierungsmatrix erstellen(def tm_3dfrm)
   7.Funktion: Erzeugung des Knotenträgheitskraftvektors(def bfvec_3dfrm)
   8.Funktion: Erstellung eines Knotentemperatur-Lastvektors(def tfvec_3dfrm)
   9.Funktion: Berechnung der Elementquerschnittskraft(def calsecf_3dfrm)
  10.Funktion: Hauptroutine(def main_3dfrm)
       (1)Rufen Sie den Namen der Eingabe- / Ausgabedatei über die Befehlszeile ab
       (2)Dateneingabe
       (3)Eingabedaten exportieren
       (4)Starre Matrixanordnung
       (5)Randbedingungsverarbeitung (Verarbeitung des bekannten Verschiebungsgrades)
       (6)Lösen Sie die Steifigkeitsgleichung (Verschiebungsberechnung)
       (7)Randbedingungsverarbeitung (Einbeziehung bekannter Verschiebungen)
       (8)Berechnung der Elementquerschnittskraft
       (9)Berechnungsergebnisse exportieren
       (10)Informationsanzeige wie Berechnungszeit
  11.Hauptroutinenausführung

Ausführungsbefehl und Eingabe- / Ausgabedaten

Analyse-Ausführungsbefehl

python3 py_fem_3dfrmx.py inp.txt out.txt
Element Beschreibung
py_fem_3dfrmx.py Python-FEM-Programmname
inp.txt Name der Eingabedatendatei (leere getrennte Daten)
out.txt Name der Ausgabedatendatei (leere getrennte Daten)

Eingabedatendateiformat

npoin  nele  nsec  npfix  nlod                   # Basic values for analysis
E  po  A  Ix  Iy  Iz  theta  alpha  gamma  gkX  gkY  gkZ
    ..... (1 to nsec) .....                      # Material properties
node_1  node_2  isec                             # Element connectivity, material set number
    ..... (1 to nele) .....                      #
x  y  z  deltaT                                  # Node coordinate, temperature change of node
    ..... (1 to npoin) .....                     #
lp  kox  koy  koz  kmx  kmy  kmz  rdis_x  rdis_y  rdis_z  rrot_x  rrot_y  rrot_z
    ..... (1 to npfix) .....                     # Restricted node and restrict conditions
lp  fp_x  fp_y  fp_z  mp_x  mp_y  mp_z           # Loaded node and loading conditions
    ..... (1 to nlod) .....                      #     (omit data input if nlod=0)
Element Beschreibung
npoin, nele, nsec Anzahl der Knoten, Anzahl der Elemente, Anzahl der Querschnittsmerkmale
npfix, nlod Anzahl der Rückhalteknoten, Anzahl der Ladeknoten
E, po, A Elementelastizitätskoeffizient, Element-Poisson-Verhältnis, Elementquerschnittsfläche
Ix, Iy, Iz Verdrehungskonstante, sekundäres Querschnittsmoment um die y-Achse, sekundäres Querschnittsmoment um die z-Achse >
Theta, Alpha Codewinkel, Elementlinien-Expansionskoeffizient
gamma, gkX, gkY, gkZ Volumenvolumeneinheit, Beschleunigung der x / y / z-Richtung (Verhältnis von g)
Knoten_1, Knoten_2, isec Knoten 1, Knoten 2, Schnittmerkmalnummer
x, y, z, deltaT Knoten x-Koordinate, Knoten y-Koordinate, Knoten z-Koordinate, Knotentemperaturänderung
lp, kox, koy, koz Beschränkungsknotennummer, x / y / z Richtungsverschiebung mit oder ohne Einschränkung (Einschränkung: 1, Freiheit: 0)
kmx, kmy, kmz x / y / z Mit oder ohne axiale Rotationsbeschränkung (Einschränkung: 1, Freiheit: 0) < / tr>
r_disx, r_disy, rdis_z x / y / z Richtung bekannte Verschiebung (geben Sie 0 ein, auch wenn keine Einschränkung vorliegt)
rrot_x, rrot_y, rrot_z x / y / z Bekannter Rotationsbetrag um die Achse (geben Sie 0 ein, auch wenn er nicht eingeschränkt ist)
lp, fp_x, fp_y, fp_z Knotennummer laden, Last in x-Richtung, Last in y-Richtung, Last in z-Richtung >
mp_x, mp_y, mp_z x-Achsenmoment (Verdrehung), y-Achsenmoment, z-Achsenmoment

Ausgabedatendateiformat

npoin  nele  nsec npfix  nlod
    (Each value of above)
sec  E      po     A    J    Iy    Iz    theta
sec  alpha  gamma  gkX  gkY  gkZ
    sec   : Material number
    E     : Elastic modulus
    po    : Poisson's ratio
    A     : Section area
    J     ; Torsional constant
    Iy    : Moment of inertia around y-axis
    Iz    : Moment of inertia around z-axis
    theta : Chord angle
    alpha : Thermal expansion coefficient
    gamma : Unit weight
    gkX   : Acceleration in X-direction (ratio to g)
    gkY   : Acceleration in Y-direction (ratio to g)
    gkZ   : Acceleration in Z-direction (ratio to g)
    ..... (1 to nsec) .....
node   x   y   z   fx   fy   fz   mx   my   mz  deltaT
    node   : Node number
    x      : x-coordinate
    y      : y-coordinate
    fx     : Load in X-direction
    fy     : Load in Y-direction
    fz     : Load in Z-direction
    mx     : Moment load around X-axis
    my     : Moment load around Y-axis
    mz     : Moment load around Z-axis
    deltaT : Temperature change of node
    ..... (1 to npoin) .....
node   kox   koy   koz   kmx   kmy   kmz   rdis_x   rdis_y   rdis_z   rrot_x   rrot_y   rrot_z
    node    : Node number
    kox     : Index of restriction in X-direction (0: free, 1: fixed)
    koy     : Index of restriction in Y-direction (0: free, 1: fixed)
    koz     : Index of restriction in Z-direction (0: free, 1: fixed)
    kmx     : Index of restriction in rotation around X-axis  (0: free, 1: fixed)
    kmy     : Index of restriction in rotation around Y-axis  (0: free, 1: fixed)
    kmz     : Index of restriction in rotation around Z-axis  (0: free, 1: fixed)
    rdis_x  : Given displacement in X-direction
    rdis_y  : Given displacement in Y-direction
    rdis_z  : Given displacement in Z-direction
    rrot_x  : Given displacement in rotation around X-axis
    rrot_y  : Given displacement in rotation around Y-axis
    rrot_z  : Given displacement in rotation around Z-axis
    ..... (1 to npfix) .....
elem     i     j   sec
    elem : Element number
    i    : Node number of start point
    j    : Node number of end point
    sec  : Material number
    ..... (1 to nele) .....
node   dis-x   dis-y   dis-z   rot-x   rot-y   rot-z
    node  : Node number
    dis-x : Displacement in X-direction
    dis-y : Displacement in Y-direction
    dis-z : Displacement in Z-direction
    rot-x : Rotation around X-axis
    rot-y : Rotation around Y-axis
    rot-z : Rotation around Z-axis
    ..... (1 to npoin) .....
elem nodei   N_i   Sy_i   Sz_i   Mx_i   My_i   Mz_i
elem nodej   N_j   Sy_j   Sz_j   Mx_j   My_j   Mz_j
    elem  : Element number
    nodei : First node nymber
    nodej : Second node number
    N_i   : Axial force of node-i
    Sy_i  : Shear force of node-i in y-direction
    Sz_i  : Shear force of node-i in z-direction
    Mx_i  : Moment of node-i around x-axis (torsional moment)
    My_i  : Moment of node-i around y-axis (bending moment)
    Mz_i  : Moment of node-i around z-axis (bending moment)
    N_j   : Axial force of node-j
    Sy_j  : Shear force of node-j in y-direction
    Sz_j  : Shear force of node-j in z-direction
    Mx_j  : Moment of node-j around x-axis (torsional moment)
    My_j  : Moment of node-j around y-axis (bending moment)
    Mz_j  : Moment of node-j around z-axis (bending moment)
    ..... (1 to nele) .....
n=(total degrees of freedom)  time=(calculation time)
Element Beschreibung
dis-x, dis-y, dis-z Verschiebung in x-Richtung, Verschiebung in y-Richtung, Verschiebung in z-Richtung
rot-x, rot-y, rot-z Rotationsbetrag der x-Achse, Rotationsbetrag der y-Achse, Rotationsbetrag der z-Achse td>
N, Sy, Sz Axialkraft, Scherkraft in y-Richtung, Scherkraft in z-Richtung
Mx, My, Mz x-Achsenmoment (Verdrehung), y-Achsen-Biegemoment, z-Achsen-Biegemoment

Programmierung

Der Hauptteil des Programms wird unten zusammen mit der "Programmstruktur" gezeigt.

1. Laden Sie das Modul

# ======================
# 3D Frame Analysis
# ======================
import numpy as np
from scipy.sparse.linalg import spsolve
from scipy.sparse import csr_matrix
import sys
import time

2. Funktion: Dateneingabe (def inpdata_3dfrm)

def inpdata_3dfrm(fnameR,nod,nfree):
    f=open(fnameR,'r')
    text=f.readline()
    text=text.strip()
    text=text.split()
    npoin=int(text[0]) # Number of nodes
    nele =int(text[1]) # Number of elements
    nsec =int(text[2]) # Number of sections
    npfix=int(text[3]) # Number of restricted nodes
    nlod =int(text[4]) # Number of loaded nodes
    # array declaration
    ae    =np.zeros((12,nsec),dtype=np.float64)      # Section characteristics
    node  =np.zeros((nod+1,nele),dtype=np.int)       # node-element relationship
    x     =np.zeros((3,npoin),dtype=np.float64)      # Coordinates of nodes
    deltaT=np.zeros(npoin,dtype=np.float64)          # Temperature change of nodes
    mpfix =np.zeros((nfree,npoin),dtype=np.int)      # Ristrict conditions
    rdis  =np.zeros((nfree,npoin),dtype=np.float64)  # Ristricted displacement
    fp    =np.zeros(nfree*npoin,dtype=np.float64)    # External force vector
    # section characteristics
    for i in range(0,nsec):
        text=f.readline()
        text=text.strip()
        text=text.split()
        ae[0,i] =float(text[0])  # E     : elastic modulus
        ae[1,i] =float(text[1])  # po    : Poisson's ratio
        ae[2,i] =float(text[2])  # aa    : section area
        ae[3,i] =float(text[3])  # aix   : tortional constant
        ae[4,i] =float(text[4])  # aiy   : moment of inertia aroutd y-axis
        ae[5,i] =float(text[5])  # aiz   : moment of inertia around z-axis
        ae[6,i] =float(text[6])  # theta : chord angle
        ae[7,i] =float(text[7])  # alpha : thermal expansion coefficient
        ae[8,i] =float(text[8])  # gamma : unit weight of material
        ae[9,i] =float(text[9])  # gkX   : acceleration in X-direction
        ae[10,i]=float(text[10]) # gkY   : acceleration in Y-direction
        ae[11,i]=float(text[11]) # gkZ   : acceleration in Z-direction
    # element-node
    for i in range(0,nele):
        text=f.readline()
        text=text.strip()
        text=text.split()
        node[0,i]=int(text[0]) #node_1
        node[1,i]=int(text[1]) #node_2
        node[2,i]=int(text[2]) #section characteristic number
    # node coordinates
    for i in range(0,npoin):
        text=f.readline()
        text=text.strip()
        text=text.split()
        x[0,i]=float(text[0])    # x-coordinate
        x[1,i]=float(text[1])    # y-coordinate
        x[2,i]=float(text[2])    # z-coordinate
        deltaT[i]=float(text[3]) # Temperature change
    # boundary conditions (0:free, 1:restricted)
    for i in range(0,npfix):
        text=f.readline()
        text=text.strip()
        text=text.split()
        lp=int(text[0])              #fixed node
        mpfix[0,lp-1]=int(text[1])   #fixed in x-direction
        mpfix[1,lp-1]=int(text[2])   #fixed in y-direction
        mpfix[2,lp-1]=int(text[3])   #fixed in z-direction
        mpfix[3,lp-1]=int(text[4])   #fixed in rotation around x-axis
        mpfix[4,lp-1]=int(text[5])   #fixed in rotation around y-axis
        mpfix[5,lp-1]=int(text[6])   #fixed in rotation around z-axis
        rdis[0,lp-1]=float(text[7])  #fixed displacement in x-direction
        rdis[1,lp-1]=float(text[8])  #fixed displacement in y-direction
        rdis[2,lp-1]=float(text[9])  #fixed displacement in z-direction
        rdis[3,lp-1]=float(text[10]) #fixed rotation around x-axis
        rdis[4,lp-1]=float(text[11]) #fixed rotation around y-axis
        rdis[5,lp-1]=float(text[12]) #fixed rotation around z-axis
    # load
    if 0&lt;nlod:
        for i in range(0,nlod):
            text=f.readline()
            text=text.strip()
            text=text.split()
            lp=int(text[0])           #loaded node
            fp[6*lp-6]=float(text[1]) #load in x-direction
            fp[6*lp-5]=float(text[2]) #load in y-direction
            fp[6*lp-4]=float(text[3]) #load in z-direction
            fp[6*lp-3]=float(text[4]) #moment around x-axis
            fp[6*lp-2]=float(text[5]) #moment around y-axis
            fp[6*lp-1]=float(text[6]) #moment around z-axis
    f.close()
    return npoin,nele,nsec,npfix,nlod,ae,node,x,deltaT,mpfix,rdis,fp

3. Funktion: Eingabedaten schreiben (def prinp_3dfrm)

def prinp_3dfrm(fnameW,npoin,nele,nsec,npfix,nlod,ae,node,x,deltaT,mpfix,rdis,fp):
    fout=open(fnameW,'w')
    # print out of input data
    print('{0:&gt;5s} {1:&gt;5s} {2:&gt;5s} {3:&gt;5s} {4:&gt;5s}'.format('npoin','nele','nsec','npfix','nlod'),file=fout)
    print('{0:5d} {1:5d} {2:5d} {3:5d} {4:5d}'.format(npoin,nele,nsec,npfix,nlod),file=fout)
    print('{0:&gt;5s} {1:&gt;15s} {2:&gt;15s} {3:&gt;15s} {4:&gt;15s} {5:&gt;15s} {6:&gt;15s} {7:&gt;15s}'
    .format('sec','E','po','A','J','Iy','Iz','theta'),file=fout)
    print('{0:&gt;5s} {1:&gt;15s} {2:&gt;15s} {3:&gt;15s} {4:&gt;15s} {5:&gt;15s}'
    .format('sec','alpha','gamma','gkX','gkY','gkZ'),file=fout)
    for i in range(0,nsec):
        print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:15.7e} {7:15.7e}'
        .format(i+1,ae[0,i],ae[1,i],ae[2,i],ae[3,i],ae[4,i],ae[5,i],ae[6,i]),file=fout)
        print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e}'
        .format(i+1,ae[7,i],ae[8,i],ae[9,i],ae[10,i],ae[11,i]),file=fout)
    print('{0:&gt;5s} {1:&gt;15s} {2:&gt;15s} {3:&gt;15s} {4:&gt;15s} {5:&gt;15s} {6:&gt;15s} {7:&gt;15s} {8:&gt;15s} {9:&gt;15s} {10:&gt;15s}'
    .format('node','x','y','z','fx','fy','fz','mx','my','mz','deltaT'),file=fout)
    for i in range(0,npoin):
        lp=i+1
        print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:15.7e} {7:15.7e} {8:15.7e} {9:15.7e} {10:15.7e}'
        .format(lp,x[0,i],x[1,i],x[2,i],fp[6*lp-6],fp[6*lp-5],fp[6*lp-4],fp[6*lp-3],fp[6*lp-2],fp[6*lp-1],deltaT[i]),file=fout)
    print('{0:&gt;5s} {1:&gt;5s} {2:&gt;5s} {3:&gt;5s} {4:&gt;5s} {5:&gt;5s} {6:&gt;5s} {7:&gt;15s} {8:&gt;15s} {9:&gt;15s} {10:&gt;15s} {11:&gt;15s} {12:&gt;15s}'
    .format('node','kox','koy','koz','kmx','kmy','kmz','rdis_x','rdis_y','rdis_z','rrot_x','rrot_y','rrot_z'),file=fout)
    for i in range(0,npoin):
        if 0&lt;mpfix[0,i]+mpfix[1,i]+mpfix[2,i]+mpfix[3,i]+mpfix[4,i]+mpfix[5,i]:
            lp=i+1
            print('{0:5d} {1:5d} {2:5d} {3:5d} {4:5d} {5:5d} {6:5d} {7:15.7e} {8:15.7e} {9:15.7e} {10:15.7e} {11:15.7e} {12:15.7e}'
            .format(lp,mpfix[0,i],mpfix[1,i],mpfix[2,i],mpfix[3,i],mpfix[4,i],mpfix[5,i],rdis[0,i],rdis[1,i],rdis[2,i],rdis[3,i],rdis[4,i],rdis[5,i]),file=fout)
    print('{0:&gt;5s} {1:&gt;5s} {2:&gt;5s} {3:&gt;5s}'.format('elem','i','j','sec'),file=fout)
    for ne in range(0,nele):
        print('{0:5d} {1:5d} {2:5d} {3:5d}'.format(ne+1,node[0,ne],node[1,ne],node[2,ne]),file=fout)
    fout.close()

4. Funktion: Berechnungsergebnis exportieren (def prout_3dfrm)

def prout_3dfrm(fnameW,npoin,nele,node,disg,fsec):
    fout=open(fnameW,'a')
    # displacement
    print('{0:&gt;5s} {1:&gt;15s} {2:&gt;15s} {3:&gt;15s} {4:&gt;15s} {5:&gt;15s} {6:&gt;15s}'
          .format('node','dis-x','dis-y','dis-z','rot-x','rot-y','rot-z'),file=fout)
    for i in range(0,npoin):
        lp=i+1
        print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:15.7e}'
              .format(lp,disg[6*lp-6],disg[6*lp-5],disg[6*lp-4],disg[6*lp-3],disg[6*lp-2],disg[6*lp-1]),file=fout)
    # section force
    print('{0:&gt;5s} {1:&gt;5s} {2:&gt;15s} {3:&gt;15s} {4:&gt;15s} {5:&gt;15s} {6:&gt;15s} {7:&gt;15s}'
    .format('elem','nodei','N_i','Sy_i','Sz_i','Mx_i','My_i','Mz_i'),file=fout)
    print('{0:&gt;5s} {1:&gt;5s} {2:&gt;15s} {3:&gt;15s} {4:&gt;15s} {5:&gt;15s} {6:&gt;15s} {7:&gt;15s}'
    .format('elem','nodej','N_j','Sy_j','Sz_j','Mx_j','My_j','Mz_j'),file=fout)
    for ne in range(0,nele):
        print('{0:5d} {1:5d} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:15.7e} {7:15.7e}'
        .format(ne+1,node[0,ne],fsec[0,ne],fsec[1,ne],fsec[2,ne],fsec[3,ne],fsec[4,ne],fsec[5,ne]),file=fout)
        print('{0:5d} {1:5d} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:15.7e} {7:15.7e}'
        .format(ne+1,node[1,ne],fsec[6,ne],fsec[7,ne],fsec[8,ne],fsec[9,ne],fsec[10,ne],fsec[11,ne]),file=fout)
    fout.close()

5. Funktion: Elementsteifigkeitsmatrix erstellen (def sm_3dfrm)

\begin{equation}
\boldsymbol{[k]}=
\begin{bmatrix}
\cfrac{EA}{L} & 0 & 0 & 0 & 0 & 0 & -\cfrac{EA}{L} & 0 & 0 & 0 & 0 & 0 \\
0 & \cfrac{12EI_z}{L^3} & 0 & 0 & 0 & \cfrac{6EI_z}{L^2} & 0 & -\cfrac{12EI_z}{L^3} & 0 & 0 & 0 & \cfrac{6EI_z}{L^2} \\
0 & 0 & \cfrac{12EI_y}{L^3} & 0 & -\cfrac{6EI_y}{L^2} & 0 & 0 & 0 & -\cfrac{12EI_y}{L^3} & 0 & -\cfrac{6EI_y}{L^2} & 0 \\
0 & 0 & 0 & \cfrac{GJ}{L} & 0 & 0 & 0 & 0 & 0 & -\cfrac{GJ}{L} & 0 & 0 \\
0 & 0 & -\cfrac{6EI_y}{L^2} & 0 & \cfrac{4EI_y}{L} & 0 & 0 & 0 & \cfrac{6EI_y}{L^2} & 0 & \cfrac{2EI_y}{L} & 0 \\
0 & \cfrac{6EI_z}{L^2} & 0 & 0 & 0 & \cfrac{4EI_z}{L} & 0 & -\cfrac{6EI_z}{L^2} & 0 & 0 & 0 & \cfrac{2EI_z}{L} \\
-\cfrac{EA}{L} & 0 & 0 & 0 & 0 & 0 & \cfrac{EA}{L} & 0 & 0 & 0 & 0 & 0 \\
0 & -\cfrac{12EI_z}{L^3} & 0 & 0 & 0 & -\cfrac{6EI_z}{L^2} & 0 & \cfrac{12EI_z}{L^3} & 0 & 0 & 0 & -\cfrac{6EI_z}{L^2} \\
0 & 0 & -\cfrac{12EI_y}{L^3} & 0 & \cfrac{6EI_y}{L^2} & 0 & 0 & 0 & \cfrac{12EI_y}{L^3} & 0 & \cfrac{6EI_y}{L^2} & 0 \\
0 & 0 & 0 & -\cfrac{GJ}{L} & 0 & 0 & 0 & 0 & 0 & \cfrac{GJ}{L} & 0 & 0 \\
0 & 0 & -\cfrac{6EI_y}{L^2} & 0 & \cfrac{2EI_y}{L} & 0 & 0 & 0 & \cfrac{6EI_y}{L^2} & 0 & \cfrac{4EI_y}{L} & 0 \\
0 & \cfrac{6EI_z}{L^2} & 0 & 0 & 0 & \cfrac{2EI_z}{L} & 0 & -\cfrac{6EI_z}{L^2} & 0 & 0 & 0 & \cfrac{4EI_z}{L} \\
\end{bmatrix}
\end{equation}

def sm_3dfrm(EA,GJ,EIy,EIz,x1,y1,z1,x2,y2,z2):
    ek=np.zeros((12,12),dtype=np.float64) # local stiffness matrix
    xx=x2-x1
    yy=y2-y1
    zz=z2-z1
    el=np.sqrt(xx**2+yy**2+zz**2)
    ek[ 0, 0]= EA/el
    ek[ 0, 6]=-EA/el
    ek[ 1, 1]= 12*EIz/el**3
    ek[ 1, 5]=  6*EIz/el**2
    ek[ 1, 7]=-12*EIz/el**3
    ek[ 1,11]=  6*EIz/el**2
    ek[ 2, 2]= 12*EIy/el**3
    ek[ 2, 4]= -6*EIy/el**2
    ek[ 2, 8]=-12*EIy/el**3
    ek[ 2,10]= -6*EIy/el**2
    ek[ 3, 3]= GJ/el
    ek[ 3, 9]=-GJ/el
    ek[ 4, 2]= -6*EIy/el**2
    ek[ 4, 4]=  4*EIy/el
    ek[ 4, 8]=  6*EIy/el**2
    ek[ 4,10]=  2*EIy/el
    ek[ 5, 1]=  6*EIz/el**2
    ek[ 5, 5]=  4*EIz/el
    ek[ 5, 7]= -6*EIz/el**2
    ek[ 5,11]=  2*EIz/el
    ek[ 6, 0]=-EA/el
    ek[ 6, 6]= EA/el
    ek[ 7, 1]=-12*EIz/el**3
    ek[ 7, 5]= -6*EIz/el**2
    ek[ 7, 7]= 12*EIz/el**3
    ek[ 7,11]= -6*EIz/el**2
    ek[ 8, 2]=-12*EIy/el**3
    ek[ 8, 4]=  6*EIy/el**2
    ek[ 8, 8]= 12*EIy/el**3
    ek[ 8,10]=  6*EIy/el**2
    ek[ 9, 3]=-GJ/el
    ek[ 9, 9]= GJ/el
    ek[10, 2]= -6*EIy/el**2
    ek[10, 4]=  2*EIy/el
    ek[10, 8]=  6*EIy/el**2
    ek[10,10]=  4*EIy/el
    ek[11, 1]=  6*EIz/el**2
    ek[11, 5]=  2*EIz/el
    ek[11, 7]= -6*EIz/el**2
    ek[11,11]=  4*EIz/el
    return ek

6. Funktion: Koordinatentransformationsmatrix erstellen (def tm_3dfrm)

\begin{gather}
\ell=\sqrt{(x_j-x_i)^2+(y_j-y_i)^2+(z_j-z_i)^2}\\
l=\cfrac{x_j-x_i}{\ell} \qquad m=\cfrac{y_j-y_i}{\ell} \qquad n=\cfrac{z_j-z_i}{\ell}
\end{gather}

\begin{equation}
\boldsymbol{[t_1]}=
\begin{bmatrix}
1 & 0 & 0 \\
0 & \cos\theta & \sin\theta \\
0 & -\sin\theta & \cos\theta
\end{bmatrix}
\qquad \text{($\theta$ : chord angle)}
\end{equation}

\begin{align}
&\boldsymbol{[t_2]}=
\begin{bmatrix}
0 & 0 & n \\
n & 0 & 0 \\
0 & 1 & 0
\end{bmatrix}
& &\text{Wenn die Elementlängsachse (Elementkoordinatensystem x-Achse) parallel zur Gesamtkoordinatensystem-Z-Achse ist}\\
&\boldsymbol{[t_2]}=
\begin{bmatrix}
l & m & n \\
-\cfrac{m}{\sqrt{l^2+m^2}} & \cfrac{l}{\sqrt{l^2+m^2}} & 0 \\
-\cfrac{l \cdot n}{\sqrt{l^2+m^2}} & -\cfrac{m \cdot n}{\sqrt{l^2+m^2}} & \sqrt{l^2+m^2}
\end{bmatrix}
& &\text{Wenn die Elementlängsachse (Elementkoordinatensystem x-Achse) nicht parallel zur Gesamtkoordinatensystem-Z-Achse ist}
\end{align}

\begin{equation}
\boldsymbol{[t]}=\boldsymbol{[t_1]}\cdot\boldsymbol{[t_2]}
\end{equation}

\begin{equation}
\boldsymbol{[T]}=
\begin{bmatrix}
\boldsymbol{[t]} & \boldsymbol{0}   & \boldsymbol{0}   & \boldsymbol{0} \\
\boldsymbol{0}   & \boldsymbol{[t]} & \boldsymbol{0}   & \boldsymbol{0} \\
\boldsymbol{0}   & \boldsymbol{0}   & \boldsymbol{[t]} & \boldsymbol{0} \\
\boldsymbol{0}   & \boldsymbol{0}   & \boldsymbol{0}   &\boldsymbol{[t]}
\end{bmatrix}
\end{equation}

def tm_3dfrm(theta,x1,y1,z1,x2,y2,z2):
    tt=np.zeros((12,12),dtype=np.float64) # transformation matrix
    t1=np.zeros((3,3),dtype=np.float64)
    t2=np.zeros((3,3),dtype=np.float64)
    xx=x2-x1
    yy=y2-y1
    zz=z2-z1
    el=np.sqrt(xx**2+yy**2+zz**2)
    t1[0,0]=1
    t1[1,1]= np.cos(theta)
    t1[1,2]= np.sin(theta)
    t1[2,1]=-np.sin(theta)
    t1[2,2]= np.cos(theta)
    ll=(x2-x1)/el
    mm=(y2-y1)/el
    nn=(z2-z1)/el
    if x2-x1==0.0 and y2-y1==0.0:
        t2[0,2]=nn
        t2[1,0]=nn
        t2[2,1]=1.0
    else:
        qq=np.sqrt(ll**2+mm**2)
        t2[0,0]=ll
        t2[0,1]=mm
        t2[0,2]=nn
        t2[1,0]=-mm/qq
        t2[1,1]= ll/qq
        t2[2,0]=-ll*nn/qq
        t2[2,1]=-mm*nn/qq
        t2[2,2]=qq
    t3=np.dot(t1,t2)
    tt[0:3,0:3]  =t3[0:3,0:3]
    tt[3:6,3:6]  =t3[0:3,0:3]
    tt[6:9,6:9]  =t3[0:3,0:3]
    tt[9:12,9:12]=t3[0:3,0:3]
    return tt

7. Funktion: Erzeugung eines Knotenträgheitskraftvektors (def bfvec_3dfrm)

\begin{equation*}
\{\boldsymbol{F_{b}}\}=\{\boldsymbol{f_{b}}\}=\frac{\gamma\cdot A \cdot L}{2}
\begin{Bmatrix}
k_x \\
k_y \\
k_z \\
0   \\
0   \\
0   \\
k_x \\
k_y \\
k_z \\
0   \\
0   \\
0
\end{Bmatrix}
\end{equation*}

def bfvec_3dfrm(A,gamma,gkX,gkY,gkZ,x1,y1,z1,x2,y2,z2):
    # Body force vector in global coordinate system
    bfe_g=np.zeros(12,dtype=np.float64)
    xx=x2-x1
    yy=y2-y1
    zz=z2-z1
    el=np.sqrt(xx**2+yy**2+zz**2)
    bfe_g[0]=0.5*gamma*A*el*gkX
    bfe_g[1]=0.5*gamma*A*el*gkY
    bfe_g[2]=0.5*gamma*A*el*gkZ
    bfe_g[6]=0.5*gamma*A*el*gkX
    bfe_g[7]=0.5*gamma*A*el*gkY
    bfe_g[8]=0.5*gamma*A*el*gkZ
    return bfe_g

8. Funktion: Erstellung des Knotentemperatur-Lastvektors (def tfvec_3dfrm)

\begin{equation*}
\{\boldsymbol{f_{t}}\}=EA\cdot\alpha\cdot\Delta T
\begin{Bmatrix}
-1 \\
 0 \\
 0 \\
 0 \\
 0 \\
 0 \\
 1 \\
 0 \\
 0 \\
 0 \\
 0 \\
 0
\end{Bmatrix}
\end{equation*}

def tfvec_3dfrm(EA,alpha,tem):
    # Thermal load vector  in local coordinate system
    tfe_l=np.zeros(12,dtype=np.float64)
    tfe_l[0]=-EA*alpha*tem
    tfe_l[6]= EA*alpha*tem
    return tfe_l

9. Funktion: Berechnung der Elementquerschnittskraft (def calsecf_3dfrm)

$ N_1 '$ und $ N_2' $ sind Axialkraftkorrekturen aufgrund von Temperaturänderungen.

\begin{align}
&\{\boldsymbol{f_{sec}}\}=[\boldsymbol{k}]\cdot ([\boldsymbol{T}]\{\boldsymbol{U}\})\\
&N_1'=N_1+EA\cdot\alpha\cdot\Delta T \\
&N_2'=N_2-EA\cdot\alpha\cdot\Delta T
\end{align}

def calsecf_3dfrm(EA,GJ,EIy,EIz,theta,alpha,tem,dis,x1,y1,z1,x2,y2,z2):
    ek=sm_3dfrm(EA,GJ,EIy,EIz,x1,y1,z1,x2,y2,z2) # Stiffness matrix in local coordinate
    tt=tm_3dfrm(theta,x1,y1,z1,x2,y2,z2) # Transformation matrix
    secf=np.dot(ek,np.dot(tt,dis))
    secf[0]=secf[0]+EA*alpha*tem
    secf[6]=secf[6]-EA*alpha*tem
    return secf

10. Funktion: Hauptroutine

Hauptroutinenverarbeitungsverfahren

Element Beschreibung
npoin Gesamtzahl der Knoten
nfree 1 Knotenfreiheit (6 hier)
mpfix Array mit Einschränkungsbedingungen für jeden Knoten (= 1: Einschränkung, = 0: frei)
fp Knotenlastvektor
rdis Verschiebungsbetrag des Knotens, der die Verschiebung einschränkt
gk Gesamtsteifigkeitsmatrix
disg Lösung (Verschiebung) der Gesamtsteifigkeitsgleichung

  • Lösen der Steifigkeitsgleichung (Verschiebungsberechnung) </ b>: Nach der CSR-Konvertierung der Gesamtsteifigkeitsmatrix gk </ var> werden die simultanen linearen Gleichungen mit scipy.sparse.linalg.spsolve () </ var> gelöst.

  • Randbedingungsverarbeitung (Einbeziehung bekannter Verschiebungen) </ b>: Vergessen Sie nach dem Lösen der simultanen linearen Gleichungen nicht, die bekannte Verschiebung in dem der bekannten Verschiebung entsprechenden Freiheitsbegriff erneut einzugeben.

  • Berechnung der Elementquerschnittskraft </ b>: Führen Sie für jedes Element eine Querschnittskraftberechnung mit der Funktion calsecf_3dfrm </ em> durch

  • Berechnungsergebnisse exportieren </ b>: Berechnungsergebnis, das von der Funktion prout_3dfrm </ em> ausgegeben wird

  • Informationsanzeige wie Berechnungszeit </ b>: Die Differenz zwischen der Abschlusszeit der Berechnung und der Startzeit wird wie folgt als Verarbeitungszeit ausgegeben.

import time

start=time.time()
......
dtime=time.time()-start

Hauptroutineprogramm

def main_3dfrm():
    start=time.time()
    args = sys.argv
    fnameR=args[1] # input data file
    fnameW=args[2] # output data file
    nod=2              # Number of nodes per element
    nfree=6            # Degree of freedom per node
    # data input
    npoin,nele,nsec,npfix,nlod,ae,node,x,deltaT,mpfix,rdis,fp=inpdata_3dfrm(fnameR,nod,nfree)
    # print out of input data
    prinp_3dfrm(fnameW,npoin,nele,nsec,npfix,nlod,ae,node,x,deltaT,mpfix,rdis,fp)
    # array declaration
    ir=np.zeros(nod*nfree,dtype=np.int)   # Work vector for matrix assembly
    gk=np.zeros((nfree*npoin,nfree*npoin),dtype=np.float64)   # Global stiffness matrix
    # assembly of stiffness matrix & load vectors
    for ne in range(0,nele):
        i=node[0,ne]-1
        j=node[1,ne]-1
        m=node[2,ne]-1
        x1=x[0,i]; y1=x[1,i]; z1=x[2,i]
        x2=x[0,j]; y2=x[1,j]; z2=x[2,j]
        ee   =ae[0,m]  # elastic modulus
        po   =ae[1,m]  # Poisson's ratio
        aa   =ae[2,m]  # section area
        aix  =ae[3,m] # tortional constant
        aiy  =ae[4,m] # moment of inertia around y-axis
        aiz  =ae[5,m] # moment of inertia around z-axis
        theta=np.radians(ae[6,m]) # chord angle
        alpha=ae[7,m] # thermal expansion coefficient
        gamma=ae[8,m] # unit weight of material
        gkX  =ae[9,m]   # seismic coefficient in X-direction
        gkY  =ae[10,m]  # seismic coefficient in Y-direction
        gkZ  =ae[11,m]  # seismic coefficient in Z-direction
        A=aa  # section area
        EA=ee*aa
        GJ=ee/2/(1+po)*aix
        EIy=ee*aiy
        EIz=ee*aiz
        tem=0.5*(deltaT[i]+deltaT[j])                # average temperature change
        tt   =tm_3dfrm(theta,x1,y1,z1,x2,y2,z2)         # Transformation matrix
        ek   =sm_3dfrm(EA,GJ,EIy,EIz,x1,y1,z1,x2,y2,z2) # Stiffness matrix in local coordinate
        ck   =np.dot(np.dot(tt.T,ek),tt)             # Stiffness matrix in global coordinate
        tfe_l=tfvec_3dfrm(EA,alpha,tem)              # Thermal load vector in local coordinate
        tfe  =np.dot(tt.T,tfe_l)                     # Thermal load vector in global coordinate
        bfe  =bfvec_3dfrm(A,gamma,gkX,gkY,gkZ,x1,y1,z1,x2,y2,z2) # Body force vector in global coordinate
        ir[11]=6*j+5; ir[10]=ir[11]-1; ir[9]=ir[10]-1; ir[8]=ir[9]-1; ir[7]=ir[8]-1; ir[6]=ir[7]-1
        ir[5] =6*i+5; ir[4] =ir[5]-1 ; ir[3]=ir[4]-1 ; ir[2]=ir[3]-1; ir[1]=ir[2]-1; ir[0]=ir[1]-1
        for i in range(0,nod*nfree):
            it=ir[i]
            fp[it]=fp[it]+tfe[i]+bfe[i]
            for j in range(0,nod*nfree):
                jt=ir[j]
                gk[it,jt]=gk[it,jt]+ck[i,j]
    # treatment of boundary conditions
    for i in range(0,npoin):
        for j in range(0,nfree):
            if mpfix[j,i]==1:
                iz=i*nfree+j
                fp[iz]=0.0
    for i in range(0,npoin):
        for j in range(0,nfree):
            if mpfix[j,i]==1:
                iz=i*nfree+j
                fp=fp-rdis[j,i]*gk[:,iz]
                gk[:,iz]=0.0
                gk[iz,iz]=1.0
    # solution of simultaneous linear equations
    #disg = np.linalg.solve(gk, fp)
    gk = csr_matrix(gk)
    disg = spsolve(gk, fp, use_umfpack=True)
    # recovery of restricted displacements
    for i in range(0,npoin):
        for j in range(0,nfree):
            if mpfix[j,i]==1:
                iz=i*nfree+j
                disg[iz]=rdis[j,i]
    # calculation of section force
    dis=np.zeros(12,dtype=np.float64)
    fsec  =np.zeros((nod*nfree,nele),dtype=np.float64)  # Section force vector
    for ne in range(0,nele):
        i=node[0,ne]-1
        j=node[1,ne]-1
        m=node[2,ne]-1
        x1=x[0,i]; y1=x[1,i]; z1=x[2,i]
        x2=x[0,j]; y2=x[1,j]; z2=x[2,j]
        ee   =ae[0,m]  # elastic modulus
        po   =ae[1,m]  # Poisson's ratio
        aa   =ae[2,m]  # section area
        aix  =ae[3,m] # tortional constant
        aiy  =ae[4,m] # moment of inertia around y-axis
        aiz  =ae[5,m] # moment of inertia around z-axis
        theta=np.radians(ae[6,m]) # chord angle
        alpha=ae[7,m] # thermal expansion coefficient
        EA=ee*aa
        GJ=ee/2/(1+po)*aix
        EIy=ee*aiy
        EIz=ee*aiz
        tem=0.5*(deltaT[i]+deltaT[j]) # average temperature change
        dis[0]=disg[6*i]  ; dis[1] =disg[6*i+1]; dis[2]= disg[6*i+2]
        dis[3]=disg[6*i+3]; dis[4] =disg[6*i+4]; dis[5]= disg[6*i+5]
        dis[6]=disg[6*j]  ; dis[7] =disg[6*j+1]; dis[8]= disg[6*j+2]
        dis[9]=disg[6*j+3]; dis[10]=disg[6*j+4]; dis[11]=disg[6*j+5]
        fsec[:,ne]=calsecf_3dfrm(EA,GJ,EIy,EIz,theta,alpha,tem,dis,x1,y1,z1,x2,y2,z2)
    # print out of result
    prout_3dfrm(fnameW,npoin,nele,node,disg,fsec)
    # information
    dtime=time.time()-start
    print('n={0}  time={1:.3f}'.format(nfree*npoin,dtime)+' sec')
    fout=open(fnameW,'a')
    print('n={0}  time={1:.3f}'.format(nfree*npoin,dtime)+' sec',file=fout)
    fout.close()

11. Hauptroutinenausführung

#==============
# Execution
#==============
if __name__ == '__main__': main_3dfrm()

das ist alles

Recommended Posts