PFG · Justificación teórica

Diseño generativo

📄 Descargar PDF original

1.1 Diseño paramétrico

Como se explica en el libro (Tedeschi y Andreani 2014) El arquitecto italiano Luigi Moretti Inventó en 1939 el término arquitectura paramétrica. A lo largo del siglo XIX Surgió un nuevo enfoque para el diseño que sustituía el proceso aditivo de geometrías “form-making” por el enfoque en búsqueda de formas “form-finding” pioneros como Gaudí (1852-1926), Isler (1926-2009), Otto (1925-) y Musmeci (1926-1981) en la naturaleza procesos de autoformación como una manera de organizar las construcciones.

En el ámbito del diseño asistido por ordenador (DAO) se pueden crear diseños mediante dibujo a mano en un plano virtual o directamente mediante un programa informático que defina los parámetros de la geometría. A esto último se lo conoce como script. Hoy en día todos los programas comerciales de diseño permiten la interacción de scripts en uno o varios lenguajes informáticos. Los más extendidos son Visual Basic, Python, C-Sharp.

A medio camino entre las dos técnicas antes descritas se encuentran los entornos gráficos de programación como Grashopper para Rhino3D, Dynamo de Autodesk, Houdini o Blender. En ellas se diseña el propio código con diagramas de flujo. O en herramientas como Houdini, se integras varias capas del flujo de trabajo en el mismo diagrama, capa de geometría, capa de materiales, capa de renderizados.

1.2 Método de los elementos finitos

El método de los elementos finitos, a partir de ahora MEF, nació en los años 50 como una generalización de la idea básica del cálculo matricial de estructuras. En el cálculo matricial dividimos la estructura en barras dentro de las cuales se conoce la solución exacta en función de unos coeficientes que se hacen coincidir con los movimientos de los nodos de los extremos. Para la resolución de un problema por MEF, conoceremos las fuerzas exteriores aplicadas en algunos nodos, {F} se conoce como vector de cargas y aglutina las fuerzas exteriores y las reacciones que actúan sobre el sólido. {u}, denominado vector de desplazamientos nodales representa de forma ordenada los desplazamientos de los nodos. Como resultado obtendremos la ecuación matricial {F} = [K]{u} donde [K] es la matriz global de rigidez del sistema. La obtención de dicha matriz se hace mediante el ensamblaje de las distintas matrices de rigidez de los elementos que constituyen la estructura donde los nodos comunes entre partes discretizadas se suman en el mismo puesto de la matriz global.

El carácter general del método ha permitido su aplicación en más ámbitos que la elasticidad, transmisión de calor, mecánica de fluidos, electromagnetismo, entre otros.

Como se explica en (Beltran 1999) “El MEF puede entenderse también desde el punto de vista matemático como un procedimiento de obtener soluciones aproximadas al problema planteado de forma débil. La idea es aproximar los espacios 𝒮\mathcal{S} y 𝒱\mathcal{V}, que son de dimensión infinita, por espacios 𝒮h\mathcal{S}^{h} y 𝒱h\mathcal{V}^{h}de dimensión finita.

𝒮h⊂𝒮\mathcal{S}^{h}\mathcal{\subset S} y 𝒱h⊂𝒱\mathcal{V}^{h}\mathcal{\subset V}

Una vez definidos estos dos espacios 𝒮h\mathcal{S}^{h} y 𝒱h\mathcal{V}^{h} mediante sus correspondientes bases de funciones, se resuelve la forma débil del problema dentro ellos. Es decir, se busca una función uh∈𝒮hu^{h} \in \mathcal{S}^{h} que haga que se cumpla:

∫01w,xhu,xhdx=∫01whfdx+wh(0)h∀wh∈𝒱h\int_{0}^{1}\mspace{2mu} w_{,x}^{h}u_{,x}^{h}dx = \int_{0}^{1}\mspace{2mu} w^{h}fdx + w^{h}(0)h\ \ \ \ \ \ \ \ \ \ \ \ \ \ \forall\ w^{h} \in \mathcal{V}^{h}

Como en el caso del método de Galerkin, la aportación del MEF moderno dentro de esta interpretación matemática del mismo consiste en una forma sistemática de construir los espacios 𝒮h\mathcal{S}^{h} y 𝒱h\mathcal{V}^{h}”

Ahondar más en la teoría y bases matemáticas del método sale de los objetivos del presente trabajo, a modo de resumen. Tras discretizar la geometría a estudiar en una nube de puntos que se conoce como malla, de la cual se conocen las coordenadas de todos sus puntos. Deben conocerse las propiedades del material y establecer un modelo de comportamiento, en nuestro caso el modelo constitutivo es el elástico deformable que necesitará de parámetros del material como módulo de elasticidad, coeficiente de Poisson. Aplicando el Principio de los Trabajos Virtuales, podemos establecer que, para unos desplazamientos virtuales de los nudos del elemento, {uδ}\left\{ u^{\delta} \right\}, el trabajo realizado por las fuerzas reales que actúan en los nudos del elemento, {F}\{ F\}, debe ser igual al trabajo realizado por las tensiones internas {σ}\{\sigma\} reales que actúan dentro del elemento teniendo en cuenta el campo de deformaciones virtuales {εδ}\left\{ \varepsilon^{\delta} \right\} que aparecería en el interior del elemento.

Volumen tetraédrico con vértices identificados de manera consistente. En este caso, visto desde p, ordenados en el sentido inverso a las agujas de un reloj. ijmp

El sistema de referencia elegido, plano cartesiano, cilíndrico, tridimensional cartesiano influirá en la matriz de rigidez local que rige el comportamiento de cada elemento finito. Las tensiones en el interior del elemento se podrían obtener a partir de las deformaciones utilizando las ecuaciones constitutivas. A continuación particularizamos el problema FEM para el caso tridimensional isotrópico tal como se describe en el capítulo 6 de (Zienkiewicz, Taylor y Zhu 2013)

El desplazamiento de cada punto se expresa en forma de vector en coordenadas cartesianas:

𝒖={uvw}\mathbf{u} = \left\{ \begin{aligned} & u \\ & v \\ & w \end{aligned} \right\}

Expresando los desplazamientos en forma de polinomio

u=α1+α2x+α3y+α4zu = \alpha_{1} + \alpha_{2}x + \alpha_{3}y + \alpha_{4}z

Igualando los desplazamientos en los nodos tendremos cuatro ecuaciones del tipo

ui=α1+α2xi+α3yi+α4ziu_{i} = \alpha_{1} + \alpha_{2}x_{i} + \alpha_{3}y_{i} + \alpha_{4}z_{i}

Escribiendo la solución en forma de determinantes tendremos

u=16V[(ai+bix+ciy+diz)ui+(aj+bjx+cjy+djz)uj+(am+bmx+cmy+dmz)um+(ap+bpx+cpy+dpz)up]\begin{aligned} u = & \frac{1}{6V}\left\lbrack \left( a_{i} + b_{i}x + c_{i}y + d_{i}z \right)u_{i} + \left( a_{j} + b_{j}x + c_{j}y + d_{j}z \right)u_{j} \right.\ \\ & \left. \ + \left( a_{m} + b_{m}x + c_{m}y + d_{m}z \right)u_{m} + \left( a_{p} + b_{p}x + c_{p}y + d_{p}z \right)u_{p} \right\rbrack \end{aligned}

donde

6V=det⁡|1xiyizi1xjyjzj1xmymzm1xpypzp|6V = \det\left| \begin{matrix} 1 & x_{i} & y_{i} & z_{i} \\ 1 & x_{j} & y_{j} & z_{j} \\ 1 & x_{m} & y_{m} & z_{m} \\ 1 & x_{p} & y_{p} & z_{p} \end{matrix} \right|

en el que, por cierto, el valor VV representa el volumen del tetraedro. Al expandir los otros determinantes relevantes en sus cofactores, tenemos

ai=det⁡|xjyjzjxmymzmxpypzp|bi=−det⁡|1yjzj1ymzm1ypzp|a_{i} = \det\left| \begin{matrix} x_{j} & y_{j} & z_{j} \\ x_{m} & y_{m} & z_{m} \\ x_{p} & y_{p} & z_{p} \end{matrix} \right|\ \ b_{i} = - \det\left| \begin{matrix} 1 & y_{j} & z_{j} \\ 1 & y_{m} & z_{m} \\ 1 & y_{p} & z_{p} \end{matrix} \right|

ci=−det⁡|xj1zjxm1zmxp1zp|di=−det⁡|xjyj1xmym1xpyp1|c_{i} = - \det\left| \begin{matrix} x_{j} & 1 & z_{j} \\ x_{m} & 1 & z_{m} \\ x_{p} & 1 & z_{p} \end{matrix} \right|\ \ d_{i} = - \det\left| \begin{matrix} x_{j} & y_{j} & 1 \\ x_{m} & y_{m} & 1 \\ x_{p} & y_{p} & 1 \end{matrix} \right|

con las otras constantes definidas por el intercambio cíclico de los subíndices en el orden ii, j,m,pj,m,p.

El orden de los números nodales i,j,m,pi,j,m,p debe seguir una regla de "mano derecha" obvia de la Fig. 6.1. En esto, los tres primeros nodos se numeran en sentido contrario a las agujas del reloj cuando se ven desde el último.

El desplazamiento de cada elemento está definido por las 12 componentes de los nodos como

𝒂e={𝒂i𝒂j𝒂m𝒂p}\mathbf{a}^{e} = \left\{ \begin{array}{r} \mathbf{a}_{i} \\ \mathbf{a}_{j} \\ \mathbf{a}_{m} \\ \mathbf{a}_{p} \end{array} \right\} con 𝒂i={uiviwi}\mathbf{a}_{i} = \left\{ \begin{array}{r} u_{i} \\ v_{i} \\ w_{i} \end{array} \right\} 𝒂j={ujvjwj}\ \ \mathbf{a}_{j} = \left\{ \begin{array}{r} u_{j} \\ v_{j} \\ w_{j} \end{array} \right\} 𝒂m={umvmwm}\ \ \mathbf{a}_{m} = \left\{ \begin{array}{r} u_{m} \\ v_{m} \\ w_{m} \end{array} \right\} 𝒂p={upvpwp}\ \ \mathbf{a}_{p} = \left\{ \begin{array}{r} u_{p} \\ v_{p} \\ w_{p} \end{array} \right\}

Como el polinomio con el que definimos los desplazamientos es lineal, ahora podemos conocer el desplazamiento de cualquier punto interior a un elemento con

𝒖=[𝑰Ni,𝑰Nj,𝑰Nm,𝑰Np]𝒂e=𝑵𝒂e\mathbf{u} = \left\lbrack \mathbf{I}N_{i},\mathbf{I}N_{j},\mathbf{I}N_{m},\mathbf{I}N_{p} \right\rbrack\mathbf{a}^{e} = \mathbf{Na}^{e}

Con la función de forma definida como Ni=ai+bix+ciy+diz6VN_{i} = \frac{a_{i} + b_{i}x + c_{i}y + d_{i}z}{6V}

Nj=aj+bjx+cjy+djz6VN_{j} = \frac{a_{j} + b_{j}x + c_{j}y + d_{j}z}{6V}

Nm=am+bmx+cmy+dmz6VN_{m} = \frac{a_{m} + b_{m}x + c_{m}y + d_{m}z}{6V}

Np=ap+bpx+cpy+dpz6VN_{p} = \frac{a_{p} + b_{p}x + c_{p}y + d_{p}z}{6V}

Y siendo I la matriz identidad de 3x3.

Si construimos el vector de deformaciones a partir de las deformaciones lineales y angulares del tensor de deformaciones tridimensional de la siguiente forma:

ε={εxεyεzγxyγyzγzx}={∂u∂x∂v∂y∂w∂z∂u∂y+∂v∂x∂v∂z+∂w∂y∂w∂x+∂u∂z}\varepsilon = \left\{ \begin{array}{r} \varepsilon_{x} \\ \varepsilon_{y} \\ \varepsilon_{z} \\ \gamma_{xy} \\ \gamma_{yz} \\ \gamma_{zx} \end{array} \right\} = \left\{ \begin{array}{r} \frac{\partial u}{\partial x} \\ \frac{\partial v}{\partial y} \\ \frac{\partial w}{\partial z} \\ \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \\ \frac{\partial v}{\partial z} + \frac{\partial w}{\partial y} \\ \frac{\partial w}{\partial x} + \frac{\partial u}{\partial z} \end{array} \right\}

Podemos definir una matriz de deformaciones [S] tal que 𝜺=𝑺𝒖\mathbf{\varepsilon} = \mathbf{Su}, en la cual se comprueba que

𝜺=𝑺𝑵𝒂e=𝑩𝒂e=[𝑩i,𝑩j,𝑩m,𝑩p]𝒂e\mathbf{\varepsilon} = \mathbf{SNa}^{e} = \mathbf{Ba}^{e} = \left\lbrack \mathbf{B}_{i},\mathbf{B}_{j},\mathbf{B}_{m},\mathbf{B}_{p} \right\rbrack\mathbf{a}^{e}

Con

𝑩i=[∂Ni∂x,0,00,∂Ni∂y,00,0,∂Ni∂z∂Ni∂y,∂Ni∂x,00,∂Ni∂z,∂Ni∂y∂Ni∂z,0,∂Ni∂x]=16V[bi,0,00,ci,00,0,dici,bi,00,di,cidi0,bi]\mathbf{B}_{i} = \begin{bmatrix} \frac{\partial N_{i}}{\partial x}, & 0, & 0 \\ 0, & \frac{\partial N_{i}}{\partial y}, & 0 \\ 0, & 0, & \frac{\partial N_{i}}{\partial z} \\ \frac{\partial N_{i}}{\partial y}, & \frac{\partial N_{i}}{\partial x}, & 0 \\ 0, & \frac{\partial N_{i}}{\partial z}, & \frac{\partial N_{i}}{\partial y} \\ \frac{\partial N_{i}}{\partial z}, & 0, & \frac{\partial N_{i}}{\partial x} \end{bmatrix} = \frac{1}{6V}\begin{bmatrix} b_{i}, & 0, & 0 \\ 0, & c_{i}, & 0 \\ 0, & 0, & d_{i} \\ c_{i}, & b_{i}, & 0 \\ 0, & d_{i}, & c_{i} \\ d_{i} & 0, & b_{i} \end{bmatrix}

𝑩j=[∂Nj∂x,0,00,∂Nj∂y,00,0,∂Nj∂z∂Nj∂y,∂Nj∂x,00,∂Nj∂z,∂Nj∂y∂Nj∂z,0,∂Nj∂x]=16V[bj,0,00,cj,00,0,djcj,bj,00,dj,cjdj0,bj]\mathbf{B}_{j} = \begin{bmatrix} \frac{\partial N_{j}}{\partial x}, & 0, & 0 \\ 0, & \frac{\partial N_{j}}{\partial y}, & 0 \\ 0, & 0, & \frac{\partial N_{j}}{\partial z} \\ \frac{\partial N_{j}}{\partial y}, & \frac{\partial N_{j}}{\partial x}, & 0 \\ 0, & \frac{\partial N_{j}}{\partial z}, & \frac{\partial N_{j}}{\partial y} \\ \frac{\partial N_{j}}{\partial z}, & 0, & \frac{\partial N_{j}}{\partial x} \end{bmatrix} = \frac{1}{6V}\begin{bmatrix} b_{j}, & 0, & 0 \\ 0, & c_{j}, & 0 \\ 0, & 0, & d_{j} \\ c_{j}, & b_{j}, & 0 \\ 0, & d_{j}, & c_{j} \\ d_{j} & 0, & b_{j} \end{bmatrix}

𝑩m=[∂Nm∂x,0,00,∂Nm∂y,00,0,∂Nm∂z∂Nm∂y,∂Nm∂x,00,∂Nm∂z,∂Nm∂y∂Nm∂z,0,∂Nm∂x]=16V[bm,0,00,cm,00,0,dmcm,bm,00,dm,cmdm0,bm]\mathbf{B}_{m} = \begin{bmatrix} \frac{\partial N_{m}}{\partial x}, & 0, & 0 \\ 0, & \frac{\partial N_{m}}{\partial y}, & 0 \\ 0, & 0, & \frac{\partial N_{m}}{\partial z} \\ \frac{\partial N_{m}}{\partial y}, & \frac{\partial N_{m}}{\partial x}, & 0 \\ 0, & \frac{\partial N_{m}}{\partial z}, & \frac{\partial N_{m}}{\partial y} \\ \frac{\partial N_{m}}{\partial z}, & 0, & \frac{\partial N_{m}}{\partial x} \end{bmatrix} = \frac{1}{6V}\begin{bmatrix} b_{m}, & 0, & 0 \\ 0, & c_{m}, & 0 \\ 0, & 0, & d_{m} \\ c_{m}, & b_{m}, & 0 \\ 0, & d_{m}, & c_{m} \\ d_{m} & 0, & b_{m} \end{bmatrix}

𝑩p=[∂Np∂x,0,00,∂Np∂y,00,0,∂Np∂z∂Np∂y,∂Np∂x,00,∂Np∂z,∂Np∂y∂Np∂z,0,∂Np∂x]=16V[bp,0,00,cp,00,0,dpcp,bp,00,dp,cpdp0,bp]\mathbf{B}_{p} = \begin{bmatrix} \frac{\partial N_{p}}{\partial x}, & 0, & 0 \\ 0, & \frac{\partial N_{p}}{\partial y}, & 0 \\ 0, & 0, & \frac{\partial N_{p}}{\partial z} \\ \frac{\partial N_{p}}{\partial y}, & \frac{\partial N_{p}}{\partial x}, & 0 \\ 0, & \frac{\partial N_{p}}{\partial z}, & \frac{\partial N_{p}}{\partial y} \\ \frac{\partial N_{p}}{\partial z}, & 0, & \frac{\partial N_{p}}{\partial x} \end{bmatrix} = \frac{1}{6V}\begin{bmatrix} b_{p}, & 0, & 0 \\ 0, & c_{p}, & 0 \\ 0, & 0, & d_{p} \\ c_{p}, & b_{p}, & 0 \\ 0, & d_{p}, & c_{p} \\ d_{p} & 0, & b_{p} \end{bmatrix}

En consecuencia, el vector de tensiones tendrá la forma:

𝝈={σxσyσzτxyτyzτzx}=𝑫(ε−ε0)+𝝈0\mathbf{\sigma} = \left\{ \begin{array}{r} \sigma_{x} \\ \sigma_{y} \\ \sigma_{z} \\ \tau_{xy} \\ \tau_{yz} \\ \tau_{zx} \end{array} \right\} = \mathbf{D}\left( \varepsilon - \varepsilon_{0} \right) + \mathbf{\sigma}_{0}

Y recordando que suponemos comportamiento isótropo podemos definir la matriz [D] en función del módulo de Young y el ratio de Poisson

𝑫=E(1+ν)(1−2ν)[1−ν,ν,ν,0,0,01−ν,ν,0,0,01−ν,0,0,0(1−2ν)/2,0,0 Sym. (1−2ν)/2,0(1−2ν)/2]\mathbf{D} = \frac{E}{(1 + \nu)(1 - 2\nu)}\begin{bmatrix} 1 - \nu, & \nu, & \nu, & 0, & 0, & 0 \\ & 1 - \nu, & \nu, & 0, & 0, & 0 \\ & & 1 - \nu, & 0, & 0, & 0 \\ & & & (1 - 2\nu)/2, & 0, & 0 \\ & \text{ Sym. } & & & (1 - 2\nu)/2, & 0 \\ & & & & & (1 - 2\nu)/2 \end{bmatrix}

La matriz de rigidez definida por la relación general

δ𝒂eT𝒒e=δ𝒂eT(∫Ve𝑩T𝝈d(vol)−∫Ve𝑵T𝒃d(vol))\delta\mathbf{a}^{eT}\mathbf{q}^{e} = \delta\mathbf{a}^{eT}\left( \int_{V^{e}}^{}\mspace{2mu}\mathbf{B}^{T}\mathbf{\sigma}d(vol) - \int_{V^{e}}^{}\mspace{2mu}\mathbf{N}^{T}\mathbf{b}d(vol) \right) (aplicación PTV al VeV^{e})

ahora se puede integrar explícitamente ya que los componentes de deformación y tensión son constantes dentro del elemento.

La submatriz ij de la matriz global de rigidez será una matriz de tres por tres definida como

𝑲ije=𝑩iT𝑫𝑩jVe\mathbf{K}_{ij}^{e} = \mathbf{B}_{i}^{T}\mathbf{D}\mathbf{B}_{j}V^{e}

Donde VeV^{e} representa el volumen del tetraedro

Las fuerzas nodales debidas a las acciones serán de la forma:

𝒇ie=−𝑩iT𝑫ε0Ve\mathbf{f}_{i}^{e} = - \mathbf{B}_{i}^{T}\mathbf{D}\varepsilon_{0}V^{e}

Nótese que el FEM solo tiene en cuenta la fuerza aplicada en los nodos de la malla, cuando en un programa informático de análisis establecemos una acción tipo fuerza distribuida, a nivel interno de cálculo, esta fuerza solo se aplicará como fuerzas puntuales.

3.3 Algoritmos recursivos para el diseño

Citando a Francisco Perales López, “Se dice que un sistema es recursivo cuando está parcial o completamente definido en términos de sí mismo.” En el capítulo 3 de su libro se explica que la definición recurrente de un conjunto muestra el modo de construir cualquier elemento de él mismo a partir de ciertos elementos base y de ciertas reglas específicas de construcción. Esto es extrapolable a muchos campos de conocimiento. Como ejemplo más que conocido de tipología estructural, el domo de Richard Buckminster Fuller, que se define como:

  • Base: Icosaedro

  • Recurrencia: Proyección esférica de los vértices y puntos equidistantes de las aristas

Es un ejemplo perfecto de implementación de scripts para el diseño. Se le denomina iteración u orden a cada repetición del algoritmo.

1.3 Diseño generativo

Dado un sólido elástico isótropo con una forma conocida. Conociendo en módulo de Young y el de Poisson del material. Los programas de simulación mediante la técnica de análisis por elementos finitos (FEM) pueden calcular, para unas restricciones y solicitaciones dadas, el desplazamiento y el estado de tensiones de cada punto de la discretización del sólido. Como se explica en (Zienkiewicz, Taylor y Zhu 2013), para el caso de solidos en tres dimensiones se discretiza la geometría a un conjunto de tetraedros. Conocida la matriz de rigidez del elemento se puede construir la matriz global de la estructura. El orden de la matriz global es de decenas de miles de filas y la tarea de invertir dicha matriz sólo es posible mediante herramientas de computación. Para materiales isótropos, actualmente, el flujo de trabajo desde la discretización de la figura a la resolución del problema para unas solicitaciones concretas puede automatizarse íntegramente.

El diseño generativo es un sistema recursivo que podemos definir como una base y una recurrencia.

La base es un sistema que consta de:

  • Geometrías a conectar.

  • Geometrías a evitar.

  • Acciones (conjunto de pares solicitaciones mecánicas y restricciones).

  • Base de datos de propiedades mecánicas de materiales isótropos.

  • Objetivos.

La recurrencia

  • Paso 0: (solo se ejecuta en la primera iteración) Se rellena todo el volumen entre las geometrías a conectar y se resta las geometrías a evitar.

  • Paso 1: Se resuelve el problema de por FEM de cada acción por separado y se mapean a un sumatorio de estado total de deformación y de tensión de las distintas acciones.

  • Paso 2: Se elimina material de donde está menos solicitado, o se añade si se excede un valor crítico según un criterio de agotamiento elegido.

El algoritmo se repite un número indeterminado de veces hasta que converge, es decir, hasta que todo el material que queda está trabajando por causa de alguna de las acciones. El proceso se repite para todos los materiales que se han especificado en la base de datos de la base.

Nótese que en el paso dos de la recurrencia hemos sido deliberadamente ambiguos. El criterio de agotamiento que suponemos para averiguar la solicitación mínima crítica de cada punto una vez resueltas las acciones no se especifica en la herramienta que usamos en el presente trabajo (Fusion 360). Lo usual en materiales dúctiles es el criterio de Von Mises y Hencky. Citando la definición de (Rodríguez-Avial Llardent 2013) “El agotamiento se produce cuando la energía interna de distorsión por unidad de volumen alcanza el valor correspondiente a la fluencia en el ensayo de tracción”. Pero al hacer una evaluación de varios materiales, podría darse el caso que para cada material escogiese un criterio especifico. Esto forma parte del funcionamiento intrínseco del programa y no lo podemos modificar o saber, sólo intuir.

Lo que podemos modificar es el objetivo. La herramienta nos da a elegir entre dos objetivos a conseguir:

  • Mínima masa.

  • Máxima rigidez.

Y un coeficiente de seguridad que por defecto es 2.

2- Caso de estudio 1: Generación del brazo del chasis

2.1.Generalidades dron acrobático

Los drones acrobáticos, como cualquier UA (aeronave no tripulada) están sujetos a la normativa de REGLAMENTO DE EJECUCIÓN (UE) 2019/947 DE LA COMISIÓN relativo a las normas y los procedimientos aplicables a la utilización de aeronaves no tripuladas. En esta norma se especifican 7 clases de drones en función de sus capacidades operativas y una serie de atributos, peso máximo de despegue (MTOW), geoconsciencia, energía de impacto, entre otras. La clase C0 será la menos restrictiva pero la masa máxima será de 250g. Los modelos comercializados antes de enero de 2021 no están sujetos a esa norma. Si se excede de esa masa, pertenecen a la clase C1 con un peso máximo de 900g. En el presente proyecto partiremos de las especificaciones de un dron acrobático comercializado en 2019 y montado por piezas con un peso de 360 g sin batería. Nótese que no incluimos la batería en nuestros cálculos porque es intercambiable, son del orden de 200 g

El conjunto se compone de:

Denominación Repetición componente Masa [g] Masa total [g] %masa total
Chasis 1 180 180 0,50
Propulsor 4 27 108 0,30
Palas 4 4 16 0,04
Controladora de vuelo 1 10 10 0,03
Receptor emisora 1 8 8 0,02
Emisora de video 1 8 8 0,02
Controlador electrónico de velocidad 4 8 32 0,09
Masa total 362

Dron acrobático montado por piezas.

El chasis de este modelo sigue la tipología estructural más extendida que mezcla laminado compuesto de fibra de carbono y epoxi de alta tenacidad con columnas de aluminio. Los cuatro brazos de construyen como una pieza solidaria de material compuesto a la que se atornilla un rigidizador en forma de cruz, del mismo material, en cada lado. El cuerpo principal actúa como un panel sándwich, donde la tapa o la base hacen de extremos y las columnas, construidas en aluminio y de sección circular, con uniones hiperestáticas a ambos lados hacen de núcleo. El resultado es un chasis muy ligero que a lo largo de su vida útil puede caer varias veces desde alturas considerables sin resultar dañado en exceso. A continuación, se muestran las vistas de este chasis.

Diagrama, Dibujo de ingeniería Descripción generada automáticamente

2.2.Definición geometrías

En el módulo de Generative Design de Fusion 360 abrimos un nuevo estudio estructural. A través de la herramienta de dibujo creamos el boceto de las piezas a unir. En nuestro caso es medio cuerpo del dron y un tubo de pared gruesa que hará de sostén del motor. El ESC, o controlador electrónico de velocidad, debe situarse lo más cerca posible del motor para reducir las pérdidas de disipación eléctrica. Es por ello, que lo definimos como geometría a evitar y lo situamos en el brazo del dron. En el plano 1 Datos geométricos de partida representamos esta figura.

Definición de base de la recurrencia

A continuación, añadimos un cilindro 124,5 mm al soporte del motor como geometría a evitar. Con esto evitaremos que se generen partes de la estructura que interfieran con el paso de las aspas.

2.3.Definición de solicitaciones.

Justificación de secciones del brazo 1

Justificación de secciones del brazo 2

Si situamos un eje coordenado según las reglas de la mano derecha en el centro de gravedad de la aeronave, suponiendo que el avance longitudinal es en el eje X , el giro de la máquina en los ejes X, Y, Z se les denomina alabeo cabeceo y guiñada respectivamente. Las solicitaciones que tendrá que aguantar la estructura del chasis se deben a las reacciones que los propulsores generan al actuar, los momentos que estas fuerzas generan.

Supóngase un eje coordenado local del brazo. Como los propulsores actúan en la dirección vertical, resulta poco intuitivo comprender por qué el eje de máximo momento de inercia de la sección del brazo no es vertical sino horizontal. Esto es así porque a medida que el brazo flexiona, al estar el motor en la punta, cambia la dirección de aplicación de la fuerza, como los brazos tienen una disposición simétrica las componentes horizontales de las fuerzas se anulan entre sí. Conformando un sistema dinámico estable. Además, la controladora de vuelo evita una diferencia de par en los propulsores muy grande para conservar la estabilidad en vuelo. Por otro lado, el alabeo en la dirección del eje del brazo resulta en un sistema de retroalimentación positiva o inestable.

En el método del diseño generativo las secciones de máximo momento de inercia se regulan a través de las solicitaciones de la distintas acciones. En nuestro caso serán:

Acción Fx [N] Fy [N] Fz [N] Mz [N m]
1 Censurado Censurado Censurado Censurado
2 Censurado Censurado Censurado Censurado
3 Censurado Censurado Censurado Censurado
4 Censurado Censurado Censurado Censurado
5 Censurado Censurado Censurado Censurado

Siguiendo las indicaciones del manual de Autodesk, es aconsejable incluir el peso propio en todas las acciones y la diferencia en el orden de magnitud en entre los momentos y las fuerzas deben ser del orden de las indicadas aquí. Por esto las unidades de momento por defecto en el módulo son de Nmm.

El punto de aplicación será la superficie interior del tubo soporte motor. Como restricciones impondremos que el bloqueo en todos los grados de libertad del cuerpo.

2.4.Definición de objetivos y materiales

Los objetivos del criterio de diseño son minimizar la masa de la estructura con un factor mínimo de seguridad de 2. Habilitamos la obtención de formas compatible con fresado de 5 ejes y procesos de fabricación aditiva.

Material Módulo de Young [GPa] Módulo de Poisson Tensión de fluencia [MPa] Tipo
Aluminio AlSI10Mg 71 0,33 240 metal
Rilsan Invent Natural - PA 11 1,5 0,44 43 plástico
Resina Termoplástica 3,3 0,36 57,2 plástico
Polietileno de alta densidad 0,91 0,39 20,7 plástico
Resina fenólica 6,89 0,34 51,7 plástico
Resina Acetal 2,9 0,39 68,2 plástico
Plástico ABS 2,24 0,38 20 plástico

Tabla materiales estudio 1

La elección de los materiales del estudio responde al espíritu exploratorio de este proyecto, el programa solo nos deja elegir 7 materiales por cada estudio. Se seleccionaron los de menor densidad disponibles en la librería del programa, en su mayoría materiales poliméricos. En los capítulos 9 y 12 Nótese la diferencia en orden de magnitud del modulo de Young y de la tensión de fluencia en el ensayo de tracción. En el capitulo 4.6 se analizará la repercusión de estos parámetros en las geometrías obtenidas.

Presionamos generar y se manda el problema a internet, al tiempo se pueden explorar las geometrías solución.

2.5.Resultados

El programa genera cuatro geometrías solución para cada material. Estas se corresponden con un sistema de fabricación cada una y están ordenadas de izquierda a derecha de la siguiente manera:

  • Fabricación aditiva con eje x de cambio de capa

  • Fabricación aditiva con eje y de cambio de capa

  • Fabricación aditiva con eje z de cambio de capa

  • Fabricación mediante fresado de 5 ejes

Nótese que la geometría se encuentra en el primer cuadrante según la regla de la mano derecha.

El número de iteraciones necesarias para cada una es variable con un mínimo de 4 iteraciones y un máximo de 33. A continuación presentamos las soluciones ordenadas por material y el tipo de fabricación antes mencionado.

Resultados estudio 1 Aluminio AlSI10Mg

Resultados estudio 1 Rilsan Invent Natural - PA 11

Resultados estudio 1 Resina termoplástica

Resultados estudio 1 Polietileno de alta densidad

Resultados estudio 1 Resina fenólica

Resultados estudio 1 Resina acetal o polioximetilenos

Resultados estudio 1 Plástico ABS

2.6.Discusión de resultados

El sistema nos deja ordenar los resultados en una gráfica, elegimos representar la masa frente a la energía de distorsión.

Resultados de diseño generativo Masa frente a Energía de distorsión

Llegados a este punto conviene matizar el funcionamiento de la recurrencia en aras de entender porque algunas geometrías solución convergen en una pieza aceptable y otras no. En cada paso de la iteración se calculan para todo el dominio de la pieza las tensiones por el método de los elementos finitos. Acto seguido el programa calcula el coeficiente de seguridad de cada subdominio que depende del criterio de plastificación.

Citando el libro de (Blanco Claraco 2016) “Uno de los objetivos principales de la resistencia de materiales es el de "vigilar" que en ningún punto de la estructura se superen unos determinados niveles de tensión predefinidos. No obstante, como se ha explicado, el concepto de tensión en cada punto es complejo porque el vector que la define es diferente en cada uno de los infinitos planos que lo atraviesan. Para salvar esta circunstancia, diversos autores han establecido criterios simplificados mediante los cuales, a partir de las tensiones principales, se puede verificar que en ningún otro plano que pase por el punto se superan los umbrales admisibles”. A todas luces, el criterio que usa el programa es el de Von Mises y Hencky, que puede formularse en función de las tensiones principales:

σequiv =12[(σ1−σ3)2+(σ1−σ2)2+(σ2−σ3)2]\sigma_{\text{equiv }} = \sqrt{\frac{1}{2}\left\lbrack \left( \sigma_{1} - \sigma_{3} \right)^{2} + \left( \sigma_{1} - \sigma_{2} \right)^{2} + \left( \sigma_{2} - \sigma_{3} \right)^{2} \right\rbrack}

O en función de todas las componentes del tensor de tensiones:

σequiv =12[(σx−σy)2+(σy−σz)2+(σz−σx)2+6(τxy2+τyz2+τzx2)]\sigma_{\text{equiv }} = \sqrt{\frac{1}{2}\left\lbrack \left( \sigma_{x} - \sigma_{y} \right)^{2} + \left( \sigma_{y} - \sigma_{z} \right)^{2} + \left( \sigma_{z} - \sigma_{x} \right)^{2} + 6\left( \tau_{xy}^{2} + \tau_{yz}^{2} + \tau_{zx}^{2} \right) \right\rbrack}

Como el coeficiente de seguridad máximo admisible es dos, en cada iteración el algoritmo ejecuta para cada subdominio:

si: σequiv ×2<σfluencia \sigma_{\text{equiv }} \times 2 < \ \sigma_{\text{fluencia }} elimina material

si: σequiv ×2=σfluencia \sigma_{\text{equiv }} \times 2 = \ \sigma_{\text{fluencia }} no hace nada

si: σequiv ×2>σfluencia \sigma_{\text{equiv }} \times 2 > \ \sigma_{\text{fluencia }} genera material

Si comparamos las soluciones de manera visual y establecemos un criterio subjetivo muy básico:

  • Bueno

  • Malo

  • Normal

Establecemos que las soluciones para los materiales estudiados son:

  • Buenas: Aluminio AlSI10Mg.

  • Malas: Polietileno de alta densidad, Plástico ABS.

  • Normales: el resto.

Quizás el lector se esté preguntando, qué sentido tiene comparar en el mismo estudio, para las mismas solicitaciones e idénticas variables geométricas materiales tan dispares como Aluminio o Polietileno de alta densidad. La respuesta es simple, no estamos comparando los materiales en sí, sino cómo afectan los parámetros de los materiales a la evolución de la recurrencia. Y la pregunta que es necesario hacerse es, ¿Por qué algunos materiales presentan una evolución divergente? En la tabla de materiales del estudio 1 observamos que, para módulo de Young y ratio de Poisson similares, una tensión de fluencia baja será determinante para que la geometría solución se vea sobrepasada por las solicitaciones a lo largo de las iteraciones finales.

En capítulos posteriores se demostrará que las diferencias acusadas en las propiedades de los materiales del estudio puede aprovecharse como un sistema de calibración de las entradas del algoritmo de generación de geometría. De hecho, en los capítulos finales optimizaremos los parámetros del algoritmo para aleaciones y para materiales poliméricos.

Los resultados más prometedores por tener menos masa son los resultados 24-Acetal, 8-Rislan, 12-Termopástico, 4-aluminio

Si la geometría de aluminio es esbelta y aloja de manera satisfactoria la ESC, la de menor peso y menor energía de deformación es 24-Acetal

A continuación, presentamos un gráfico que representa los desplazamientos máximos frente a la masa de la pieza para tener una medida de la rigidez de las soluciones obtenidas.

Resultados de diseño generativo Masa frente a desplazamiento máximo

Como era de suponer, mayores módulos de Young acarreará menos desplazamiento máximo.

A priori, dado que el propulsor necesitará espacio debajo de él para que pase el aire, las geometrías solución deben ser más esbeltas.

La herramienta de exploración de resultados nos permite observar cada solución como una sucesión de instantáneas de todas sus iteraciones, gracias a esto podemos intuir cómo funciona el programa a nivel interno y también detectar geometría frontera que no estén funcionando correctamente. En nuestro caso el cilindro que representa la pala es rellenado por encima de esta con material. Por esto, se descartan las soluciones y reformulamos las geometrías base.

Iteración primera de todos los materiales del estudio 1

En la foto de la primera iteración, que por ser la geometría base común para todos los materiales, es igual para todos ellos. Se puede apreciar que se genera material por encima del rotor. Aunque esa zona está poco solicitada y acabará desapareciendo, se desperdicia esfuerzo computacional como el número de iteraciones posibles tiene un máximo, no permite que se alcancen geometrías óptimas.

Bibliografía

BELTRAN, F., 1999. Teoria General del Metodo de los Elementos Finitos. S.l.: s.n.

BLANCO CLARACO, J.L., 2016. Resistencia de materiales: resumen de teoría y problemas resueltos. Almería: Editorial Universidad de Almería. ISBN 978-84-16642-36-6.

RODRÍGUEZ-AVIAL LLARDENT, M., 2013. Elasticidad y resistencia de materiales II. Madrid: UNED - Universidad Nacional de Educación a Distancia. ISBN 978-84-362-6287-2.

TEDESCHI, A. y ANDREANI, S., 2014. AAD, Algorithms-aided Design: Parametric Strategies Using Grasshopper. S.l.: Le Penseur Publisher. ISBN 978-88-95315-30-0.

ZIENKIEWICZ, O., TAYLOR, R. y ZHU, J.Z., 2013. The Finite Element Method: its Basis and Fundamentals: Seventh Edition. 7th ed. Waltham, MA, Kidlington, Oxford, UK: Butterworth-Heinemann. ISBN 978-1-85617-633-0.