resolução de area de um circulo por MEF no MATLAB

Prévia do material em texto

UNIVERSIDADE DO ESTADO DE MATO GROSSO
CAMPUS UNIVERSITÁRIO DE TANGARÁ DA SERRA
BACHARELADO EM ENGENHARIA CIVIL
USO DO MATLAB E MEF PARA CÁLCULO DA ÁREA DO CÍRCULO
Acadêmico: Jhonatta Andrade
Prof.º Renê Quispe Rodriguez
Tangará da Serra – MT
Setembro de 2017
1 INTRODUÇÃO
	O método dos elementos finitos (MEF ou FEM), consiste em uma ferramenta utilizada na resolução de problemas através de métodos numéricos, onde não é possível utilização de métodos analíticos. Este método foi desenvolvido com o objetivo de resolver problemas complexos em mecânica dos sólidos, porém o seu uso não se restringe à essa área, sendo utilizado em diversos problemas físicos complexos, tais como análise de fluxos de fluídos, entre outras aplicações.
	O MEF resolve problemas descritos por equações diferencias de maneira aproximada, utilizando para isso a aproximação do tipo polinomial nodal em subdomínios, através de discretização.
Polinomial: já que a aproximação se dá por meio de polinômios, que são capazes de representar o comportamento físico analisado em uma dada situação;
Nodal: a análise do comportamento é feita através de nós;
Subdomínio: é uma fração do domínio que será analisada por discretização;
Discretização: subdivisão do domínio de interesse em elementos finitos interconectados por nós;
Cada elemento finito conserva as propriedades do meio original;
	A resolução de um problema físico utilizando MEF segue, de maneira simplificada os passos a seguir: 
Passo 0 – Entender o problema físico a ser resolvido
Consiste em definir a Equação diferencial que descreve o comportamento físico, o fenômeno físico em si e analisar as condições de contorno do mesmo.
Passo 1 – Discretização e escolha do tipo de elemento
Neste passo a malha é gerada (elementos e nós), também são definidos o tamanho e tipos de elementos, que influem na complexidade e consequente precisão da aproximação, e o tipo de elemento a ser utilizado (1D, 2D e 3D).
Passo 2 – Escolha do modelo de Aproximação
É neste passo que a função de aproximação é definida, geralmente polinomiais ou trigonométricas, ou seja é escolhida a forma de distribuição da grandeza desconhecida (deslocamentos, temperaturas, pressão, etc.), de maneira a minimizar o erro na interpolação.
Passo 3 – definir o gradiente e as relações constitutivas 
Em mecânica por exemplo, pode-se necessitar das relações entre as deformações e deslocamentos e/ou tensões e deformações. 
Passo 4 – Obter as equações dos elementos 
As obtenções das equações seguem 2 enfoques principais: 
Métodos residuais 
Métodos baseados em energia
Sendo o segundo mais utilizado em aplicações em mecânica estrutural.
2 PROBLEMA
	Com o objetivo de utilizar de maneira prática o MEF, foi solicitado a criação de um Script em MATLAB para resolução de barras de uma forma geral. O problema proposto foi: 
Pede-se: 
Item 1: Deslocamentos nodais
Item 2: Tensões nos Elementos 
3 MÉTODO DE RESOLUÇÃO
Para Resolução do problema e criação do Script no Matlab foram seguidos os passos propostos em aula para a implementação computacional. 
Foram desenvolvidos 2 scripts, um para armazenar as informações referentes a barra, tais como (coordenadas dos nós, propriedades dos elementos, condições de contorno, incidência e Forças nodais) e outro com a rotina para determinação dos deslocamentos nodais e tensões em cada elemento.
Script 1 – Dados_barra
Neste Script temos diversas matrizes que determinam nosso problema.
[COORD] = matriz das coordenadas dos nós. [COORD] = [x,y]
[PGEOM] = matriz que armazena as propriedades geométricas dos elementos da barra. [PGEOM] = [A1, A2, A3]. 
[PMAT] = matriz que armazena as propriedades materiais dos elementos da barra.
[PMAT] = [E1, E2, E3] .
[INCI] = matriz de incidência do elemento. Esta matriz relaciona o elemento, seu tipo geométrico, tipo de material e os nós que conformam o elemento. 
[INCI] = [Elemento, Tipo Geométrico, Tipo Material, Nó 1, Nó 2]. Obs.: por se tratar de barras apenas 2 nós conformam cada elemento.
[CONT] = Matriz que armazena as condições de contorno de cada nó. Relaciona o nó com seu respectivo G.D.L. restrito (adotando-se 1 para restrição no eixo x, e 2 para restrição no eixo y).
[CONT] = [Nó, GDL Restrito]
[LOADS] = Matriz que armazena as forças nodais. Relaciona o nó, a direção da força (1 para x e 2 para y) e o valor da força. 
[LOADS] = [Nó, direção da Força, Valor da Força]
Importante: As unidades utilizadas devem ser equivalentes, para que haja consistência nos resultados apresentados pelo programa. Para tal utilizou-se as seguintes unidades de medidas.
Área: m²
Forças: N
E: N/m² (ou Pa)
Coordenadas dos nós: m
Para essas unidades os resultados apresentados para deslocamento e tensões são dados em metros e N/m² respectivamente.
Script 2 – programa_barras
Para realização dos cálculos referentes aos problemas de barras, é necessário algumas informações que serão usadas na criação das matrizes utilizadas na resolução do problema. Esses valores são: Número de Nós (NNOS) e Número de Elementos (NELEMENTOS). Esses valores podem ser quantificados de forma bastante simples no Matlab, bastando para tanto acessar a informação do número de linhas das matrizes [COORD] e [INCI], respectivamente, para tal utiliza-se o comando size(COORD,1) e size(INCI, 1)
O próximo passo seria determinar a matriz [ID] (matriz identificador array), essa matriz numera os GDL dos nós. Para o caso de barras são considerados apenas deslocamentos horizontais, desta forma a matriz [ID] possuirá apenas uma linha e o número de colunas será igual ao número de nós. 
 Para determinação da Matriz [ID] foram usados 3 passos: 
Criar uma matriz de uma linha e com o número de colunas igual ao número de nós, totalmente preenchida com 1. Comando: ID = ones(1,NNOS).
Trocar o valor da matriz [ID] por zero nos nós onde o GDL é Restrito. Utiliza-se um laço de repetição (for) que identifica o valor do nó restrito na matriz [CONT] e altera seu valor para zero, desta forma o nó restrito não irá entrar na análise dos deslocamentos.
Numerar sequencialmente os deslocamentos livres dos nós na matriz [ID]. Utiliza-se um laço de repetição (for) que numera de forma sequencial todos os nós cujo valor na matriz [ID] eram iguais a 1. 
Com a Matriz [ID] finalizada passou-se ao cálculo das matrizes de rigidez local e, de forma simultânea, através de três laços de repetição interligados. O primeiro deles determina a matriz de rigidez local de cada elemento [Ke] e o vetor de conectividade [LM], os dois últimos laços armazenam cada matriz de rigidez local em sua respectiva localização dentro da matriz de rigidez global [KG], somando os valores de cada posição caso haja sobreposição. 
4 RESULTADOS
Através da aplicação da fórmula supracitada no MATLAB, criou-se um script (Anexo A), que através de iterações analisa a convergência da área em função do número de subdivisões do círculo em triângulos isósceles (‘n’), desta forma pode-se observar, conforme o gráfico abaixo, a forma que essa convergência ocorre para uma circunferência de raio um. Em todos os gráficos apresentados abaixo, analisou-se, conforme solicitado, a convergência para que o erro seja inferior a 1%.
 Figura 2 - Convergência da área calculada em função do número de divisões (n)
 Fonte: Autor
Conforme pode-se observar no gráfico acima, a área converge de forma logarítmica, observando-se boa convergência com valores de n maiores que vinte. No caso acima foi utilizada uma circunferência de raio um, porém a convergência se comporta da mesma forma com circunferências de qualquer raio.
O gráfico abaixo relaciona o erro (%) em função de ‘n’ onde se pode observar que o erro se torna relativamente pequeno quando o valor de ‘n’ é maior que vinte. 
Figura 3- Erro percentual em função de n
 Fonte: Autor
A tabela a seguir faz uma relação entre os valores de N, área calculada, área real e o erro percentual, onde podemos concluirque para que o erro seja inferior a 1%, N deve ser igual ou superior a 26, o que fica claro quando se executa o script (Anexo A). 
Tabela 1 – Relação entre N, área calculada, área real e erro (%)
	N
	Área Calculada
	Área Real
	Erro (%)
	3
	1,299038106
	3,141592654
	58,65033284
	4
	2
	3,141592654
	36,33802276
	5
	2,377641291
	3,141592654
	24,31732714
	6
	2,598076211
	3,141592654
	17,30066569
	7
	2,736410189
	3,141592654
	12,89735843
	8
	2,828427125
	3,141592654
	9,968368384
	9
	2,892544244
	3,141592654
	7,927457104
	10
	2,938926261
	3,141592654
	6,451071621
	11
	2,973524496
	3,141592654
	5,349775611
	12
	3
	3,141592654
	4,507034145
 Continuação (...)
	13
	3,020700618
	3,141592654
	3,848113
	14
	3,037186174
	3,141592654
	3,323361469
	15
	3,050524823
	3,141592654
	2,898779077
	16
	3,061467459
	3,141592654
	2,55046416
	17
	3,070554163
	3,141592654
	2,261225398
	18
	3,07818129
	3,141592654
	2,018446395
	19
	3,084644957
	3,141592654
	1,812701468
	20
	3,090169944
	3,141592654
	1,636835692
	21
	3,094929331
	3,141592654
	1,48533968
	22
	3,099058125
	3,141592654
	1,353916087
	23
	3,102662868
	3,141592654
	1,239173552
	24
	3,105828541
	3,141592654
	1,138407053
	25
	3,10862359
	3,141592654
	1,049437902
	26
	3,111103636
	3,141592654
	0,97049558
Fonte: Autor
5 REFERÊNCIAS
Azevedo, A. F. M. Método dos elementos finitos. Faculdade de Engenharia da universidade do porto: Portugal, 2013.
Becker, A. J. et. al. Noções Básicas de programação em MATLAB. Pet matemática ufsm. Santa Maria, 2010.
6 ANEXO A – Script utilizado no MATLAB para resolução do problema proposto
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%% PROGRAMA CÁLCULO DE ÁREA DE UMA CIRCUNFERÊNCIA %%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
clear %limpa os valores das variáveis
clc % comando para limpar a tela
erro_adm = input('entre com o erro admissível (erro admissível > 0) = '); %recebe o valor para o erro admissível
clc; % comando para limpar a tela
while erro_adm <= 0 % laço de repetição que impede que o erro admissível seja menor ou igual a zero
erro_adm = input('entre com o erro admissível (erro admissível > 0) = '); %recebe o valor para o erro admissível
clc; % comando para limpar a tela
end
raio = input('entre com o raio (raio > 0)= ') ; %recebe o valor do raio
clc % comando para limpar a tela
while raio <= 0 % laço de repetição que impede que o raio seja menor ou igual a zero
 raio = input('entre com o raio (raio > 0)= ') ; %recebe o valor do raio
 clc % comando para limpar a tela
end
n = 3 ; % n mínimo = 3 (obs.: não foi dada a opção para o usuário setar o valor de n)
area_calc = raio^2*sin(pi/n)*cos(pi/n)*n; %calcula a área em função de n e r
area_real = pi*raio^2; %calcula a área real em função de r
erro = (area_real - area_calc)/area_real; % calcula o erro entre a área real e a área calculada
while erro >= erro_adm; %laço de repetição que incrementa n com passo 1 até que o erro seja menor ou igual ao erro admissível
 er(n-2) = erro % vetor que armazena o erro em função de n
 ac(n-2) = area_calc; % vetor que armazena a área calculada em função de n
 n = n + 1; % incrementa 1 em n
 area_calc = raio^2*sin(pi/n)*cos(pi/n)*n; % calcula a área em função de n e r
 erro = (area_real - area_calc)/area_real; % calcula o novo erro em função do novo n 
end
for c = 1:n-2 % laço para criação do vetor com a área real para posterior plotagem
 ar(c) = pi*raio^2 % área real da circunferência
end
ac(n-2) = area_calc; % adiciona o ultimo termo do vetor da área calculada
er(n-2) = erro % adiciona o ultimo termo do vetor do erro
t = 3:1:n; % vetor t, será usado para armazenar os valores de n
clc % comando para limpar a tela
'n : '
disp(n); % mostra o valor de n
'área real: '
disp(area_real); % mostra a área real
'área Calculada: '
disp(area_calc); % mostra a área calculada
'erro(%): '
disp(erro*100); % mostra o erro
 
%manipulação Gráfica e plotagem
 
subplot(2,2,1), plot(t, ar,'-*r');
xlabel('N');
ylabel('área Real');
grid minor
grid on
axis([0 n+5 0 area_real+5]);
 
subplot(2,2,2), plot(t, ac,'-*b');
xlabel('N');
ylabel('Área Calculada');
grid minor
grid on
axis([0 n+5 0 area_calc+5]);
 
subplot(2,2,3), plot(t, er*100, '-*k');
xlabel('N');
ylabel('Erro(%)');
grid minor
grid on
axis([0 n+5 0 100]);
 
subplot(2,2,4), plot(t,ac,'-*b',t,er,'-*k',t,ar,'-*r') 
xlabel('N');
legend('Área Calculada','Erro', 'Área Real', 0);
grid minor
grid on
axis([0 n+5 0 area_calc+5]);