LinearTriangleEl elemento LinearTriangle es un elemento finito bidimensional con coordenadas locales y globales, caracterizado por una función de forma lineal. Puede ser utilizado para problemas de esfuerzo y deformación plana. Este elemento tiene un módulo de elasticidad $E$, una relación de Poisson $\nu$ y un espesor $t$. Cada triangulo tiene tres nodos con dos grados de libertad en el plano (ux y uy), las coordenadas globales de estos nodos se denotan por $(x_i,y_i)$, $(x_j,y_j)$ y $(x_m,y_m)$ (como se muestra en la figura). El orden de los nodos para cada elemento es importante y deben listarse en sentido antihorario comenzando desde cualquier nodo.
La matriz de rigidez por elemento viene dada por:
$$ [k] = tA[B]^T[D][B] $$Donde $A$ es el área del elemento dada por:
$$ A = \frac{1}{2} \left( x_i(y_j-y_m) + x_j(y_m - y_i) + x_m(y_i - y_j) \right) $$y $[B]$ es la matriz dada por:
$$ [B] = \frac{1}{2A} \begin{bmatrix} \beta_i & 0 & \beta_j & 0 & \beta_m & 0 \\ 0 & \gamma_i & 0 & \gamma_j & 0 & \gamma_m \\ \gamma_i & \beta_i & \gamma_j & \beta_j & \gamma_m & \beta_m \\ \end{bmatrix} $$Donde $\beta_i$, $\beta_j$, $\beta_m$, $\gamma_i$, $\gamma_j$ y $\gamma_m$ están dados por:
$$ \beta_i = y_j - y_m $$$$ \beta_j = y_m - y_i $$$$ \beta_m = y_i - y_j $$$$ \gamma_i = x_m - x_j $$$$ \gamma_j = x_i - x_m $$$$ \gamma_m = x_j - x_i $$Para el caso de esfuerzo plano la matriz $D$ viene dada por:
$$ [D] = \frac{E}{1-\nu^2} \begin{bmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \frac{1-\nu}{2} \\ \end{bmatrix} $$Los esfuerzos en cada elemento se calculan mediante:
$$ \{\sigma\} = [D][B]\{u\} $$Donde $\sigma$ es el vector de esfuerzos en el plano, es decir:
$$ \sigma = \begin{Bmatrix} \sigma_x \\ \sigma_y \\ \tau_{xy} \end{Bmatrix} $$y $u$ el vector de desplazamientos en cada nodo del elemento:
$$ \{u\} = \begin{Bmatrix} ux_i \\ uy_i \\ ux_j \\ uy_j \\ ux_m \\ uy_m \end{Bmatrix}$$
In [1]:
%matplotlib inline
from nusa import *
Definimos algunas propiedades a utilizar:
In [2]:
E = 200e9 # Módulo de elasticidad
nu = 0.3 # Relación de Poisson
t = 0.1 # Espesor
Creamos los nodos y el elemento:
In [3]:
n1 = Node((0,0))
n2 = Node((1,0.5))
n3 = Node((0,1))
e1 = LinearTriangle((n1,n2,n3),E,nu,t)
Creamos el modelo y agregamos los nodos y elementos a este:
In [4]:
m = LinearTriangleModel("Modelo 01")
for n in (n1,n2,n3): m.add_node(n)
m.add_element(e1)
Aplicamos las condiciones de carga y restricciones:
In [5]:
m.add_constraint(n1, ux=0, uy=0)
m.add_constraint(n3, ux=0, uy=0)
m.add_force(n2, (5e3, 0))
En este punto podemos graficar el modelo con las condiciones impuestas, para esto utilizamos el método plot_model:
In [6]:
m.plot_model()
Finalmente resolvemos el modelo:
In [7]:
m.solve()
Podemos consultar los desplazamientos nodales:
In [8]:
for n in m.get_nodes():
print(n.ux,n.uy)
Además, utilizar las herramientas de postproceso para graficar el campo de desplazamientos en el elemento:
In [9]:
m.plot_nsol('ux')
In [10]:
import nusa.mesh as nmsh
md = nmsh.Modeler()
BB, ES = 1, 0.1
a = md.add_rectangle((0,0),(BB,BB), esize=ES)
nc, ec = md.generate_mesh()
x,y = nc[:,0], nc[:,1]
nodos = []
elementos = []
for k,nd in enumerate(nc):
cn = Node((x[k],y[k]))
nodos.append(cn)
for k,elm in enumerate(ec):
i,j,m = int(elm[0]),int(elm[1]),int(elm[2])
ni,nj,nm = nodos[i],nodos[j],nodos[m]
ce = LinearTriangle((ni,nj,nm),200e9,0.3,0.25)
elementos.append(ce)
m = LinearTriangleModel()
for node in nodos: m.add_node(node)
for elm in elementos: m.add_element(elm)
# Aplicando condiciones de frontera en los extremos
minx, maxx = min(x), max(x)
miny, maxy = min(y), max(y)
P = 100e3/((BB/ES)+1)
for node in nodos:
if node.x == minx:
m.add_constraint(node, ux=0, uy=0)
if node.x == maxx:
m.add_force(node, (P,0))
m.plot_model()
m.solve()
In [11]:
# Esfuerzo de von Mises
m.plot_nsol("seqv","Pa")
In [12]:
# Deformación unitaria en la dirección X
m.plot_nsol("exx", "")
In [13]:
%matplotlib inline
from nusa import *
import nusa.mesh as nmsh
md = nmsh.Modeler()
a = md.add_rectangle((0,0),(1,1), esize=0.1)
b = md.add_circle((0.5,0.5), 0.1, esize=0.05)
md.substract_surfaces(a,b)
nc, ec = md.generate_mesh()
x,y = nc[:,0], nc[:,1]
nodos = []
elementos = []
for k,nd in enumerate(nc):
cn = Node((x[k],y[k]))
nodos.append(cn)
for k,elm in enumerate(ec):
i,j,m = int(elm[0]),int(elm[1]),int(elm[2])
ni,nj,nm = nodos[i],nodos[j],nodos[m]
ce = LinearTriangle((ni,nj,nm),200e9,0.3,0.1)
elementos.append(ce)
m = LinearTriangleModel()
for node in nodos: m.add_node(node)
for elm in elementos: m.add_element(elm)
# Aplicando condiciones de frontera en los extremos
minx, maxx = min(x), max(x)
miny, maxy = min(y), max(y)
for node in nodos:
if node.x == minx:
m.add_constraint(node, ux=0, uy=0)
if node.x == maxx:
m.add_force(node, (10e3,0))
m.plot_model()
m.solve()
m.plot_nsol("seqv") # Esfuerzo de von Mises
In [14]:
# Element solution
m.plot_esol("sxx")
In [15]:
# generando la geometría
md = nmsh.Modeler()
g = md.geom # Para acceder a la clase SimpleGMSH
p1 = g.add_point((0,0))
p2 = g.add_point((1,0))
p3 = g.add_point((2,0))
p4 = g.add_point((2,1))
p5 = g.add_point((3,1))
p6 = g.add_point((3,2))
p7 = g.add_point((0,2))
p8 = g.add_point((0.7,1.4))
p9 = g.add_point((0.7,1.7), esize=0.1)
L1 = g.add_line(p1,p2)
L2 = g.add_circle(p3,p2,p4)
L3 = g.add_line(p4,p5)
L4 = g.add_line(p5,p6)
L5 = g.add_line(p6,p7)
L6 = g.add_line(p7,p1)
L7 = g.add_circle(p8,p9)
loop1 = g.add_line_loop(L1,L2,L3,L4,L5,L6) # boundary
loop2 = g.add_line_loop(L7)# hole
g.add_plane_surface(loop1,loop2)
nc, ec = md.generate_mesh()
x,y = nc[:,0], nc[:,1]
nodos = []
elementos = []
for k,nd in enumerate(nc):
cn = Node((x[k],y[k]))
nodos.append(cn)
for k,elm in enumerate(ec):
i,j,m = int(elm[0]),int(elm[1]),int(elm[2])
ni,nj,nm = nodos[i],nodos[j],nodos[m]
ce = LinearTriangle((ni,nj,nm),200e9,0.3,0.1)
elementos.append(ce)
m = LinearTriangleModel()
for node in nodos: m.add_node(node)
for elm in elementos: m.add_element(elm)
# Aplicando condiciones de frontera en los extremos
minx, maxx = min(x), max(x)
miny, maxy = min(y), max(y)
for node in nodos:
if node.x == minx:
m.add_constraint(node, ux=0, uy=0)
if node.x == maxx:
m.add_force(node, (10e3,1))
m.plot_model()
m.solve()
m.plot_nsol("seqv") # Esfuerzo de von Mises
In [ ]: