Implementación de una nueva formulación de deformación plana generalizada en un código de elementos finitos en MATLAB

i Equation Chapter 1 Section 1 Trabajo Fin de Grado Grado en Ingeniería de las Tecnologías Industriales Implementación de una nueva formulación de Deformación Plana Generalizada en un código de elementos finitos en MATLAB Autor: Manuel Seco Ruiz Tutor: Vladislav Mantič, José Reinoso Grupo de Elasticidad y Resistencia de Materiales Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2016 ii iii Trabajo Fin de Grado Grado en Ingeniería de las Tecnologías Industriales Implementación de una nueva formulación de Deformación Plana Generalizada en un código de elementos finitos en MATLAB Autor: Manuel Seco Ruiz Tutor: Vladislav Mantič José Reinoso Grupo de Elasticidad y Resistencia de Materiales Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2016 v Trabajo Fin de Grado: Implementación de una nueva formulación de Deformación Plana Generalizada en un código de elementos finitos en MATLAB Autor: Manuel Seco Ruiz Tutor: Antonio Blázquez, Vladislav Mantič El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, 2016 El Secretario del Tribunal vii A mi familia ix Agradecimientos Quiero mostrar mis agradecimientos en primer lugar a la Escuela Superior de Ingeniería de Sevilla (y, por ende, a la Universidad de Sevilla) por, a pesar de todo, haberme dado las herramientas y las posibilidades de estudiar una carrera universitaria, una ingeniería. Doy las gracias también a todas las personas que he conocido en estos años, profesores y compañeros de clase, por todo lo que me han aportado. Pero, sobre todo, quiero dar las gracias a “Esiland”, a la “Shavalería” y a los “Máquinas”, no solo compañeros de clase, sino amigos. No me quiero olvidar de los que me han acompañado en estos cinco años de caminar universitario, mi familia, los amigos del SARUS, mis experiencias erasmus en Suecia y toda la gente de los Equipos. Mencionar de manera especial a mis tutores, Antonio Blázquez y Vladislav Mantič, por haberme guiado en este largo trabajo, y a José Reinoso, por su inestimable ayuda cada vez que la necesité. Extiendo este agradecimiento a todo el Grupo de Elasticidad y Resistencia de Materiales, por todos los recursos que nos han ofrecido a los estudiantes que hemos realizado el Trabajo Fin de Grado con ellos. xvi xvii ÍNDICE DE FIGURAS Figura 1. Descripción de la geometría 5 Figura 2. Elemento serendípito de 8 nodos (●) y 9 puntos de integración (■) en el espacio natural 9 Figura 3. Funciones de forma del elemento representadas en el espacio natural 10 Figura 4.Situación deformada (-) 19 Figura 5. Elemento sometido a ensayo de tracción 19 Figura 6. Configuración deformada e indeformada del problema 6.3 22 Figura 7. Situación deformada del problema 6.9 26 Figura 8. Malla de nodos y elementos del problema 6.9 26 Figura 9. Situación deformada e indeformada del problema 6.10 27 Figura 10. Malla de nodos y elementos 27 Figura 11. Puntos de Gauss 29 Figura 12. Situación deformada e indeformada del problema 6.11 29 Figura 13. Situación deformada e indeformada del problema 6.12 31 xix ÍNDICE DE TABLAS Tabla 1. Coordenadas y pesos de los 9 puntos de Gauss del elemento 9 Tabla 2. Comparación de resultados del ensayo 6.1 20 Tabla 3. Comparación entre resultados analíticos y numéricos del problema 6.2 20 Tabla 4. Comparación entre resultados numéricos y analíticos del problema 6.3 21 Tabla 5. Comparación de resultados del problema 6.4 23 Tabla 6. Comparación de resultados del problema 6.5 23 Tabla 7. Comparación de resultados analíticos y numéricos del problema 6.6 24 Tabla 8. Comparación de resultados analíticos y numéricos del problema 6.7 25 Tabla 9. Comparación de resultados analíticos y numéricos del problema 6.8 25 Tabla 10. Resultados en desplazamientos de los nodos del contorno de los problemas 6.1 y 6.9 26 Tabla 11. Comparación de resultados analíticos y numéricos del problema 6.10 29 Tabla 12. Comparación de resultados del test de la parcela del problema 6.11 31 1 1 INTRODUCCIÓN L trabajo que se presenta consiste en la implementación de la formulación del Problema Plano Generalizado (PPG) en un código de Elementos Finitos. El elemento resultante de este trabajo es un elemento plano con capacidad de sufrir desplazamientos en la dirección normal al plano. Este elemento será útil para estudiar problemas de barras de material compuesto, en los cuales las cargas contenidas en el plano pueden provocar desplazamientos fuera de él debido a la naturaleza anisótropa del material. El código utilizado para la implementación del PPG ha sido escrito en MATLAB ®. La intención del autor es continuar el desarrollo de este elemento en un futuro TFM, traduciendo el código a lenguaje Fortran e implementándolo en una subrutina para Abaqus ®. E 3 2 ESTRUCTURA DEL TRABAJO El trabajo se estructurará en tres bloques. En el primer bloque se realizará una presentación de toda la teoría que hay detrás de este trabajo, es decir, la formulación del Problema Plano Generalizado y la formulación del elemento serendípito de 8 nodos. En el segundo bloque se desarrollará el procedimiento seguido para el cáculo de la matriz de rigidez de un elemento y se explicará el código de ensamblaje para problemas con más de un elemento. El tercer bloque consiste en la validación de resultados. Se comprobará que la matriz de rigidez de un elemento ha sido correctamente calculada sometiendo el elemento a problemas sencillos, y comparando los resultados obtenidos con la solución analítica de dichos problemas. Además, se realizará un “test de la parcela” o “patch test” para comprobar la convergencia del elemento, es decir, que al aumentar el número de elementos la solución obtenida converge a la solución exacta. 5 3 EL PROBLEMA PLANO GENERALIZADO Considérese un cuerpo homogéneo de forma cilíndrica recta y sección transversal constante. Se asume que el cuerpo tiene un comportamiento anisótropo. Sea D el dominio tridimensional completo, ∂D el contorno lateral del cilindro y ∂De las dos secciones extremas del cilindro si este es de longitud finita (e=b o t, ya sea la sección extrema inferior o superior, respectivamente). La sección transversal del cilindro D determina un dominio bidimensional Ω, con contorno ∂Ω determinado por la sección transversal de ∂D (véase Figura 1.). Se considera el caso en que este cuerpo estará sujeto a cargas en el dominio o en el contorno y restricciones en los desplazamientos de tal forma que las soluciones de tensiones y deformaciones dependan únicamente de las coordenadas del plano perpendicular al eje del cilindro. Si el eje longitudinal del cilindro es paralelo al eje 𝑥3, entonces los tensores de tensiones y deformaciones pueden expresarse como: 𝜀𝑖𝑗=𝜀𝑖𝑗(𝑥1,𝑥2) y 𝜎𝑖𝑗=𝜎𝑖𝑗(𝑥1,𝑥2) 𝑖,𝑗=1,2,3 ( 1 ) La situación descrita anteriormente es denominada Problema Plano Generalizado (PPG). Las situaciones de deformación plana, deformación plana generalizada, tensión plana generalizada, así como torsión libre pueden considerarse casos particulares del PPG. 3.1 Formulación del Problema Plano Generalizado (PPG) Se integra la relación 𝜀−𝑢 𝜀𝑖𝑗=12(𝑢𝑖,𝑗+𝑢𝑗,𝑖) ( 2 ) a lo largo de la coordenada 𝑥3, obteniéndose: D ∂Dt ∂Db Ω ∂Ω Figura 1. Descripción de la geometría Pseudocódigo 12 12 5.1.3 Evaluación de las funciones de forma y sus derivadas parciales El primer paso a seguir una vez se saben las coordenadas del punto de integración actual (xi y eta para kip) es evaluar las funciones de forma y sus derivadas parciales en dicho punto, puesto que será necesario utilizarlas a lo largo del código. La función kshapefunctions2 contiene las funciones de forma detalladas en las ecuaciones ( 11 ) y ( 12 ), y devuelve los valores de todas las funciones de forma y sus derivadas parciales en cualquier punto del elemento. Nótese que, de momento, todos estos cálculos se realizan en el espacio natural. %Evalúa las funciones de forma y sus derivadas [shapef,dshape] =kshapefunctions2(xi,eta,ndim,nnode); %-------------------------------------------------- 5.1.4 Cálculo del jacobiano A continuación, se calcula el jacobiano en el punto de integración. El jacobiano (y su inversa) será necesario para resolver las integrales del cálculo de la matriz de rigidez mediante cuadratura de Gauss (para realizar la transformación de coordenadas naturales a coordenadas reales) y para poder aplicar la regla de la cadena a la hora de expresar las derivadas de las funciones de forma en el espacio real. La función kJacobian2D calcula la matriz jacobiana en el punto donde se han evaluado previamente las funciones de forma: %-------------------------------------------------- % Calcula el Jacobiano y su inversa [xjac,xjacinv,detxjac]=kJacobian2D(ndim,nnode,shapef,dshape,nodes _coord); %-------------------------------------------------- El jacobiano es una matriz cuyos componentes son las derivadas de las coordenadas cartesianas respecto a las coordenadas naturales: 𝐽(𝜉1,𝜉2)=𝜕(𝑥,𝑦) 𝜕(𝜉1,𝜉2)= [ 𝜕𝑥 𝜕𝜉1𝜕𝑦 𝜕𝜉1 𝜕𝑥 𝜕𝜉2𝜕𝑦 𝜕𝜉2 ] ( 13 ) Para el caso de los elementos cuadriláteros, el operador de compatibilidad B (como se verá más adelante), no admite una forma algebraica concreta y debe ser evaluado a partir de las funciones de forma. De igual manera se ha de actuar con la matriz jacobiana, ya que esta se empleará en el cálculo de las componentes de la matriz B. Así, haciendo uso del concepto de discretización, los términos del jacobiano pueden expresarse: 𝑥≈∑𝑁𝐴(𝜉1,𝜉2)𝑥𝐴 𝑛𝑛 𝐴=1 ; 𝑦≈∑𝑁𝐴(𝜉1,𝜉2)𝑦𝐴 𝑛𝑛 𝐴=1 ( 14 ) 𝜕𝑥 𝜕𝜉𝑖≈∑𝜕𝑁𝐴(𝜉1,𝜉2) 𝜕𝜉𝑖𝑥𝐴 𝑛𝑛 𝐴=1 ;𝑖=1,2 𝜕𝑦 𝜕𝜉𝑖≈∑𝜕𝑁𝐴(𝜉1,𝜉2) 𝜕𝜉𝑖𝑦𝐴 𝑛𝑛 𝐴=1 ;𝑖=1,2 ( 15 ) Siendo 𝑛𝑛 el número de nodos (8) del elemento, 𝑁𝐴 la función de forma del nodo 𝐴, y (𝑥𝐴,𝑦𝐴) las coordenadas del nodo 𝐴 en el espacio real. Lo expresado en las ecuaciones ( 13 ) y ( 15 ) se implementa en las líneas de 13 13 Implementación de una nueva formulación de Deformación Plana Generalizada en un código de elementos finitos en MATLAB código de la función kJacobian2D que se muestran a continuación: function [xjac,xjacinv,detxjac]=kJacobian2D(ndim,nnode,shapef,dshape,xref) %Fase de inicialización for i=1:ndim for j=1:ndim xjac(i,j)=0.0; xjacinv(i,j)=0.0; end end % xref son las coordenadas de los nodos en el espacio real for i=1:nnode xjac(1,1) = xjac(1,1) + dshape(1,i)*xref(1,i); xjac(1,2) = xjac(1,2) + dshape(1,i)*xref(2,i); xjac(2,1) = xjac(2,1) + dshape(2,i)*xref(1,i); xjac(2,2) = xjac(2,2) + dshape(2,i)*xref(2,i); end %Calcula el determinante del jacobiano detxjac=det(xjac); if(detxjac > 0.0) %Puede formularse la inversa %Calcula la inversa xjacinv = inv(xjac); else end end 5.1.5 Cálculo del operador de compatibilidad B Una vez se han obtenido los valores de las funciones de forma (en coordenadas naturales) y la matriz jacobiana (que permite expresar estas coordenadas naturales en coordenadas reales), se procede al cálculo de la matriz B u operador de compatibilidad. En este caso, la matriz B adopta la siguiente forma: 𝐵= [ 0 0 0 0 0 0 0 0 𝑥 𝑦 1 0 0 0 0 −𝑥 0 0 0 𝑦 0 0 0 0 | | | 𝜕𝑁1 𝜕𝑥 0 0 … 𝜕𝑁8 𝜕𝑥 0 0 0𝜕𝑁1 𝜕𝑦 0 … 0 𝜕𝑁8 𝜕𝑦 0 0 0 0 … 0 0 0 0 0 𝜕𝑁1 𝜕𝑦 … 0 0 𝜕𝑁8 𝜕𝑦 0 0 𝜕𝑁1 𝜕𝑥 … 0 0 𝜕𝑁8 𝜕𝑥 𝜕𝑁1 𝜕𝑦 𝜕𝑁1 𝜕𝑥 0 … 𝜕𝑁8 𝜕𝑦 𝜕𝑁8 𝜕𝑥 0 ] ( 16 ) Siendo (𝑥,𝑦) las coordenadas del punto de Gauss en el que se está evaluando la matriz. Por tanto, hay que transformar las coordenadas (𝜉1,𝜉2) en (𝑥,𝑦) a través de la ecuación ( 14 ). Las derivadas parciales de las funciones de forma en la matriz B son respecto a las coordenadas reales, mientras que las ya conocidas y calculadas en la matriz dshape son las derivadas respecto a las coordenadas naturales. Para realizar la transformación, se aplica la regla de la cadena: Pseudocódigo 14 14 𝜕𝑁𝑖 𝜕𝑥=𝜕𝑁𝑖 𝜕𝜉1𝜕𝜉1 𝜕𝑥+𝜕𝑁𝑖 𝜕𝜉2𝜕𝜉2 𝜕𝑥 𝜕𝑁𝑖 𝜕𝑦=𝜕𝑁𝑖 𝜕𝜉1𝜕𝜉1 𝜕𝑦+𝜕𝑁𝑖 𝜕𝜉2𝜕𝜉2 𝜕𝑦 ( 17 ) Los términos 𝜕𝑁𝑖 𝜕𝜉𝑗 están almacenados en la matriz dshape, mientras que los términos 𝜕𝜉1 𝜕𝑥 son las componentes de la matriz jacobiana inversa (𝑱−1), almacenadas en la matriz xjacinv. A continuación se muestra el código de la función kBmat2DLin, en el cual se implementa todo el cálculo desarrollado en este subapartado: function [Bmat2] = kBmat2DLin(ndim,nnode,ngauss,shapef,dshape,xjac,xjacinv,xi,eta,xr ef) % xi y eta son coordenadas naturales, hay que transformarlas al espacio % real: x1_real=N1(xi,eta)*x1_nodo1+N2(xi,eta)*x1_nodo2... x1_real=0; x2_real=0; for i=1:8 x1_real=x1_real+shapef(i)*xref(1,i); x2_real=x2_real+shapef(i)*xref(2,i); end %inicialización Bmat=zeros(6,ndim*nnode); % El operador lineal B tiene el siguiente aspecto % 0 0 0 0 | Nk,x 0 0| % 0 0 0 0 | 0 Nk,y 0| % x1 x2 1 0 | 0 0 0| % 0 0 0 -x1 | 0 0 Nk,x| % 0 0 0 x2 | 0 0 Nk,y| % 0 0 0 0 | Nk,y Nk,x 0| % dN dN dxi dN deta %---- = ----*---- + ----*---- % dx dxi dx deta dx kint = 1; for i=1:nnode Bmat(1,kint) = dshape(1,i)*xjacinv(1,1)+ ... dshape(2,i)*xjacinv(1,2); Bmat(2,kint+1)= dshape(1,i)*xjacinv(2,1)+... dshape(2,i)*xjacinv(2,2); % 0 0 Nk,2| Bmat(4,kint+2) = Bmat(2,kint+1); % 0 0 Nk,1| Bmat(5,kint+2) = Bmat(1,kint); % Nk,2 Nk,1 0| 15 15 Implementación de una nueva formulación de Deformación Plana Generalizada en un código de elementos finitos en MATLAB Bmat(6,kint) = Bmat(2,kint+1); Bmat(6,kint+1) = Bmat(1,kint); kint = kint + 3; end A=zeros(6,4); A(3,1)=x1_real; A(3,2)=x2_real; A(3,3)=1; A(4,4)=-x1_real; A(5,4)=x2_real; Bmat2=[A Bmat]; 5.1.6 Cálculo de la matriz de rigidez La matriz de rigidez es el resultado del siguiente cálculo: 𝑘𝑒=∫𝑩𝑻𝑪𝑩𝑑Ω Ω𝑒 ( 18 ) Ω𝑒 es el dominio del elemento. La matriz 𝑩 contiene las derivadas de las funciones de forma y su cálculo se ha detallado en el apartado 5.1.5. La matriz 𝑪 es la matriz constitutiva del material, y contiene las constantes elásticas que caracterizan el comportamiento del material. En un primer caso, se va a tomar un material isótropo por simplicidad, por lo que en la matriz 𝑪 solo entrarán en juego dos constantes elásticas (𝐸 y 𝜈). La matriz constitutiva será de la siguiente manera: 𝑪= [ 1𝐸−𝜈 𝐸−𝜈 𝐸0 0 0 −𝜈 𝐸1𝐸−𝜈 𝐸0 0 0 −𝜈 𝐸−𝜈 𝐸1𝐸0 0 0 0 0 0 2(1+𝜈) 𝐸0 0 0 0 0 0 2(1+𝜈) 𝐸0 0 0 0 0 0 2(1+𝜈) 𝐸 ] −1 ( 19 ) 5.1.6.1 Cuadratura de Gauss Para realizar el cálculo de la integral en la ecuación ( 18 ) de una manera computacionalmente eficiente y con un excelente grado de aproximación, se utiliza la técnica de la cuadratura de Gauss, que a continuación se detalla: ∫ ∫𝐹(𝜉1,𝜉2)𝑑𝜉1𝑑𝜉2 1 −1 1 −1 = ∑∑𝑤𝑖𝑤𝑗𝐹(𝜉1𝑖,𝜉2𝑗) 𝑛𝑝𝑔2 𝑗 𝑛𝑝𝑔1 𝑖 ( 20 ) Como se aprecia en la ecuación ( 20 ), la técnica de la cuadratura de Gauss para realizar integrales consiste Pseudocódigo 16 16 en calcular dicha integral mediante la suma de una serie de productos de uno pesos asignados a los puntos de Gauss (𝑤𝑖𝑤𝑗) por el valor de la función en dicho punto de Gauss (𝐹(𝜉1𝑖,𝜉2𝑗)). La suma de este producto en cada punto de integración da como resultado la integral de 𝐹 en el dominio del elemento. La integral que da como resultado la matriz de rigidez se resuelve en el espacio isoparamétrico (𝜉1,𝜉2), ya que en él los límites de integración son muy sencillos (en el espacio isoparamétrico todos los elementos son cuadrados de lado 2, y los límites de integración van de -1 a 1, mientras que en el espacio real los límites de integración serían los contornos parabólicos del elemento real). Para expresar esta integral en el espacio isoparamétrico, hay que transformar el diferencial de área de elemento 𝑑Ω, siguiendo la siguiente relación: 𝑑Ω=𝑑𝑥𝑑𝑦=𝐽𝑑𝜉1𝑑𝜉2 ( 21 ) Siendo 𝐽 el determinante de la matriz jacobiana. El cálculo de la matriz de rigidez se realiza en las siguientes líneas de código: %-------- CmatVoigt ------------- E=210e9; %acero nu=0.3; %acero CmatVoigt=inv([1/E -nu/E -nu/E 0 0 0 -nu/E 1/E -nu/E 0 0 0 -nu/E -nu/E 1/E 0 0 0 0 0 0 2*(1+nu)/E 0 0 0 0 0 0 2*(1+nu)/E 0 0 0 0 0 0 2*(1+nu)/E]); %--------------------------------- factor=wxi*weta*h*detxjac; %wi*wj*|J| Kloop=BmatLin'*CmatVoigt*BmatLin*factor; % aportación a la matriz de rigidez del punto de integración "kip" Kelement=Kelement+Kloop; % acumulación de la matriz de rigidez del elemento Kelement es la matriz del elemento calculada por la función StiffnessMatrixISO y es el argumento de salida. Una vez que se ha implementado la rutina que calcula la matriz de rigidez de un elemento, se procede a la programación de una rutina de ensamblaje, para resolver problemas con más de un elemento. 5.2 Rutina de ensamblaje La rutina de ensamblaje tiene como objetivo calcular la matriz de rigidez global de una malla de elementos. Esta se realiza calculando la matriz de rigidez de cada elemento y añadiendo la aportación de esta matriz elemental a la matriz global en las posiciones que correspondan a los grados de libertad del elemento en cuestión. Se comienza declarando dos variables globales, ya que serán utilizadas en esta rutina y en varias subrutinas. Estas dos variables almacenarán una matriz de nodos (número de nodo y sus coordenadas) y otra matriz de elementos (número de elemento y nodos que lo forman, la conectividad), como se muestra a continuación: clear all close all global elements global nodes A continuación se procede a la generación de las matrices elements (𝑛𝑒×9) y nodes (𝑛𝑛×4). Esto se realiza importando los datos desde un fichero .inp de Abaqus®. 17 17 Implementación de una nueva formulación de Deformación Plana Generalizada en un código de elementos finitos en MATLAB Es decir, previamente se ha debido crear un malla de nodos y elementos utilizando el software de elementos finitos Abaqus®. Una vez creada, se genera el archivo .inp correspondiente que contiene, entre otras cosas, las matrices deseadas, elements y nodes. La lectura del archivo .inp y la extracción de dichas matrices es realizada por la función [nodes,elements] =readInp(''). El código de esta subrutina puede ser encontrado en el Anexo 1. Una vez obtenidas las matrices anteriores, se continúa con la generación de una serie de variables que serán necesarias a lo largo del código: [nodes,elements]=readInp('Job-2.inp.txt'); % Función que lee un inp de Abaqus y devuelve los nodos y la conectividad de los elementos numnode=length(nodes(:,1)); %número de nodos numelem=length(elements(:,1)); %número de elementos ndof=3; %número de grados de libertad de cada nodo %total number of unknowns = number of nodes of the model * number of dof per node totalUnknown = numnode*ndof+4; % +4 para A,B,C,D desp=zeros(totalUnknown,1); F_ext = zeros(totalUnknown,1); %vector de fuerzas externas F_int = zeros(totalUnknown,1); %vector de fuerzas internas Kglobal = zeros(totalUnknown,totalUnknown); % Kglobal = K; %matriz de rigidez global Fglobal=zeros(totalUnknown,1); Una vez generadas las variables anteriores, se entra en un bucle que recorre todos los elementos uno por uno. En cada pasada se calculan las posiciones que corresponden a ese elemento dentro de la matriz de rigidez global, se calcula la matriz de rigidez de ese elemento y se asignan los valores correspondientes a cada posición dentro de la matriz de rigidez global. Se muestra el código: % bucle sobre los elementos for iel = 1 : numelem %contador de elementos sctrB = assemblyPHN(iel); %función que devuelve los índices de la matriz global en función de la conectividad de los elementos Kedd = zeros(length(sctrB),length(sctrB)); Fedd = zeros(length(sctrB),1); [Kedd]=StiffnessMatrixISO(nodes,elements,iel); % Devuelve Kedd (matriz de rigidez del elemento "iel") incluyendo A,B,C,D %-------------- Assembly ----------------------- Kglobal(sctrB,sctrB) = Kglobal(sctrB,sctrB) + Kedd; F_ext(sctrB,1) = Fglobal(sctrB,1) + Fedd; end iel es el contador de elementos, que va desde el elemento 1 hasta el numelem. La variable sctrB es un vector que contiene las posiciones que corresponden al elemento iel dentro de la matriz de rigidez global. Este vector se calcula en la función assemblyPHN(iel,elements), que se detalla a continuación: function [sctrB] = assemblyPHN(e) global elements % esta variable global “elements” contiene el número del element y su conectividad sctr = elements(e,2:9); nn = length(sctr); Pseudocódigo 18 18 % sctr contiene los nodos que conforman el elemento e. A continuación se va % a construir un vector a partir de estos nodos, cuyas componentes % representan las posiciones de los 3 grados de libertad de dichos nodos en % la matriz global. Estas posiciones serán 3*nodo+4(A,B,C,D) y - 2(si es el % grado de libertad x), -1(si es el grado de libertad y) o 0 (si es elgrado de libertad z) for k = 1 : nn sctrBfem(3*k-2) = 3*sctr(k)-2+4; %3 is the number of degrees of freedom per node sctrBfem(3*k-1) = 3*sctr(k)-1+4; sctrBfem(3*k) = 3*sctr(k)+4; end sctrB=[1 2 3 4 sctrBfem]; % a los grados de libertad de los nodos, se añaden los 4 grados de libertad del nodo "maestro" % de tal forma que se acumulan en las cuatro primeras posiciones las aportaciónes de cada elemento a los grados % de libertad A, B, C y D. Una vez conocidas las posiciones dentro de la matriz de rigidez global, se calcula la matriz de rigidez del elemento iel y se introducen en la matriz de rigidez global en las posiciones que corresponden: [Kedd]=StiffnessMatrixISO(nodes,elements,iel); Kglobal(sctrB,sctrB) = Kglobal(sctrB,sctrB) + Kedd; 19 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 -0.2 0 0.2 0.4 0.6 0.8 1 Figura 4.Situación deformada (−) e indeformada (- - -) 6 VALIDACIÓN DE RESULTADOS En este apartado se va a proceder a la validación de la matriz de rigidez de un elemento para comprobar que está correctamente calculada. Esto se va a realizar mediante varios problemas “test” con comportamiento isótropo, resolviéndolos a mano y con el elemento, y comprobando que el resultado es el mismo. El elemento tiene unas dimensiones de 1 x 1 y unas constantes elásticas de 𝐸=210𝑒9 𝑃𝑎 y 𝜈=0.3. 6.1 Ensayo de tracción Se aplica una tensión de tracción de 5 𝐺𝑃𝑎 a un elemento de dimensiones 1𝑚×1𝑚 (Figura 5) en su lado superior, teniendo apoyo de bolas en el lado inferior y estando impedidos los movimientos como sólido rígido. La solución analítica es la siguiente: 𝜎𝑖𝑗=(0 0 0 5·109)⟹𝜀𝑖𝑗=(−9.2857·10−3 0 0 0.021667)⟹{𝑢𝑥=−9.2857·10−3 𝑥 𝑢𝑦=0.021667 𝑦 ( 22 ) Los resultados en desplazamientos (en metros) son: Solución numérica Solución analítica Error relativo Nodo 𝑢𝑥 𝑢𝑦 𝑢𝑥 𝑢𝑦 𝑢𝑥 𝑢𝑦 1 0 0 0 0 0 0 2 -0,00928571 0 -0,00928571 0 5,60E-16 0 3 -0,00928571 0,02166667 -0,00928571 0,02166667 3,74E-16 4,80E-16 4 0 0,02166667 0 0,02166667 0 6,41E-16 5 -0,00464286 0 -0,00464286 0 9,34E-17 0 6 -0,00928571 0,01083333 -0,00928571 0,01083333 0 4,00E-16 Figura 5. Elemento sometido a ensayo de tracción Validación de resultados 20 20 7 -0,00464286 0,02166667 -0,00464286 0,02166667 5,60E-16 3,20E-16 8 -6,26E-18 0,01083333 0 0,01083333 6,74E-16 1,60E-16 Tabla 2. Comparación de resultados del ensayo 6.1 6.2 Bitracción A un elemento de las mismas dimensiones que el del problema 6.1 se le aplica el siguiente campo de desplazamientos en los nodos: {𝑢𝑥=𝑥 𝑢𝑦=𝑦 𝑢𝑧=0⟹𝜀𝑖𝑗=(1 0 0 0 1 0 0 0 0)⟹𝜎𝑖𝑗=(4.0385·1011 0 0 0 4.0385·1011 0 0 0 2.4231·1011) ( 23 ) A coontinuación se muestra la comparación de resultados analíticos y numéricos: Nodo 1 2 3 4 5 6 7 8 𝜎𝑥𝑥 Analítica 4.0385e+11 4.0385e+11 4.0385e+11 4.0385e+11 4.0385e+11 4.0385e+11 4.0385e+11 4.0385e+11 Numérica 4.0385e+11 4.0385e+11 4.0385e+11 4.0385e+11 4.0385e+11 4.0385e+11 4.0385e+11 4.0385e+11 Diferencia 1.5113e-16 1.5113e-16 1.5113e-16 1.5113e-16 1.5113e-16 1.5113e-16 1.5113e-16 1.5113e-16 𝜎𝑦𝑦 Analítica 4.0385e+11 4.0385e+11 4.0385e+11 4.0385e+11 4.0385e+11 4.0385e+11 4.0385e+11 4.0385e+11 Numérica 4.0385e+11 4.0385e+11 4.0385e+11 4.0385e+11 4.0385e+11 4.0385e+11 4.0385e+11 4.0385e+11 Diferencia 1.5113e-16 1.5113e-16 1.5113e-16 1.5113e-16 1.5113e-16 1.5113e-16 1.5113e-16 1.5113e-16 𝜎𝑧𝑧 Analítica 2.4231e+11 2.4231e+11 2.4231e+11 2.4231e+11 2.4231e+11 2.4231e+11 2.4231e+11 2.4231e+11 Numérica 2.4231e+11 2.4231e+11 2.4231e+11 2.4231e+11 2.4231e+11 2.4231e+11 2.4231e+11 2.4231e+11 Diferencia 1.5113e-16 1.5113e-16 1.5113e-16 1.5113e-16 1.5113e-16 1.5113e-16 1.5113e-16 1.5113e-16 𝜎𝑥𝑦 Analítica 0 0 0 0 0 0 0 0 Numérica 0 0 0 0 0 0 0 0 Diferencia 0 0 0 0 0 0 0 0 𝜎𝑥𝑧 Analítica 0 0 0 0 0 0 0 0 Numérica 0 0 0 0 0 0 0 0 Diferencia 0 0 0 0 0 0 0 0 𝜎𝑦𝑧 Analítica 0 0 0 0 0 0 0 0 Numérica 2.6902e-05 2.6902e-05 2.6902e-05 2.6902e-05 2.6902e-05 2.6902e-05 2.6902e-05 2.6902e-05 Diferencia 6.6613e-17 6.6613e-17 6.6613e-17 6.6613e-17 6.6613e-17 6.6613e-17 6.6613e-17 6.6613e-17 Tabla 3. Comparación entre resultados analíticos y numéricos del problema 6.2 21 21 Implementación de una nueva formulación de Deformación Plana Generalizada en un código de elementos finitos en MATLAB 6.3 Ensayo de cortadura Se impone el siguiente campo de desplazamientos a un elemento similar al de los problemas anteriores: {𝑢𝑥=𝑦 𝑢𝑦=𝑥 𝑢𝑧=0⟹𝜀𝑖𝑗=(0 2 0 2 0 0 0 0 0)⟹𝜎𝑖𝑗=( 0 1.6154·1011 0 1.6154·1011 0 0 0 0 0) ( 24 ) Nodo 1 2 3 4 5 6 7 8 𝜎𝑥𝑥 Analítica 0 0 0 0 0 0 0 0 Numérica 7.3979e-05 7.3979e-05 7.3979e-05 7.3979e-05 7.3979e-05 7.3979e-05 7.3979e-05 7.3979e-05 Diferencia 4.5797e-16 4.5797e-16 4.5797e-16 4.5797e-16 4.5797e-16 4.5797e-16 4.5797e-16 4.5797e-16 𝜎𝑦𝑦 Analítica 0 0 0 0 0 0 0 0 Numérica -6.7254e-06 -6.7254e-06 -6.7254e-06 -6.7254e-06 -6.7254e-06 -6.7254e-06 -6.7254e-06 -6.7254e-06 Diferencia 4.1633e-17 4.1633e-17 4.1633e-17 4.1633e-17 4.1633e-17 4.1633e-17 4.1633e-17 4.1633e-17 𝜎𝑧𝑧 Analítica 0 0 0 0 0 0 0 0 Numérica 2.0176e-05 2.0176e-05 2.0176e-05 2.0176e-05 2.0176e-05 2.0176e-05 2.0176e-05 2.0176e-05 Diferencia 1.249e-16 1.249e-16 1.249e-16 1.249e-16 1.249e-16 1.249e-16 1.249e-16 1.249e-16 𝜎𝑥𝑧 Analítica 0 0 0 0 0 0 0 0 Numérica 0 0 0 0 0 0 0 0 Diferencia 0 0 0 0 0 0 0 0 𝜎𝑦𝑧 Analítica 0 0 0 0 0 0 0 0 Numérica 0 0 0 0 0 0 0 0 Diferencia 0 0 0 0 0 0 0 0 𝜎𝑥𝑦 Analítica 1.6154e+11 1.6154e+11 1.6154e+11 1.6154e+11 1.6154e+11 1.6154e+11 1.6154e+11 1.6154e+11 Numérica 1.6154e+11 1.6154e+11 1.6154e+11 1.6154e+11 1.6154e+11 1.6154e+11 1.6154e+11 1.6154e+11 Diferencia 1.8892e-16 1.8892e-16 1.8892e-16 1.8892e-16 1.8892e-16 1.8892e-16 1.8892e-16 1.8892e-16 Tabla 4. Comparación entre resultados numéricos y analíticos del problema 6.3 Validación de resultados 28 28 2 2 4,0385E+11 4,0385E+11 0,00054932 2 3 4,0385E+11 4,0385E+11 0,00146484 2 4 4,0385E+11 4,0385E+11 0,00231934 2 5 4,0385E+11 4,0385E+11 -0,00085449 2 6 4,0385E+11 4,0385E+11 -0,00286865 2 7 4,0385E+11 4,0385E+11 0,00128174 2 8 4,0385E+11 4,0385E+11 0,00079346 2 9 4,0385E+11 4,0385E+11 -0,00109863 3 1 4,0385E+11 4,0385E+11 -0,00061035 3 2 4,0385E+11 4,0385E+11 -0,00158691 3 3 4,0385E+11 4,0385E+11 0,00024414 3 4 4,0385E+11 4,0385E+11 -0,00042725 3 5 4,0385E+11 4,0385E+11 -0,00073242 3 6 4,0385E+11 4,0385E+11 -0,00018311 3 7 4,0385E+11 4,0385E+11 -0,00018311 3 8 4,0385E+11 4,0385E+11 0,00018311 3 9 4,0385E+11 4,0385E+11 0 4 1 4,0385E+11 4,0385E+11 0 4 2 4,0385E+11 4,0385E+11 0 4 3 4,0385E+11 4,0385E+11 0,00012207 4 4 4,0385E+11 4,0385E+11 6,10E-05 4 5 4,0385E+11 4,0385E+11 0 4 6 4,0385E+11 4,0385E+11 6,10E-05 4 7 4,0385E+11 4,0385E+11 -0,00036621 4 8 4,0385E+11 4,0385E+11 0 4 9 4,0385E+11 4,0385E+11 0,00012207 5 1 4,0385E+11 4,0385E+11 -6,10E-05 5 2 4,0385E+11 4,0385E+11 0,00018311 5 3 4,0385E+11 4,0385E+11 -0,00012207 5 4 4,0385E+11 4,0385E+11 0,00012207 5 5 4,0385E+11 4,0385E+11 -6,10E-05 5 6 4,0385E+11 4,0385E+11 0 29 29 Implementación de una nueva formulación de Deformación Plana Generalizada en un código de elementos finitos en MATLAB 5 7 4,0385E+11 4,0385E+11 -6,10E-05 5 8 4,0385E+11 4,0385E+11 0,00024414 5 9 4,0385E+11 4,0385E+11 0,00048828 Tabla 11. Comparación de resultados analíticos y numéricos del problema 6.10 Se puede apreciar que las diferencias en el cálculo numérico de las tensiones con respecto a los cálculos analíticos son muy pequeñas en comparación al valor de estas. Se ha elegido mostrar solo los valores de 𝜎𝑥 por ser ésta la mayor tensión y para simplificar la tabla. 6.11 Test de la parcela 3 (ensayo de cortadura) En este problema, la misma malla utilizada en el problema 6.10 es sometida al campo de desplazamientos del problema 6.3. De nuevo, se comparan los resultados en tensiones analítico y numérico en los puntos de integración y se muestra la diferencia entre ellos. De esta forma se pone en evidencia que la malla de elementos realiza el cálculo numérico de una forma fiable y precisa. Elemento Punto de Gauss 𝜎𝑥𝑦 analítica 𝜎𝑥𝑦 numérica Diferencia 1 1 1,6154E+11 1,6154E+11 0 1 2 1,6154E+11 1,6154E+11 0 1 3 1,6154E+11 1,6154E+11 0,00012207 1 4 1,6154E+11 1,6154E+11 -3,05E-05 0 0.5 1 1.5 2 0 0.5 1 1.5 2 2.5 3 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 Figura 11. Puntos de Gauss 0 0.5 1 1.5 2 2.5 0 0.5 1 1.5 2 2.5 3 Figura 12. Situación deformada e indeformada del problema 6.11 Validación de resultados 30 30 1 5 1,6154E+11 1,6154E+11 0 1 6 1,6154E+11 1,6154E+11 3,05E-05 1 7 1,6154E+11 1,6154E+11 3,05E-05 1 8 1,6154E+11 1,6154E+11 6,10E-05 1 9 1,6154E+11 1,6154E+11 -0,00018311 2 1 1,6154E+11 1,6154E+11 -0,00015259 2 2 1,6154E+11 1,6154E+11 0,00033569 2 3 1,6154E+11 1,6154E+11 0,00015259 2 4 1,6154E+11 1,6154E+11 0,00112915 2 5 1,6154E+11 1,6154E+11 -0,00012207 2 6 1,6154E+11 1,6154E+11 -0,00045776 2 7 1,6154E+11 1,6154E+11 0,00036621 2 8 1,6154E+11 1,6154E+11 -0,00021362 2 9 1,6154E+11 1,6154E+11 -0,00079346 3 1 1,6154E+11 1,6154E+11 -0,00036621 3 2 1,6154E+11 1,6154E+11 -0,00021362 3 3 1,6154E+11 1,6154E+11 -0,00079346 3 4 1,6154E+11 1,6154E+11 -0,00033569 3 5 1,6154E+11 1,6154E+11 -0,00045776 3 6 1,6154E+11 1,6154E+11 -0,00048828 3 7 1,6154E+11 1,6154E+11 -6,10E-05 3 8 1,6154E+11 1,6154E+11 -6,10E-05 3 9 1,6154E+11 1,6154E+11 0,00021362 4 1 1,6154E+11 1,6154E+11 -3,05E-05 4 2 1,6154E+11 1,6154E+11 3,05E-05 4 3 1,6154E+11 1,6154E+11 0 4 4 1,6154E+11 1,6154E+11 9,16E-05 4 5 1,6154E+11 1,6154E+11 6,10E-05 4 6 1,6154E+11 1,6154E+11 3,05E-05 4 7 1,6154E+11 1,6154E+11 -6,10E-05 4 8 1,6154E+11 1,6154E+11 -6,10E-05 4 9 1,6154E+11 1,6154E+11 0,00018311 31 31 Implementación de una nueva formulación de Deformación Plana Generalizada en un código de elementos finitos en MATLAB 5 1 1,6154E+11 1,6154E+11 0 5 2 1,6154E+11 1,6154E+11 0,00012207 5 3 1,6154E+11 1,6154E+11 -6,10E-05 5 4 1,6154E+11 1,6154E+11 -3,05E-05 5 5 1,6154E+11 1,6154E+11 -3,05E-05 5 6 1,6154E+11 1,6154E+11 0,00018311 5 7 1,6154E+11 1,6154E+11 -3,05E-05 5 8 1,6154E+11 1,6154E+11 -3,05E-05 5 9 1,6154E+11 1,6154E+11 0,00021362 Tabla 12. Comparación de resultados del test de la parcela del problema 6.11 6.12 Test de la parcela 4 (campo de desplazamientos parabólico) En este problema, la misma malla utilizada en el problema 6.10 es sometida al campo de desplazamientos del problema 6.4. De nuevo, se comparan los resultados en tensiones analítico y numérico en los puntos de Gauss y se muestra la diferencia. 0 0.5 1 1.5 2 2.5 0 0.5 1 1.5 2 2.5 3 Figura 13. Situación deformada e indeformada del problema 6.12 Validación de resultados 32 32 Elemento Punto de Gauss 𝜎𝑥 analítica 𝜎𝑥 numérica Diferencia Error relativo 1 1 1,4575E+11 1,3467E+11 1,1073E+10 0,01008333 1 2 5,5901E+11 5,6032E+11 -1304497524 -0,00118786 1 3 9,7228E+11 9,8421E+11 -1,1935E+10 -0,0108679 1 4 2,0866E+11 1,7043E+11 3,8228E+10 0,0348102 1 5 5,3712E+11 5,3712E+11 6,10E-05 5,56E-17 1 6 8,6557E+11 9,0013E+11 -3,4552E+10 -0,0314623 1 7 2,7157E+11 2,4022E+11 3,1343E+10 0,02854072 1 8 5,1522E+11 4,983E+11 1,6922E+10 0,01540871 1 9 7,5887E+11 7,7088E+11 -1,2015E+10 -0,01094043 2 1 8,3542E+11 8,3497E+11 452599768 0,00041213 2 2 9,6434E+11 9,6417E+11 174552759 0,00015895 2 3 1,0933E+12 1,0936E+12 -315681947 -0,00028746 2 4 8,5485E+11 8,5384E+11 1016074801 0,00092523 2 5 9,7529E+11 9,7529E+11 -0,00012207 -1,11E-16 2 6 1,0957E+12 1,0967E+12 -996778397 -0,00090765 2 7 8,7428E+11 8,7357E+11 714776970 0,00065087 2 8 9,8624E+11 9,8582E+11 414017976 0,000377 2 9 1,0982E+12 1,0985E+12 -305746621 -0,00027841 3 1 7,9773E+11 8,1713E+11 -1,9401E+10 -0,01766658 3 2 8,8747E+11 9,1458E+11 -2,7111E+10 -0,02468688 3 3 9,7721E+11 9,7721E+11 -0,00024414 -2,22E-16 3 4 5,1522E+11 5,3435E+11 -1,9132E+10 -0,01742103 3 5 5,3712E+11 5,3712E+11 -6,10E-05 -5,56E-17 3 6 5,5901E+11 5,3275E+11 2,6267E+10 0,02391808 3 7 2,3271E+11 2,3271E+11 0 0 3 8 1,8676E+11 1,6846E+11 1,8304E+10 0,01666723 3 9 1,4081E+11 1,2279E+11 1,8021E+10 0,01640966 4 1 2,477E+10 2,43E+10 469707104 0,00042771 4 2 1,0989E+11 1,1016E+11 -266080183 -0,00024229 4 3 1,9501E+11 1,9567E+11 -662278539 -0,00060306 33 33 Implementación de una nueva formulación de Deformación Plana Generalizada en un código de elementos finitos en MATLAB 4 4 2,2302E+10 2,0738E+10 1564089428 0,00142424 4 5 9,8942E+10 9,8942E+10 0 0 4 6 1,7558E+11 1,7715E+11 -1564089428 -0,00142424 4 7 1,9834E+10 1,9241E+10 592835188 0,00053983 4 8 8,7994E+10 8,8345E+10 -351505093 -0,00032008 4 9 1,5615E+11 1,5709E+11 -936538437 -0,0008528 5 1 2,8494E+11 2,8425E+11 691127671 0,00062933 5 2 5,0885E+11 5,0921E+11 -366584882 -0,00033381 5 3 7,3275E+11 7,3376E+11 -1007298339 -0,00091723 5 4 2,6798E+11 2,6602E+11 1956845201 0,00178188 5 5 5,0885E+11 5,0885E+11 6,10E-05 5,56E-17 5 6 7,4972E+11 7,5174E+11 -2021664461 -0,0018409 5 7 2,5101E+11 2,5006E+11 954174209 0,00086886 5 8 5,0885E+11 5,0831E+11 539971424 0,00049169 5 9 7,6668E+11 7,6723E+11 -549952385 -0,00050078 La exactitud de los cálculos en este problema disminuye debido a la distrosión de la malla. Al ser esta malla más irregular, el campo de desplazamientos parabólico no se adapta tan bien a las funciones de forma parabólicas del espacio isoparamétrico. Aun así, el máximo error relativo cometido, tomando como referencia la máxima tensión 𝜎𝑥 analítica (1098191204282,88 Pa) es del 3.48 %, lo cual es bastante aceptable. 35 7 CONCLUSIONES Y DESARROLLOS FUTUROS Una vez validados los resultados en el capítulo 6, se puede decir que la formulación del elemento está correctamente implementada y funciona. Como ya se ha visto, todo el código ha sido escrito en MATLAB ® por la facilidad que supone la depuración e identificación de errores en esta herramienta. Sin embargo, MATLAB ® no es eficiente a la hora de definir las condiciones de contorno, cargas, geometría, malla de nodos y elementos del problema. Por tanto, el siguiente paso en el futuro desarrollo de este proyecto sería realizar una traducción del código a lenguaje de programación Fortran para su implementación en un código comercial de elementos finitos como ABAQUS ®. 37 REFERENCIAS [1] A. Blázquez, V. Mantič and F. París, Application of BEM to generalized plane problems for anisotropic elastic materials in presence of contact, Engineering Analysis with Boundary Elements, vol. 30, pp. 489- 502, 2006. [2] J. Reinoso, Elementos Finitos multi-dimensionales: formulaciones 2D y 3D para sólidos elásticos lineales, Sevilla: Universidad de Sevilla, 2016. [3] O. C. Zienkiewicz, El método de los elementos finitos, Vol. 1, Madrid: McGraw-Hill, 1994. Anexo 3 44 44 shapef(5) = 1/2*(1-X)*(1+X)*(1-Y); dshape(1,5) = X*(-1+Y); dshape(2,5) = 1/2*(1+X)*(-1+X); shapef(6) = 1/2*(1+X)*(1+Y)*(1-Y); dshape(1,6) = -1/2*(1+Y)*(-1+Y); dshape(2,6) = -Y*(1+X); shapef(7) = 1/2*(1-X)*(1+X)*(1+Y); dshape(1,7) = -X*(1+Y); dshape(2,7) = -1/2*(1+X)*(-1+X); shapef(8) = 1/2*(1-X)*(1+Y)*(1-Y); dshape(1,8) = 1/2*(1+Y)*(-1+Y); dshape(2,8) = Y*(-1+X); end La función kJacobian2D calcula el jacobiano de las deformaciones: %========================================================== function [xjac,xjacinv,detxjac]=kJacobian2D(ndim,nnode,shapef,dshape,xref) %initialization phase for i=1:ndim for j=1:ndim xjac(i,j) = 0.0; xjacinv(i,j) =0.0; end end % xref son las coordenadas de los nodos del elemento en el espacio real for i=1:nnode xjac(1,1) = xjac(1,1) + dshape(1,i)*xref(1,i); xjac(1,2) = xjac(1,2) + dshape(1,i)*xref(2,i); xjac(2,1) = xjac(2,1) + dshape(2,i)*xref(1,i); xjac(2,2) = xjac(2,2) + dshape(2,i)*xref(2,i); end %COMPUTE DETERMINANT OF THE JACOBIAN detxjac=det(xjac); if(detxjac > 0.0) %It can be formulated the inverse %COMPUTE THE INVERSE xjacinv = inv(xjac); else end 45 45 Implementación de una nueva formulación de Deformación Plana Generalizada en un código de elementos finitos en MATLAB end %========================================================== La función kBmat2DLin calcula el operador de compatibilidad 𝑩: function [Bmat2] = kBmat2DLin(ndim,nnode,ngauss,shapef,dshape,xjac,xjacinv,xi,eta,xref) % xi y eta son coordenadas naturales, hay que transformarlas al espacio % real: x1_real=N1(xi,eta)*x1_nodo1+N2(xi,eta)*x1_nodo2... x1_real=0; x2_real=0; for i=1:8 x1_real=x1_real+shapef(i)*xref(1,i); x2_real=x2_real+shapef(i)*xref(2,i); end %initialization Bmat=zeros(6,ndim*nnode); % Linear B-operator looks like that % 0 0 0 0 | Nk,x 0 0| % 0 0 0 0 | 0 Nk,y 0| % x1 x2 1 0 | 0 0 0| % 0 0 0 -x1 | 0 0 Nk,x| % 0 0 0 x2 | 0 0 Nk,y| % 0 0 0 0 | Nk,y Nk,x 0| % dN dN dxi dN deta %---- = ----*---- + ----*---- % dx dxi dx deta dx kint = 1; for i=1:nnode Bmat(1,kint) = dshape(1,i)*xjacinv(1,1)+ ... dshape(2,i)*xjacinv(1,2); Bmat(2,kint+1)= dshape(1,i)*xjacinv(2,1)+... dshape(2,i)*xjacinv(2,2); % 0 0 Nk,2| Bmat(4,kint+2) = Bmat(2,kint+1); % 0 0 Nk,1| Bmat(5,kint+2) = Bmat(1,kint); % Nk,2 Nk,1 0| Bmat(6,kint) = Bmat(2,kint+1); Bmat(6,kint+1) = Bmat(1,kint); kint = kint + 3; end Anexo 3 46 46 A=zeros(6,4); A(3,1)=x1_real; A(3,2)=x2_real; A(3,3)=1; A(4,4)=-x1_real; A(5,4)=x2_real; Bmat2=[A Bmat]; %==========================================================