Analyse de la structure du squelette en trois dimensions avec Python

Gros correctif

Aperçu

J'ai créé un programme d'analyse de structure de squelette en trois dimensions que j'ai toujours voulu en Python. Dans le programme d'analyse de squelette de plan créé précédemment, la matrice de rigidité et la matrice de conversion de coordonnées sont principalement modifiées pour prendre en charge les éléments de poutre avec 1 nœud et 6 degrés de liberté.

~~ Le test n'est pas encore suffisant, mais ~~ J'ai été inspiré par l'article ci-dessous et je l'ai posté. http://qiita.com/sasaco/items/917429a497c6e9a08401

~~ Puisque le programme est créé sur Jupyter, cet article est également décrit selon la description de mon Jupyter. Le programme d'analyse de Python est divisé en trois parties. ~~ Je suis toujours inquiet, mais la partie exportation de fichiers est trop encombrée. Tu ne peux pas mieux écrire? .. ..

~~ Le calcul est en cours. .. .. ~~

L'environnement de programmation est Machine: MacBook Pro (Retina, 13-inch, Mid 2014) OS: macOS Sierra Python version: Python 3.6.1

Théorie de l'analyse de la structure du squelette en trois dimensions

Questions de base

Les conditions suivantes sont incluses dans le développement et le programme de la théorie.

Dans ce programme, en supposant que le modèle sera à grande échelle, un traitement matriciel épars est introduit lors de la résolution d'équations linéaires simultanées.

Équation de rigidité des éléments

La forme générale de l'équation de la rigidité de l'élément est la même que dans l'analyse bidimensionnelle des contraintes.

Equation pour trouver le déplacement

\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{(Matrice de rigidité des éléments)}\\
&\{\boldsymbol{f}\}=\int_A [\boldsymbol{N}]^T\{\boldsymbol{S}\}dA  & & \text{(Vecteur de force externe nodale)}\\
&\{\boldsymbol{f_t}\}=\int_V[\boldsymbol{B}]^T[\boldsymbol{D_e}]\{\boldsymbol{\epsilon_0}\}dV  & & \text{(Vecteur de distorsion du nœud initial)}\\
&\{\boldsymbol{f_b}\}=\gamma\left(\int_V [\boldsymbol{N}]^T [\boldsymbol{N}]dV\right)\{\boldsymbol{w_{nd}}\} & & \text{(Vecteur de force d'inertie nodale)}
\end{align}

Formule de calcul des contraintes internes de l'élément

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

\begin{align}
&\{\boldsymbol{\epsilon}\}  & &\text{: Vecteur de déformation de l'élément obtenu à partir du déplacement obtenu en résolvant l'équation de rigidité}\\
&\{\boldsymbol{\epsilon_0}\} & &\text{: Déformation initiale telle que déformation thermique}\\
&[\boldsymbol{D_e}]          & &\text{: Matrice de relation contrainte-déformation}
\end{align}

Matrice de rigidité des éléments de la structure tridimensionnelle

La matrice de rigidité des éléments de la structure de squelette tridimensionnelle est une matrice carrée de 12 x 12. Cela peut être dérivé en étendant celui du squelette plat, sauf pour considérer la torsion. Le processus de dérivation est omis et seuls les résultats sont décrits ici.

Équation de rigidité d'élément (expression la plus simple)

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

Vecteur de déplacement nodal

\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{: Déplacement sur l'axe X} & &\theta_x& &\text{: Quantité de rotation de l'axe X (angle de torsion)}\\
&v & &\text{: Déplacement sur l'axe Y} & &\theta_y& &\text{: Quantité de rotation de l'axe Y}\\
&w & &\text{: Déplacement sur l'axe Z} & &\theta_z& &\text{: Quantité de rotation de l'axe Z}
\end{align}

Vecteur de force nodale

\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{: Force axiale (force de direction de l'axe x)} & &M_x& &\text{: Moment de l'axe X (moment de torsion)}\\
&S_y& &\text{: Force de cisaillement sur l'axe Y}   & &M_y& &\text{: Moment de flexion sur l'axe Y}\\
&S_z& &\text{: Force de cisaillement sur l'axe Z}   & &M_z& &\text{: Moment de flexion sur l'axe Z}
\end{align}

Matrice de rigidité des éléments

\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{: Longueur de l'élément}        & &J& &\text{: Constante de torsion}\\
&E& &\text{: Coefficient d'élasticité}      & &I_y& &\text{: Moment secondaire de la section transversale de l'axe Y}\\
&G& &\text{: Coefficient d'élasticité au cisaillement}& &I_z& &\text{: Moment secondaire de la section transversale de l'axe Z}\\
&A& &\text{: La Coupe transversale}        & &
\end{align}

Vecteur de force d'inertie nodale de l'élément de squelette d'avion

Considérez la force d'inertie comme la force de l'objet et donnez "l'accélération" comme vecteur de nœud. Ici, la méthode la plus simple consiste simplement à répartir la force d'inertie agissant sur l'élément vers la force nodale. La force d'inertie totale agissant sur l'élément est la suivante en tant que valeur dans tout le système de coordonnées.

\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}

Par conséquent, la force nodale $ \ {\ boldsymbol {f_ {b}} } $ due à la force d'inertie peut être déterminée comme suit.

\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*}

Ici, $ \ gamma $ est le poids volumique unitaire de l'élément, $ A $ est l'aire de la section transversale de l'élément, $ L $ est la longueur de l'élément et $ g $ est l'accélération de la gravité.

Vecteur de charge de température nodale de l'élément de squelette d'avion

On suppose que le changement de température n'affecte que la force axiale de l'élément. En supposant que l'élévation de température est un changement de température positif, le vecteur de force nodale dans le système de coordonnées d'élément dû au changement de température est le suivant.

\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*}

Ici, $ EA $ est la rigidité axiale du membre, $ \ alpha $ est le coefficient de dilatation linéaire du membre, et $ \ Delta T $ est la quantité de changement de température du membre.

En transformant cela en coordonnées, le vecteur de force nodale dans le système de coordonnées global est déterminé.

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

Ici, $ [\ boldsymbol {T}] ^ T $ est la matrice de translocation de la matrice de transformation de coordonnées $ [\ boldsymbol {T}] $ (dans ce cas, égale à la matrice inverse).

Vecteur de force externe nodale de l'élément de cadre plan

En ce qui concerne le vecteur de force externe du nœud, dans le programme présenté ici, la force externe du nœud équivalent est directement entrée à partir du fichier de données d'entrée, et aucun traitement de calcul spécial n'est effectué. La charge est la force horizontale, la force verticale et le moment dans le système de coordonnées global.

Matrice de conversion de coordonnées

Les équations de rigidité d'élément obtenues jusqu'à présent sont basées sur le "système de coordonnées d'élément" fixé à l'élément. Dans une structure réelle, les éléments sont orientés dans différentes directions, il est donc nécessaire de trouver la relation avec le «système de coordonnées complet» défini de manière unifiée.

Ici, si la matrice de conversion de coordonnées qui convertit le déplacement total du système de coordonnées en déplacement du système de coordonnées de l'élément est $ [\ boldsymbol {T}] $, La matrice de rigidité $ [\ boldsymbol {K}] $ dans le système de coordonnées global et la matrice de rigidité $ [\ boldsymbol {k}] $ dans le système de coordonnées d'élément ont la relation suivante.

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

Cette relation peut être comprise à partir de la relation entre le déplacement du système de coordonnées d'élément $ \ {\ boldsymbol {u} } $ et le déplacement global du système de coordonnées $ \ {\ boldsymbol {U} } $ comme suit. ..

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

Ici, $ [\ boldsymbol {T}] $ est une matrice de conversion de coordonnées qui convertit le déplacement / force nodale du système de coordonnées total en déplacement / force nodale du système de coordonnées de l'élément.

Représentation spécifique de la matrice de conversion de coordonnées

Cosinus directionnel de l'axe du membre (axe x)

Soit les nœuds qui composent l'élément $ i $ et $ j $, et que leurs coordonnées soient $ (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}

Sous-matrice de transformation de coordonnées

$ \ boldsymbol {[t]} $ est une petite matrice (3 x 3) de la matrice de transformation de coordonnées (12 x 12) qui transforme l'ensemble du système de coordonnées en système de coordonnées membre.

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

\begin{align}
&\boldsymbol{\{x\}}=\{x,y,z\}^T  & &\text{: Coordonnée locale}\\
&\boldsymbol{\{X\}}＝\{X,Y,Z\}^T & &\text{: Coordonnée globale}
\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}

Matrice de conversion de coordonnées

\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}

Relation entre le système de coordonnées de membre (x-y-z) et le système de coordonnées global (X-Y-Z) où l'angle de code est 0

Équation de rigidité globale

D'après ce qui précède, l'équation de rigidité dans le système de coordonnées global est exprimée comme suit.

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

\begin{align}
&\{\boldsymbol{U}\}                          & & \text{Vecteur de déplacement de nœud dans le système de coordonnées global}\\
&\{\boldsymbol{F}\}                          & & \text{Nœud vecteur de force externe dans le système de coordonnées global}\\
&\{\boldsymbol{F_b}\}=\{\boldsymbol{f_b}\}   & & \text{Vecteur de force d'inertie du nœud dans le système de coordonnées global}\\
&\{\boldsymbol{F_t}\}=[\boldsymbol{T}]^T\{\boldsymbol{f_t}\}  & & \text{Vecteur de charge de température de nœud dans le système de coordonnées global}\\
&[\boldsymbol{K}]=[\boldsymbol{T}]^T[\boldsymbol{k}][\boldsymbol{T}]   & & \text{Matrice de rigidité dans le système de coordonnées global}\\
&[\boldsymbol{k}]                            & & \text{Matrice de rigidité dans le système de coordonnées des éléments}\\
&\{\boldsymbol{f_b}\}                          & & \text{Nœud vecteur de force d'inertie dans le système de coordonnées d'élément}\\
&\{\boldsymbol{f_t}\}                          & & \text{Vecteur de charge de température de nœud dans le système de coordonnées d'élément}\\
&[\boldsymbol{T}]                            & & \text{Matrice de conversion de coordonnées}
\end{align}

En résolvant l'équation ci-dessus, le déplacement du nœud $ \ {\ boldsymbol {U} } $ dans le système de coordonnées global peut être obtenu.

Force transversale de l'élément

Pour la force transversale de l'élément, le déplacement total du système de coordonnées $ \ {\ boldsymbol {U} } $ obtenu ci-dessus est converti en déplacement du système de coordonnées de l'élément par la matrice de conversion de coordonnées $ [\ boldsymbol {T}] $. Calculé en multipliant la matrice de rigidité des éléments $ [\ boldsymbol {k}] $. Autrement dit, la force transversale de l'élément $ \ {\ boldsymbol {f_ {sec}} } $ est calculée comme suit.

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

Notez que la force axiale calculée calculée ci-dessus n'inclut pas l'effet du changement de température, il est donc nécessaire de corriger l'effet de température de la force axiale comme suit.

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

Ici, $ N_1 '$ et $ N_2' $ sont les forces axiales aux nœuds 1 et 2 en considérant le changement de température, et $ N_1 $ et $ N_2 $ sont les forces axiales calculées à partir du déplacement comme solution des équations simultanées. $ EA $ est la rigidité axiale du membre, $ \ alpha $ est le coefficient de dilatation linéaire du membre, et $ \ Delta T $ est la quantité de changement de température du membre.

Programme d'analyse de la structure du squelette en trois dimensions

Aperçu du programme

Explication des conditions aux limites qui peuvent être gérées

Élément Description
Force externe nodale Spécifiez le nombre de nœuds chargés, le nombre de nœuds chargés et la valeur de charge
Changement de température nodale Spécifiez la quantité de changement de température à tous les nœuds. Soit l'élévation de température la valeur positive
Force d'inertie de l'élément Dans les données de propriétés du matériau, spécifiez l'accélération dans la direction x / y / z sous forme de rapport à l'accélération de gravité
Déplacement forcé nodal Spécifiez le nombre de nœuds, le nombre nodal et la valeur de déplacement pour donner le déplacement forcé

Configuration du programme (nom du programme: py_fem_3dfrmx.py)

   1.Chargement du module
   2.Fonction: saisie de données(def inpdata_3dfrm)
       (1)Lecture numérique de base
       (2)Déclaration de tableau
       (3)Lire les données de propriétés des matériaux
       (4)Lire le numéro de nœud de composition de l'élément et le numéro de matériau
       (5)Lire les coordonnées des nœuds et le changement de température des nœuds
       (6)Lecture des conditions aux limites
       (7)Lire les données de charge du nœud
   3.Fonction: exporter les données d'entrée(def prinp_3dfrm)
   4.Fonction: Exporter le résultat du calcul(def prout_3dfrm)
   5.Fonction: créer une matrice de rigidité d'élément(def sm_3dfrm)
   6.Fonction: créer une matrice de conversion de coordonnées(def tm_3dfrm)
   7.Fonction: création de vecteur de force d'inertie de nœud(def bfvec_3dfrm)
   8.Fonction: création de vecteur de charge de température de nœud(def tfvec_3dfrm)
   9.Fonction: calcul de la force transversale de l'élément(def calsecf_3dfrm)
  10.Fonction: routine principale(def main_3dfrm)
       (1)Obtenir le nom du fichier d'entrée / sortie à partir de la ligne de commande
       (2)Saisie des données
       (3)Exporter les données d'entrée
       (4)Assemblage de matrice rigide
       (5)Traitement des conditions aux limites (traitement du degré de déplacement connu)
       (6)Résoudre l'équation de rigidité (calcul du déplacement)
       (7)Traitement des conditions aux limites (incorporation du déplacement connu)
       (8)Calcul de la force transversale de l'élément
       (9)Exporter les résultats des calculs
       (10)Affichage d'informations telles que le temps de calcul
  11.Exécution de routine principale

Commande d'exécution et données d'entrée / sortie

Commande d'exécution d'analyse

python3 py_fem_3dfrmx.py inp.txt out.txt
Élément Description
py_fem_3dfrmx.py nom du programme Python FEM
inp.txt Nom du fichier de données d'entrée (données séparées par des blancs)
out.txt Nom du fichier de données de sortie (données séparées par des blancs)

Format de fichier de données d'entrée

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)
Élément Description
npoin, nele, nsec Nombre de nœuds, nombre d'éléments, nombre de caractéristiques transversales
npfix, nlod Nombre de nœuds contraints, nombre de nœuds chargés
E, po, A Coefficient d'élasticité de l'élément, coefficient de Poisson de l'élément, aire de la section transversale de l'élément
Ix, Iy, Iz Constante de torsion, moment secondaire de la section autour de l'axe y, moment secondaire de la section autour de l'axe z >
thêta, alpha Angle de code, coefficient de dilatation de la ligne d'élément
gamma, gkX, gkY, gkZ Poids volumique unitaire, accélération direction x / y / z (rapport de g)
node_1, node_2, isec Node 1, Node 2, Numéro de caractéristique de section
x, y, z, deltaT Coordonnée noeud x, coordonnée noeud y, coordonnée noeud z, changement de température du noeud
lp, kox, koy, koz Numéro de nœud de contrainte, x / y / z Déplacement directionnel Avec ou sans contrainte (contrainte: 1, liberté: 0)
kmx, kmy, kmz x / y / z Avec ou sans contrainte de rotation axiale (contrainte: 1, liberté: 0) < / tr>
r_disx, r_disy, rdis_z direction x / y / z déplacement connu (entrez 0 même si sans contrainte)
rrot_x, rrot_y, rrot_z x / y / z Quantité de rotation connue autour de l'axe (entrez 0 même si sans contrainte)
lp, fp_x, fp_y, fp_z Numéro de nœud de charge, charge dans la direction x, charge dans la direction y, charge dans la direction z >
mp_x, mp_y, mp_z Moment axe x (torsion), moment axe y, moment axe z

Format de fichier de données de sortie

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)
Élément Description
dis-x, dis-y, dis-z déplacement direction x, déplacement direction y, déplacement direction z
rot-x, rot-y, rot-z quantité de rotation axe x, quantité de rotation axe y, quantité de rotation axe z td>
N, Sy, Sz Force axiale, force de cisaillement dans la direction y, force de cisaillement dans la direction z
Mx, My, Mz Moment axe x (torsion), moment fléchissant axe y, moment fléchissant axe z

la programmation

Le corps principal du programme est présenté ci-dessous avec la «structure du programme».

1. Module de charge

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

2. Fonction: saisie de données (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. Fonction: écrire les données d'entrée (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. Fonction: résultat du calcul d'exportation (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. Fonction: Créer une matrice de rigidité d'élément (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. Fonction: Créer une matrice de transformation de coordonnées (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{Lorsque l'axe longitudinal de l'élément (axe x du système de coordonnées de l'élément) est parallèle à l'axe Z du système de coordonnées global}\\
&\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{Lorsque l'axe longitudinal de l'élément (axe x du système de coordonnées de l'élément) n'est pas parallèle à l'axe Z du système de coordonnées global}
\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. Fonction: création de vecteur de force d'inertie de nœud (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. Fonction: création de vecteur de charge de température de nœud (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. Fonction: calcul de la force transversale de l'élément (def calsecf_3dfrm)

$ N_1 '$ et $ N_2' $ sont des corrections de force axiale dues aux changements de température.

\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. Fonction: routine principale

Procédure de traitement de routine principale

Élément Description
npoin Nombre total de nœuds
nfree 1 liberté de nœud (6 ici)
mpfix Tableau contenant les conditions de contrainte pour chaque nœud (= 1: contrainte, = 0: libre)
fp Vecteur de charge nodale
rdis Quantité de déplacement du nœud qui contraint le déplacement
gk Matrice de rigidité globale
disg Solution (déplacement) de l'équation de rigidité globale

  • Résolution de l'équation de rigidité (calcul du déplacement) </ b>: Après conversion CSR de la matrice de rigidité totale gk </ var>, les équations linéaires simultanées sont résolues par scipy.sparse.linalg.spsolve () </ var>.

  • Traitement des conditions aux limites (incorporation du déplacement connu) </ b>: Après avoir résolu les équations linéaires simultanées, n'oubliez pas de rentrer le déplacement connu dans le terme de liberté correspondant au déplacement connu.

  • Calcul de la force transversale de l'élément </ b>: Effectuer le calcul de la force transversale pour chaque élément par la fonction calsecf_3dfrm </ em>

  • Exporter les résultats des calculs </ b>: Sortie du résultat du calcul par la fonction prout_3dfrm </ em>

  • Affichage des informations telles que le temps de calcul </ b>: La différence entre l'heure de fin du calcul et l'heure de début est sortie comme temps de traitement comme suit.

import time

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

Programme de routine principal

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. Exécution de la routine principale

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

c'est tout

Recommended Posts