Diferencia entre revisiones de «Optimización de una placa plana»

Revisión del 09:39 20 sep 2014

En este artículo se presenta el enunciado y posterior metodología para la resolución de un problema de optimización.

1 Enunciado y características del problema

Placa plana con agujero objeto del problema

Se define una placa metálica rectangular de pequeño espesor dispuesta en un plano vertical fijada en su borde inferior pudiendo deformarse libremente el resto de los puntos y sometida a cargas en su borde superior. El objeto del problema es optimizar la forma de un agujero con forma elíptica en el interior de la placa sometida a tensión plana de forma que manteniendo constante el área del agujero se alcancen las menores tensiones en la misma.

Los parámetros respecto a los cuales se realiza la optimización son el desplazamiento horizontal del agujero respecto del centro de la placa, el desplazamiento vertical del agujero respecto del centro también y la excentricidad de la elipse. Para resolver el problema se diseña un programa en lenguaje Matlab que de forma automática a partir de la carga actuante en el borde superior y de una configuración inicial cualquiera itera para obtener los valores de los tres parámetros que permiten obtener las tensiones mínimas. El cálculo de las tensiones existentes en la placa se realiza mediante el método de los elementos finitos empleando elementos triangulares CST. La tensión que se adopta como objetivo para minimizar es el máximo de las tensiones de Von Mises que se producen en toda la placa para la hipótesis de carga considerada.

2 Casos

3 Metodología

Se describe a continuación un resumen del proceso seguido para llegar a la resolución final del problema mediante la programación en lenguaje Matlab así como los distintos programas utilizados previamente para intentar calcular la solución.

Inicialmente se decidió resolver el problema mediante el empleo de software libre, para ello se comenzó a plantear el problema con MAT-fem. Este programa se basa en el empleo del programa GID para realizar el preprocesado (mallado, cargas, características de los elementos) y postprocesado (representación de deformaciones, tensiones, desplazamientos) del modelo, resolviendo el sistema lineal de ecuaciones del modelo con Octave. Independientemente de las ventajas de este programa, la optimización por medio de un programa que de forma automática iterase para obtener la solución óptima requiere poder programar con comandos todas y cada uno de los pasos del proceso, lo cual parecía complicado en este caso al requerir guardar los datos del mallado en un archivo en un formato determinado para después poder resolverlo en Octave. Una gran dificultad radicaba en poder ejecutar todos los procesos en el programa GID mediante comandos para poder automatizar todo el programa.

Posteriormente se planteó el problema en el software libre FreeFem++, que es un software orientado a la resolución de ecuaciones diferenciales. La definición de las características del problema, como con GID, era sencillo; sin embargo el principal problema que se encontró en FreeFem++ es que era difícil exportar los resultados para su tratamiento posterior una vez resuelto el problema (obtención de deformaciones, tensiones a partir de desplazamientos y representación de los mismos de forma gráfica). El lenguaje del programa es C++ aunque la ejecución de comandos no funcionaba adecuadamente en este lenguaje.

Debido a los problemas que se encontraron con los programas utilizados y el tiempo que se empleó para intentar solucionarlos, se opta por plantear la resolución del problema mediante el empleo del programa Matlab. Matlab tiene un módulo específico para la resolución de ecuaciones diferenciales, cuyos comandos se utilizaron para resolver el problema particularizados para el caso de tensión plana.

La resolución en Matlab comenzó calculando la tensión de Von Mises máxima que se genera en la placa para una determinada configuración del agujero, excentricidad y carga aplicada con el objetivo de aprender los comandos básicos necesarios a utilizar. A partir de este cálculo inicial se crea la función que para una configuración dada (desplazamiento horizontal y vertical y excentricidad de la elipse) con las condiciones de contorno del problema permite obtener la tensión de Von Mises máxima en la placa; función que se utiliza para el proceso posterior de optimización.

El primer paso consistió en la resolución de un problema concreto con unas determinadas condiciones de contorno y características geométricas en la placa para poder resolver posteriormente un problema con cualquier condición inicial. La resolución implica la obtención del campo vectorial de desplazamientos de cada nodo (dos componentes al ser una placa plana) y a partir de éstos calcular las tensiones normales existentes en la placa para ese estado según el eje X (horizontal) y eje Y (vertical); y las tensiones tangenciales. Una vez obtenidos estos campos escalares es posible calcular en cada punto de la placa la tensión de Von Mises existente.

Entre las mayores dificultades encontradas en el proceso se encuentran la Geometry Description Matrix y la matriz de condiciones de contorno. Los comandos para la creación por pasos de la matrices necesarias para la generación de la Geometry Description Matrix no funcionaban adecuadamente por lo que se opta por escribir la matriz directamente deducida la estructura de un caso similar al del problema, placa rectangular con agujero elíptico, creado en Matlab con ese único objeto. También es necesario definir las condiciones de contorno de cada borde de la región sobre la que se resuelve en una matriz de condiciones de contorno. A pesar de haber acabado comprendiendo la metodología correcta de creación de la matriz, expresar para cada columna las condiciones Neumann o Dirichlet existentes, la utilización de la misma en los cálculos se hace a través de un fichero que es cargado durante la ejecución que contiene dicha matriz generada mediante el módulo de Matlab para la resolución de ecuaciones diferenciales.

%Programa que realiza la optimizacion de 
%la posicion y forma del agujero elíptico
%de una placa sometida a tensión plana
clear all
clc
%comenzamos con una terna cualquiera
x=0.1; y=0.1; e=0.2;
%tamaño de los triángulos
h=0.03;
%paso de la derivada
hg=0.01;
epsilon=0.05;
a0=[x,y,e];
Fa0=VM_max(a0(1),a0(2),a0(3),h);
Fdib=[Fa0];
Fa1=Fa0+1;
gradiente0=[0,0,0];
gradiente0(1)=(VM_max(a0(1)+hg,a0(2),a0(3),h)-Fa0)/hg;
gradiente0(2)=(VM_max(a0(1),a0(2)+hg,a0(3),h)-Fa0)/hg;
gradiente0(3)=(VM_max(a0(1),a0(2),a0(3)+hg,h)-Fa0)/hg;
while(norm(gradiente0)>1e-6)
%calculamos los elementos del gradiente

a1=a0-gradiente0*epsilon;
Fa1=VM_max(a1(1),a1(2),a1(3),h);
if(Fa1<Fa0)
    epsilon=epsilon*1.1;
    a0=a1
    Fa0=Fa1
    Fdib=[Fdib,Fa0];
    gradiente0(1)=(VM_max(a0(1)+hg,a0(2),a0(3),h)-Fa0)/hg;
    gradiente0(2)=(VM_max(a0(1),a0(2)+hg,a0(3),h)-Fa0)/hg;
    gradiente0(3)=(VM_max(a0(1),a0(2),a0(3)+hg,h)-Fa0)/hg;
else
    epsilon=epsilon*0.5;
end
if epsilon<1.e-8
    break
end

end

function [ vm ] = VM_max( x,y,e,h)

%las características del rectángulo deben de ser las 
% mismas en la funcion y en el fichero principal original
B=1; H=2; %base y altura
% centrada la elipse en el origen x0=0; y0=0; 
xmin=-0.5*B; xmax=B*0.5; ymin=-H*0.5; ymax=H*0.5;

%obtenemos ahora los valores de a y b de la 
%elipse a partir de la excentricidad y sabiendo
% que mantiene el mismo área que la orignal
a_original=0.3; b_original=0.25;
area_elipse=3.1415*a_original*b_original;
%resolvemos el sistema conocida la excentricidad
% y el área
b=((1-e^2)*area_elipse^2/3.1415)^(1/4);
a=area_elipse/(3.1415*b);

%comenzamos a crear el rectángulo
dl=[2 2 2 2 4 4 4 4;xmax xmax xmin xmin -a 0 a 0;
    xmax xmin xmin xmax 0 a 0 -a;ymin ymax ymin ymin 0 -b 0 b;
    ymax ymax ymax ymin -b 0 b 0;ymax ymax 0 ymax 0 0 0 0;
    0 0 ymax 0 ymax ymax ymax ymax;0 0 0 0 0 0 0 0;
    0 0 0 0 0 0 0 0;0 0 0 0 a a a a;
    0 0 0 0 b b b b;0 0 0 0 0 0 0 0];
%mallamos la geometría
[p,e,t]=initmesh(dl);
%damos ahora el desplazamiento al agujero según queramos
for k=1:length(e)
        if 4<e(5,k) &&  e(5,k)<9 %los que están en la frontera interior
        j=e(1,k);
        
        p(1,j)=p(1,j)+x;%las x
        p(2,j)=p(2,j)+y;%las y
        end
end
%arreglamos la malla
p1=jigglemesh(p,e,t);

%cargamos el fichero con las condiciones de contorno
load matrizb
% Creamos el vector c para resolver el problema de la
% elasticidad en tensión plana
% características de la placa
E= 200e9; nu= 0.2;
c=zeros(10,1);
c(1)=E/(1-nu^2);
c(2)=0;
c(3)=(E*(1-nu))/(2*(1-nu^2));
c(4)=0;
c(5)=(E*(1-nu))/(2*(1-nu^2));
c(6)=E*nu/(1-nu^2);
c(7)=0;
c(8)=(E*(1-nu))/(2*(1-nu^2));
c(9)=0;
c(10)=E/(1-nu^2);
a=0;
f=[0;0];
% resolvemos
u = assempde(b, p1, e, t, c, a, f);

[ux,uy]=pdegrad(p,t,u);
% deformaciones horizontales
epsilonx=ux(1,:);
%deformaciones verticales
epsilony=uy(2,:);
% deformaciones angulares
gammaxy=ux(2,:)+uy(1,:);
%calculo del vector de tensiones en la placa
sigmax=(E/(1-nu^2))*(epsilonx+nu*epsilony);
%
sigmay=(E/(1-nu^2))*(nu*epsilonx+epsilony);
%
tauxy=(E/(1-nu^2))*(1-nu)*0.5*gammaxy;
%calculo de la tension de von mises
von_mises=sqrt(sigmax.^2-sigmax.*sigmay+sigmay.^2+3*tauxy.^2);
%obtenemos el maáximo de la tension de von mises
vm=max(von_mises);



end

--Gonzalo (discusión) 11:28 29 ago 2014 (CEST)