Cet article est écrit en Java en référence à Première méthode des éléments finis de Java publiée par Corona. C'est un programme réalisé avec python. De plus, cette fois, c'est au programme de créer la matrice de rigidité des éléments.
La méthode des éléments finis est utilisée pour simuler la contrainte et la déformation de la dynamique des solides. En mécanique du solide, un objet est divisé en plusieurs petites zones (éléments), et l'état de chaque élément est calculé pour déterminer l'état de déformation de la substance entière. L'état d'un objet est exprimé par trois équations différentielles (équation de contrainte parallèle, expression relationnelle déformation-déplacement, expression relationnelle contrainte-déformation (expression de construction)), mais dans la méthode des éléments finis, ces équations différentielles sont exprimées par des expressions intégrales équivalentes ( En se substituant au principe du travail virtuel), le déplacement et la déformation sont approximativement obtenus. La représentation intégrale remplace les équations différentielles par les équations rigides, qui sont des équations linéaires simultanées. De plus, le flux de calcul de la méthode des éléments finis consiste à diviser l'objet en chaque élément, à établir une équation de rigidité (matrice de rigidité d'élément) pour chaque élément, puis à intégrer les équations de rigidité de tous les éléments pour faire l'équation de rigidité globale (globale). Matrice de rigidité) est assemblée et une solution approximative est calculée.
Pour simplifier, un point final d'un objet avec une section transversale de $ A \ \ rm {[L ^ {2}]} $ et une longueur de $ l \ \ rm {[L]} $ est fixe et l'autre point final Considérons le cas où une charge de traction unidimensionnelle et uniforme $ F \ \ rm {[F]} $ fonctionne. (Ici, laissez $ \ rm {[L]} $ représenter la dimension de la longueur et $ \ rm [F] $ la dimension de la force. En utilisant une telle notation, par exemple, la pression est une aire unitaire. Puisqu'il s'agit d'une force de frappe, elle peut être écrite comme la dimension de force divisée par le carré de la dimension de longueur, c'est-à-dire $ \ rm {[FL ^ {-2}]} $. La dimension du paramètre est clarifiée. Cela permet de comprendre plus facilement ce que ce paramètre représente et comment il se rapporte à d'autres paramètres.) De plus, cette substance est transformée en un corps élastique (une substance qui se déforme lorsqu'une force est appliquée et reprend sa forme d'origine lorsque la force est arrêtée), et son rapport de Young est $ E \ \ rm {[FL ^ {-2}]}. Supposons que ce soit $. Le rapport de Young est le coefficient de proportionnalité appliqué à la déformation $ \ varepsilon $ [dimension] de la loi de Hook selon laquelle "la contrainte est proportionnelle à la déformation" qui se maintient dans un corps élastique. (Similaire à la relation entre la charge (ou la contrainte) et le déplacement lorsqu'une force est appliquée à un ressort, mais sachez que les dimensions sont différentes.)
En supposant cela, l'équilibre des forces dans cet objet peut être écrit comme suit, en supposant que la contrainte générée dans l'objet est $ \ sigma \ \ rm {[FL ^ {-2}]} $.
$ F = \ sigma A $ (équilibre de force, équation d'équilibre de contrainte) (1 ')
De plus, en supposant que la contrainte est $ \ sigma $, la loi de Hook peut être écrite comme suit en utilisant le taux de Young $ E $ et la déformation $ \ varepsilon $.
$ \ Sigma = E \ \ varepsilon $ (propriétés du matériau, équation contrainte-déformation) (2 ')
La déformation $ \ varepsilon $ peut être exprimée comme suit, où l'allongement de l'objet sous charge de traction est $ u \ rm {[L]} $.
$ \ Varepsilon = u / L $ (changement de forme, relation déformation-déplacement) (3 ')
Ces équations (1 ') - (3') sont des cas particuliers d'équations différentielles de base. Compte tenu de la situation plus générale, les trois équations différentielles de base suivantes sont dérivées.
$ d (\ sigma A) / dx = 0 $ (équilibre de force, équation d'équilibre de contrainte) (1) $ \ Sigma = E \ \ varepsilon $ (propriétés du matériau, équation contrainte-déformation) (2) $ \ Varepsilon = du / dx $ (changement de forme, expression relationnelle déformation-déplacement) (3)
Afin de déterminer analytiquement la déformation générée dans un corps élastique, il est nécessaire de résoudre les équations différentielles des équations (1) - (3).
Le principe du travail virtuel est expliqué comme suit. Veuillez vous référer aux références car l'explication détaillée sera longue.
Lorsque l'objet est déformé et dans un état équilibré, le travail virtuel effectué par la force externe est égal au travail virtuel effectué par la force interne.
Le principe du travail virtuel peut être exprimé par une formule mathématique comme suit.
À propos du déplacement virtuel arbitraire et de la déformation virtuelle avec la force équilibrée
\int_V \sigma \delta\varepsilon dV = F\delta u Est établi.
Où $ \ delta \ varepsilon $ est la déformation virtuelle et $ \ delta u $ est le déplacement virtuel.
Maintenant, la méthode des éléments finis est prête. Dans la méthode des éléments finis, l'objet est d'abord divisé en petites zones, «éléments». Ensuite, appliquez le principe du travail virtuel à chaque élément et créez une équation de rigidité d'élément (matrice de rigidité d'élément) pour chaque élément. Considérons le principe du travail virtuel pour chaque élément lorsqu'un objet unidimensionnel est divisé. Puisqu'un élément a deux nœuds (joints entre les éléments), l'un est le nœud 1 (coordonnée $ x_1 $) et l'autre est le nœud 2 (coordonnée $ x_2 $), et la force interne agissant sur chacun est $ f_1. Soit les déplacements de points respectivement $ u_1 $ et $ u_2 $ pour $ et $ f_2 $. Si la contrainte ici est $ \ sigma $ et le déplacement virtuel de chaque nœud est $ \ delta u_1 $, $ \ delta u_2 $, le principe du travail virtuel de cet élément est
$\int_{x_1}^{x_2} \sigma \delta\varepsilon A dx = f_1 \cdot \delta u_1 + f_2 \cdot \delta u_2 $ (4)
Peut être écrit.
Supposons maintenant qu'à l'intérieur de chaque élément, le déplacement $ u $ change avec une fonction linéaire de coordonnées $ x $.
Ce déplacement doit être égal à $ u_1 $, $ u_2 $ aux nœuds 1 et 2, donc
En résumant ceci pour $ a $ et $ b $,
Ce sera. Où $ l = x_2-x_1 $ est la longueur de l'élément. En utilisant (8) et (9), (5) est
De plus, en réécrivant l'équation (10),
Par conséquent, il était possible de réécrire l'équation avec les déplacements nodaux $ u_1 $ et $ u_2 $ comme nombres inconnus. De plus, $ N_1 (x) $ et $ N_2 (x) $ sont appelés fonctions de forme. De plus, l'équation (11) peut être exprimée dans une matrice,
$u(x) = {\bf N}^{{\rm T}} {\bf U} $ (12)
Peut être écrit. Où $ {\ bf N} $ et $ {\ bf U} $ sont
{\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)
C'est un vecteur.
En utilisant l'équation de déplacement (11) et l'équation de relation déformation-déplacement (3), la déformation $ \ varepsilon $ est
\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)
Peut être écrit. C'est aussi une notation matricielle,
$\varepsilon = {\bf B}^{{\rm T}} {\bf U} $ (15)
Peut être réécrit comme. Ici, $ {\ bf B} $ est appelé la matrice $ {\ bf B} $, ou matrice de déformation-déplacement, et peut s'écrire comme suit.
{\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)
La contrainte $ \ sigma $ est calculée à partir des équations de relation contrainte-déformation (2) et (15).
$\sigma = E \varepsilon = E{\bf B}^{{\rm T}} {\bf U} $ (17)
Peut être exprimé comme.
A partir des équations (12) et (15), le déplacement virtuel $ \ delta u $ et la déformation virtuelle $ \ delta \ varepsilon $ sont exprimés en notation matricielle.
$\delta u = {\bf N}^{{\rm T}} \delta{\bf U} $ (18) $\delta \varepsilon = {\bf B}^{{\rm T}} \delta{\bf U} $ (19)
Peut être écrit comme.
Considérons maintenant l'équation (4) du principe du travail virtuel.
$\int_{x_1}^{x_2} \sigma \delta\varepsilon A dx = f_1 \cdot \delta u_1 + f_2 \cdot \delta u_2 $ (4)
Tout d'abord, réécrire l'intérieur de l'intégrale sur le côté gauche de l'équation (4) avec la notation matricielle,
Par conséquent, si cela est appliqué au côté gauche de l'équation (4),
$\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)
Considérant que l'aire transversale $ A $ et le taux jeune $ E $ sont constants dans les éléments de l'équation (21),
$\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)
Où $ V $ est $ A \ int_ {x_1} ^ {x_2} dx = El = V $, qui représente le volume de l'élément. De plus, si vous définissez $ {\ bf K} = E V {\ bf B} {\ bf B} ^ {{\ rm T}} $, vous pouvez le réécrire comme suit.
$\int_{x_1}^{x_2} \sigma \delta\varepsilon A dx = {\bf U}^{{\rm T}} {\bf K} \delta{\bf U} $ (23)
$ {\ bf K} $ est appelé la matrice de rigidité des éléments.
Ensuite, si le côté droit de l'équation (4) est exprimé dans une matrice,
Où $ {\ bf F} $ est:
{\bf F} = \left(
\begin{array}{c}
f_1 \\
f_2
\end{array}
\right)
(25)
De cette manière, l'équation (4) du principe du travail virtuel peut être réécrite comme suit à l'aide des équations (23) et (24).
De plus, si la formule ci-dessus est résumée par $ \ delta {\ bf U} $,
Ici, puisque le principe du travail virtuel vaut pour tout déplacement virtuel $ \ delta {\ bf U} $,
Si vous prenez cette inversion,
$ {\bf K}{\bf U} = {\bf F}$ (29)
Ce sera. L'équation (29) est une équation appelée équation de rigidité de l'élément. Comme mentionné au début, après avoir formulé cette équation de rigidité d'élément pour chaque élément, les équations de rigidité de tous les éléments sont intégrées pour construire l'équation de rigidité globale (matrice de rigidité globale) et calculer la solution approximative. Ce sera.
Eh bien, j'ai créé un programme pour calculer une équation de rigidité d'élément unidimensionnelle avec python. La bibliothèque utilisée est numpy et le code est ci-dessous. A partir des conditions données, vous pouvez calculer la matrice $ {\ bf B} $ et la matrice de rigidité d'élément $ {\ bf K} $, contrainte pour un élément.
import numpy as np
class Element:
#Nombre de contacts
NumberOfNodeInElement = 2
def __init__(self):
#Longueur d'élément
self.length = 0
#Section transversale de l'élément
self.area = 0
#Module d'Young
self.young = 0
#Numéros de contact qui composent l'élément
self.ix = np.zeros(self.NumberOfNodeInElement)
#Coordonner les valeurs des contacts qui composent l'élément
self.xe = np.zeros(self.NumberOfNodeInElement)
#Déplacement de contact de l'élément
self.ue = np.zeros(self.NumberOfNodeInElement)
#Bmatrix
self.Bmatrix = np.zeros(self.NumberOfNodeInElement)
#Matrice de rigidité des éléments
self.Kmatrix = np.zeros((self.NumberOfNodeInElement,self.NumberOfNodeInElement))
#Calcul Bmatrix
def makeBmatrix(self):
self.length = self.xe[1]-self.xe[0]
self.Bmatrix = np.array([[-1/self.length, 1/self.length]])
#Calcul de la matrice de rigidité des éléments
def makeElementStiffness(self):
self.volume = self.area*self.length
self.Kmatrix = self.young*self.volume*np.dot(self.Bmatrix.T,self.Bmatrix)
#Calcul du stress
def stressCalculation(self):
return self.young*np.vdot(self.Bmatrix,self.ue)
Par exemple, le calcul est effectué pour un élément dont les coordonnées du nœud 1 et du nœud 2 sont $ x_1 = 50 $ (mm), $ x_2 = 150 $ (mm), aire de la section transversale 100 $ mm ^ 2 $ et taux de Young 200 GPa. Quand
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()
Le résultat du calcul dans ce cas est le suivant.
#Bmatrix
[[-0.01 0.01]]
#Bmatrix
[[ 200000. -200000.]
[-200000. 200000.]]
#stress
30.000000000000004
Cette fois, j'ai expliqué jusqu'à la matrice de rigidité des éléments unidimensionnels de la méthode des éléments finis et j'ai introduit le programme de fabrication de matrices de rigidité des éléments. La prochaine fois, nous assemblerons et calculerons l'équation de rigidité globale (matrice de rigidité globale).
Shigeru Nagagi (2010) "La première méthode des éléments finis de Java" Corona
Recommended Posts