Nuevo planeta

Hemos creado un nuevo planeta de blogs de estructuras:

www.planetaki.com/soloestructuras

Desde ahí podréis leer los últimos artículos que subamos desde los diferentes blogs que se vayan agregando.
De momento hay pocos blogs suscritos. Todo se andará.
Si tenéis un blog para agregar, lo podéis sugerir como comentario a este post.

Introducción a Octave

Octave es un lenguaje de programación de alto nivel destinado al análisis numérico.

Apareció en 1988 como apoyo para un libro de texto sobre diseño de reactores químicos, pero en los años 90 se extendió a otros ámbitos, superando sus limitaciones iniciales y convirtiéndose en una herramienta potente, de fácil uso y  rápido aprendizaje

Se trata de un lenguaje interpretado, con las ventajas que supone por su facilidad de depuración y la independencia de la máquina, pero también con sus inconvenientes,en especial su lentitud. En ocasiones, se emplean complementos en lenguaje compilado, generalmente C o C++ (el propio octave está escrito en C++), que permiten acelerar procesos críticos. No es adecuado para grandes proyectos, ya que resulta lento en su ejecución, no permite una buena organización de datos y algoritmos  complejos y su capacidad gráfica es muy limitada. Su mayor virtud reside en su enorme biblioteca de funciones matemáticas.

Octave trabaja en modo consola y permite el uso de scripts, para lo que incorpora un sencillo editor.  Para facilitar la entrada de datos se pueden utilizar herramientas externas, como qtoctave.

En esta página de la Oficina de Software Libre de la Universidad de Zaragoza se pueden encontrar una síntesis del programa, enlaces de descarga del propio programa y de varios complementos.

El programa y El manual de referencia se encuentra en la propia web del programa. Existe también una wiki, en inglés, que cuenta con un apartado en español. En español se pueden encontrar en la red manuales muy breves o algo más completos

El objetivo de este blog no es ofrecer otro manual de octave, sino descubrir poco a poco sus aplicaciones en cálculo de estructuras.

Tensor de deformaciones. Diagonalizar matrices con Maxima

Empezamos la colección de ejercicios resueltos con este problema que aparece en el libro "Resistencia de materiales" de Luis Ortiz Berrocal, en la página 51, ejercicio 1.6:

Enunciado

En un determinado punto P de la superficie de un sólido elástico sometido a carga se midieron las siguientes deformaciones: Un alargamiento longitudinal unitario de 0.0001 en la dirección x; un acortamiento longitudinal unitario de 0.0005 en la dirección de y, perpendicular a x y una deformación angular de -2·Raíz(7)·10^-4 rad.

1. Calcular analítica y gráficamente las deformaciones principales, así como las direcciones correspondientes.
2. Determinar el valor de la deformación angular máxima.
3. Hallar las componentes intrínsecas del vector deformación unitaria correspondiente a una dirección contenida en el plano xy que forma un ángulo de 45º con la dirección positiva del eje x, contado en sentido antihorario
   

Usaremos

• matrix([fila1],[fila2]...): Para definir una variable como matriz. El argumento se compone por las filas, pasadas como listas, separadas por comas. Y cada fila, que es una lista, se limita por los corchetes con los elementos separados por comas. Es más difícil contarlo que hacerlo...
• diagmatrix(rango,valor): Genera una matriz diagonal del rango que definamos y con el valor en la diagonal que definamos
• sqrt(num): Raíz cuadrada de un número
• determinant(matriz): Calcula el determinante de la matriz
• .: El punto es el operador para el producto de dos matrices. (El asterisco * es el operador del producto miembro a miembro.
• Matriz[fila,columna]: Para referirnos al término de la matriz de las filas y columnas indicadas.
  

Resolución

No vamos a hacer la resolución gráfica (círculo de Mohr) porque el ejercicio es un pretexto para ver las matrices en Maxima. Buscaremos las deformaciones principales y las direcciones principales.
Escribimos la matriz. Como todo está multiplicado por 10^-4 lo podemos obviar: Los resultados estarán también multiplicados por 10^-4
La matriz M:

M:matrix([1,-sqrt(7)],[-sqrt(7),-5]);

Al resultado anterior (%) le restamos una matriz diagonal (de rango 2 y de valor e):

%-diagmatrix(2,e);

Calculamos el determinante;

determinant(%);

Y resolvemos la ecuación de el determinante igual a cero

solve(%=0, e);

La dos soluciones de e son las deformaciones principales (recordemos que está multiplicado por 10^-4).

Las direcciones de los vectores principales las dejo para otro día se calculan sustituyendo las soluciones en la matriz de (%o2) y multiplicando por el vector de coordenadas incógnita. Hay dos soluciones que nos definen dos matrices:

M1:subst(e=2,%o2);
M2:subst(e=-6,%o2);

Definimos el vector (columna):

V:matrix([x],[y]);

Y las matrices resultantes de multiplicar M1 y M2 por el vector V.

M1V:M1.V;
M2V:M2.V;

Las coordenadas del vector resultante deben ser cero. Resolvemos el sistema (los dos sistemas)

solve([M1V[1,1]=0,M1V[2,1]=0], [x,y]);
solve([M2V[1,1]=0,M2V[2,1]=0], [x,y]);

Y el resultado es paramétrico (los parámetros son %r1 y %r2). Claro, buscamos una dirección, no un vector. Cualquier valor de los parámetros nos dará un vector de cualquier módulo, pero siempre la misma dirección)

Seguro que el último paso se puede hacer más sencillo, pero tengo que saber más Maxima

Viga de un solo vano: Script en Máxima

En la viga de un solo vano empezamos a tratar Máxima como una calculadora. Cada vez que nos daba una respuesta, le pedíamos que trabajase con ella. Ahora, lo que vamos a intentar es modificar lo que hicimos para que sirva para casi cualquier cosa. Es decir, hacer un script más genérico...

Usaremos

Lo que ya aprendimos en la viga de un solo vano: Integrate, subst, solve...

• first(A): Para acceder al primer elemento de la lista A (igual funcionan second(A), third(A)...)
  
• part(A,i): Para acceder al elemento en posición i de la lista A
• rhs: accede a lo que hay a la derecha del igual en una ecuación
• map(función, lista): Aplica una determinada función a cada elemento de la lista
  

Para mayor claridad (y para tener que subir menos imágenes) vamos a enviar las órdenes en grupos.

Integrando

Esta parte es casi igual que en el artículo anterior, pero asignamos e integramos en una sola orden:

v:integrate(q,x)+v0;
m:integrate(v,x)+m0;
c:m/EI;
g:integrate(c,x)+g0;
f:integrate(g,x)+f0;

Podríamos modificar q, y en vez de ser un valor constante, podría ser cualquier función continua.

Ecuaciones de borde

Aquí si hay modificaciones y aclaraciones:
La respuesta de solve() es una lista con todas las soluciones posibles. Es decir, si el sistema no fuera lineal habría más de una solución posible. Cada solución sería un elemento de la lista.
Y si pedimos que nos resuelva más de una incógnita (como es el caso) entonces cada elemento de la lista de soluciones sería, a su vez, una lista de incógnitas resueltas. Solve nos da una lista de listas...
Sabemos que la lista de soluciones sólo tendrá una solución (es un sistema lineal), así que sólo nos interesa el primer elemento: Asignamos el primer elemento (first) a la variable soluciones:

/*Viga biempotrada: los apoyos ni bajan ni giran: 4 ecuaciones, 4 incógnitas */
ec1:subst(x=0,f)=0;
ec2:subst(x=l,f)=0;
ec3:subst(x=0,g)=0;
ec4:subst(x=l,g)=0;
solve([ec1,ec2,ec3,ec4], [v0,m0,g0,f0]);
soluciones:first(%);

Definición de las funciones de la viga

Una vez asignadas los distintos valores de las distintas incógnitas a la lista "soluciones", las metemos en las funciones con la orden subst. Comento la primera línea, el resto son similares:
Definimos la funcion v sustituyendo (subst) v0 por el valor que hemos obtenido:
De la lista "soluciones" necesitamos el primer valor y lo cogemos con part(soluciones,1).
Todos los valores de la lista "soluciones" son igualdades. Lo que queremos es lo que está a la derecha del signo "=" y lo obtenemos con la función rhs().
En el resto de sustituciones, hay que sustituir más de un valor, por lo que metemos la lista de valores a sustituir entre corchetes (por que es una lista)

/*Conocidas las constantes de integracion, escribo las ecuaciones completas */
"Cortante";
v:subst(v0=rhs(part(soluciones,1)),v);
"Momento";
m:subst([v0=rhs(part(soluciones,1)), m0=rhs(part(soluciones,2))],m);
"Giro";
g:subst([v0=rhs(part(soluciones,1)), m0=rhs(part(soluciones,2)), g0=rhs(part(soluciones,3))],g);
"Flecha";
f:subst([v0=rhs(part(soluciones,1)), m0=rhs(part(soluciones,2)), g0=rhs(part(soluciones,3)), f0=rhs(part(soluciones,4))],f);

Análisis de las funciones

Una vez definidas las funciones, podemos buscar los ceros, los valores en los extremos y los valores máximos. Como cada función es la integral de la anterior, los ceros de unas nos indicarán los máximos de otras.
Buscamos las x que hacen que la función sea cero con solve. y las guardamos en una lista: xv0 sería la lista de valores de x que hacen que el cortante sea cero.
Ahí donde el cortante es cero, el momento será máximo.
Y para localizar los valores máximos usamos la función map, que aplicará la función sobre todos los elementos de la lista:
La función a aplicar lambda la definimos con la lista de parámetros que necesita ([i]) y como se aplica a la finción.
Al estar anidadas es un poco más confuso de ver, pero si nos fijamos donde se abren y se cierran los paréntesis, es muy aclaratorio.
Definimos listas con nombre .ext para los valores en los extremos, y listas con nombre .max para valores máximos.

/*busco los ceros y miro los máximos y los valores en los extremos*/
"Cortante extremos";
vext:map(lambda([i], subst(x=i,v)),[0,l]);
"Cortante=0 en"; xv0:solve(v=0,x);
"Momento max, Momento extremos";
mmax:map(lambda([i], subst(x=rhs(i),m)),xv0);
mext:map(lambda([i], subst(x=i,m)),[0,l]);
"Momento=0 en"; xm0:solve(m=0,x);
"Giro max, giro extremo";
gmax:map(lambda([i], subst(x=rhs(i),g)),xm0);
gext:map(lambda([i], subst(x=i,g)),[0,l]);
"Giro=0 en"; xg0:solve(g=0,x);
"Flecha max, flecha en extremos";
fmax:map(lambda([i], subst(x=rhs(i),f)),xg0);
fext:map(lambda([i], subst(x=i,f)),[0,l]);
"Flecha=0 en"; xf0:solve(f=0,x);

(Pulsa en la imagen, porque si no, te quedarás miope)

Varios finales alternativos

¿Para que vale el script?
Para adaptarlo a los diferentes casos que se nos ocurran:

1. Podemos modificar la ley de cargas, y estudiar casos en los que la carga no es constante.
2. Podemos modificar las condiciones de borde, y estudiar sin mucho trabajo otros casos típicos o no tan típicos
3. Podemos convertir algunas condiciones de borde en compatibilidad de deformaciones y equilibrio.

Lo más difícil ya está hecho. Ahora sólo se trata de jugar con ello... (otro día)

Viga de un solo vano

Vamos a resolver una viga de un vano con maxima:

Las ecuaciones de equilibrio

Las ecuaciones de equilibrio en una viga, si la carga es una ley continua, son dos integrales:

1. El cortante (v) tiene que equilibrar la suma de las fuerzas que actúan sobre la viga: Integrando la ley de cargas obtenemos la ley de cortantes
2. De la misma forma, integrando la ley de cortantes, obtenemos la ley de momentos (m).
   

Recordemos que la integral de una función no es una función inequívoca, si no una familia de funciones. O dicho de otra manera, aparece la constante de integración, que es el valor de la integral en el origen (v0, m0).

Las condiciones del material

La relación tensión deformación en una viga a flexión se traduce en una ecuación muy sencilla:

c=m/EI

Siendo c la curvatura (el inverso del radio de curvatura), I la inercia del perfil y E el módulo de Young del material.

El problema geométrico

Ahora todo se reduce a geometría: Si tenemos la curvatura en todos los puntos de la viga, deberíamos poder deducir la flecha en cada punto:

1. El giro (g) (o pendiente) es la integral de la curvatura
2. La flecha (f) es la integral de la pendiente

Las condiciones de borde

Cada vez que integramos generamos una constante de integración, que en principio es una incógnita. Tenemos por tanto 4 incógnitas, para las que necesitamos 4 ecuaciones. Nos las darán las condiciones de apoyo de una forma muy sencilla (Cada extremo, dos ecuaciones):

1. ¿Está impedido el descenso?
   
   1. Sí: El descenso en ese punto vale cero: f=0
      
   2. No: El cortante en ese punto vale cero: v=0
   
   
2. ¿Está impedido el giro?
1. Sí: El giro en ese punto vale cero: g=0
   
2. No: El momento en ese punto vale cero: m=0

Breve introducción a Maxima

Maxima es un programa que nació en los 70. Tiene una interfaz gráfica que se llama wxMaxima.
Este blog no pretende ser un manual exhaustivo (que ya está escrito), ni uno de primeros pasos (que también). Explicaré todas las funciones que use, pero recomiendo que acudáis al manual de Mario Rodríguez Riotorto, que es con el que estoy aprendiendo yo...

Usaremos:

• ; :Punto y coma para terminar la orden (y esperar respuyesta)
• : :Dos puntos para asignar un valor (o función) a una variable.
• %: Para referirnos a la última respuesta
  
• Integrate(A,B): Integra la función A respecto de B
• Subst([A1=B1, A2=B2...], C): Sustituye por B todos los valores A que encuentre en la expresion C
• Solve([lista de acuaciones], [lista de incógnitas]): Resuelve un sistema de acuaciones
• /*...*/: Para añadir comentarios

Con CTRL+INTRO ejecutamos el comando

Empezamos

Empezamos con una ley de carga constante, q, aunque podríamos haber empezado con cualquier ley continua. La carga repartida constante es, probablemente, la más habitual.
El principio sólo hay que integrar:

(%i1) v:integrate(q, x)+v0;

(%i2) m:integrate(%, x)+m0;

(%i3) c:%/EI;

(%i4) g:integrate(%, x)+g0;

(%i5) f:integrate(%, x)+f0;

Llegado este punto, procedemos a resolver el sistema de 4 ecuaciones con 4 incógnitas que nos definirá las constantes de integración. Por mucho que haya términos a la cuarta, es un problema lineal, porque todas las condiciones nos darán un valor de x que debemos sustituir:

(%i6) /*Viga biempotrada: los apoyos ni bajan ni giran: 4 ecuaciones, 4 incógnitas */
ec1:subst(x=0,f)=0;
ec2:subst(x=l,f)=0;
ec3:subst(x=0,g)=0;
ec4:subst(x=l,g)=0;

(%i10) solve([ec1,ec2,ec3,ec4], [v0,m0,g0,f0]);
/*tenemos todas las constantes de integración*/

(%i11) fbiempotrada:subst([v0=-l*q/2, m0=l^2*q/12, g0=0, f0=0],f);
/*Y las utilizamos para la función flecha*/

(%i12)subst(x=l/2,%);
/*La flecha en el centro del vano (x=l/2)*/

Que es lo que ya sabíamos...
Otro día haremos una viga biapoyada, una apoyada-empotrada, un voladizo, o una viga de varios vanos en continuidad...

Presentación del blog

Para no acabar como el Capitán McCrea...

Entre los muchos buenos propósitos para el año que comienza estaban tres que son los que dan origen a este blog:

1. Aprender un lenguaje de programación (he escogido Python)
2. Aprender alguna herramienta de cálculo numérico (he escogido Octave)
3. Aprender una herramienta de cálculo simbólico (he escogido Maxima)

Los tres son programas libres, y Octave y Maxima son además de código abierto y licencia GNU GPL

¿Por qué?

El uso de los ordenadores nos sirve para trabajar con una eficacia mucho mayor.
Pero también nos está limitando el pensamiento: Utilizamos programas, que han sido programados por otros, y que nos eximen de enfrentarnos a problemas, que son precisamente los que pretende resolver el programa.
Si le pedimos al programa algo que se salga un poco de su rango de aplicación, no lo hará correctamente. Es decir, o no lo hace, o lo hace de forma equívoca...
Si alguna vez supimos resolverlo a mano, lo vamos olvidando.
Si nunca supimos, nunca lo aprenderemos. Tenemos que entregar en un plazo corto, que no nos deja que nos recreemos en la investigación y el aprendizaje.
Para no saber cada vez un poquito menos, tengo que ponerme un pequeño reto: Éste que da origen al blog: Aprender herramientas, y obligarme a contar como se usan, aplicado al cálculo de estructuras.

¿Por qué más?

Además, y por los mismos motivos, doy clases de resistencia de materiales.
Sirva este blog para recopilar ejercicios resueltos y explicaciones parciales...

El nombre

Evidentemente hace referencia a uno de los parámetros fundamentales en la resistencia de materiales (Tensión), y a una de las herramientas que quería aprender (Máxima). (No he estado muy ocurrente)