Trabalho_06

Prévia do material em texto

Acadêmico: Geovane Augusto Haveroth - RA: 164573
Disciplina: Elementos Finitos I - IM381
Professor: Marco Lúcio Bittencourt
Trabalho 06 - Programa de barra com carga senoidal usando integração numérica.
Considere uma barra de comprimento L = 1m apoiada nas suas extremidades e submetida a uma carga
distribúıda com intensidade qx(x) definida por
qx(x) =
π
L
sin
(πx
L
)
(1)
conforme ilustrado na figura abaixo. Para este problema adotaremos EA = 1× 104.
Questão A: Solucionar este problema analiticamente (análogo ao trabalho 03).
Para um material elástico linear, homogêneo e isotrópico, tal problema é descrito pelo seguinte problema de
valor de contorno (BVP):
EA
d2ux(x)
dx2
= −qx(x), 0 < x < 1 (2)
com condições de contorno de Dirichlet
ux(0) = ux(1) = 0. (3)
Solucionamos o problema integrando duas vezes a eq. (2),
EAux(x) =
L
π
sin
(πx
L
)
+ c1x+ c2 (4)
e aplicando as condições de contorno (3) resulta em
ux(x) =
L
EAπ
sin
(πx
L
)
. (5)
Logo, a força normal é descrita como
Nx(x) = EA
dux(x)
dx
= cos
(πx
L
)
, (6)
com o diagrama de esforço normal dado pela figura abaixo.
Questão B: Solucionar este problema pelo método de elementos finitos (MEF) usando malhas de 2, 4, 8 e 16
elementos lineares e quadráticos.
Conforme apresentado no Trabalho 03, este problema se resume em solucionar o sistema de equações
Ka = f , (7)
1
sendo K a matriz de rigidez, a o vetor de incógnitas e f o vetor de forças globais. As matrizes e vetores globais
são obtidas aplicando o operador montagem A sobre a contribuição feita por cada elemento
K =
nelem
A
e=1
Ke e f =
nelem
A
e=1
fe. (8)
Na sequência apresentamos o procedimento adotado para avaliar Ke e fe numericamente.
Consideramos um conjunto de P+1 pontos de coordenadas ξp com p = (0, 1, . . . P ) no sistema de coordenada
local ξ. As funções de forma locais Φep(ξ) são definidas como o polinômio de Lagrange LPp (ξ) de ordem P
associado a coordenada ξp. Tais funções são expressas pelo produtório
Φep(ξ) := LPp (ξ) =
P∏
q=0,q 6=p
(ξ − ξq)
(ξp − ξq)
. (9)
Assim, escrevemos as funções de forma para um elemento linear de dois nós por
Φe0(ξ) =
1− ξ
2 e Φ
e
1(ξ) =
1 + ξ
2 , (10)
enquanto que para um elemento quadrático de três nós, estas funções são da forma
Φe0(ξ) = −
1
2ξ(1− ξ), Φ
e
1(ξ) = 1− ξ2 e Φe2(ξ) =
1
2ξ(1 + ξ). (11)
As derivadas de Φep(ξ) com relação à ξ em (10) e (11) são dadas respectivamente por
Φe0,ξ(ξ) = −
1
2 e Φ
e
1,ξ(ξ) =
1
2 , (12)
Φe0,ξ(ξ) = ξ −
1
2 , Φ
e
1,ξ(ξ) = −2ξ e Φe2,ξ(ξ) = ξ +
1
2 . (13)
Uma vez definido as funções de forma e suas respectivas derivadas, os coeficientes da matriz de rigidez do
elemento no sistema local de coordenadas ξ são obtidos por
Keij =
EA
he
∫ +1
−1
Φei,ξ(ξ) Φej,ξ(ξ)dξ, 0 ≤ i, j ≤ P, (14)
enquanto que os coeficientes do vetor de forças é descrito por
fei =
he
2
∫ +1
−1
qx
(
P∑
k=0
Φek(ξ)xek
)
Φei (ξ)dξ, 0 ≤ i ≤ P, (15)
onde xep é a coordenada do nó p.
Podemos integrar numericamente as expressões (14) e (15) utilizando Q pontos de integração ξm e penal-
izações wm. Estes pontos podem ou não coincidir com as coordenadas dos nós. As integrações numéricas dos
coeficientes da matriz de rigidez e forças do elemento são dadas respectivamente por
Keij ≈
EA
he
Q−1∑
m=0
Φei,ξ(ξm) Φej,ξ(ξm) wm, e fei ≈
he
2
Q−1∑
m=0
qx
(
P∑
k=0
Φek(ξm)xek
)
Φei (ξm) wm. (16)
Para solucionar o BVP dado inicialmente via MEF, devemos alterar no código a avaliação de Ke e fe, que
passa a ser numérico. Para tal, fazemos uso da função GetPointsWeight1D para extrair os pontos de integração
r e os respectivos pesos w :
[r,w] = GetPointsWeight1D(alpha,beta,Q,’GJ’,’NL’);
Adotamos ainda, por opção do acadêmico, a quadratura de Gauss-Jacobi com α = β = 0 e Q = 10. A opção de
tomar Q = 10 se deve ao fato de que Q pontos de integração são suficientes para integrar exatamente polinômios
de ordem 2Q−1 não é aplicável neste problema, visto que a função que estamos integrando qx é trigonométrica.
Vale salientar que para obter Ke exatamente, considerando elemento quadrático, bastaria tomar Q = 2. Na
sequência, as funções de forma Φep(ξm) e suas derivadas Φep,ξ(ξm) são avaliadas utilizando as seguintes linhas de
comando
2
if intorder==1 % Linear interpolation
dN = [-0.5+0*r 0.5+0*r];
N = [0.5*(1-r) 0.5*(1+r)];
elseif intorder==2 % Quadratic interpolation
dN = [-0.5+r -2*r 0.5+r];
N = [-0.5*r.*(1-r) 1-r.ˆ2 0.5*r.*(1+r)];
end
ou, generalizando a ordem de interpolação utilizando as seguintes funções
N = GetPolynomialValue(intorder,r,linspace(-1,1,Elnodes))’;
dN = Get1stDerivative(intorder,r,linspace(-1,1,Elnodes))’;
Uma vez calculado N e dN, avaliamos a matriz de rigidez do elemento conforme descrito anteriormente na
eq. (16) da forma
Ke = zeros(Elnodes,Elnodes);
for i = 1:Elnodes
for j = 1:Elnodes
for k = 1:length(w)
Ke(i,j) = Ke(i,j) + EA*dN(k,i)*dN(k,j)*w(k)*dJinv;
end
end
end
bem como o vetor de forças
Fe = zeros(Elnodes,1);
for i = 1:Elnodes
for k = 1:length(w)
Fe(i) = Fe(i) + (pi/Lb)*sin(pi*N(k,:)*xl/Lb)*N(k,i)*w(k)*dJac;
end
end
onde o termo dJac é o determinante do jacobiano definido por he/2, e dJinv sua inversa.
Questão C: Calcular os erros na norma de energia em cada caso e plotar um gráfico em escala log x log do erro
em função do tamanho h dos elementos das malhas. Avaliar o coeficiente angular das curvas de erros obtidos
para os elementos lineares e quadráticos.
A medida de erro considerada neste trabalho é avaliada com base na norma de energia ||.||A:
er =
||u||A − ||un||A
||u||A
, (17)
com a norma da aproximação un dada por
||un||A =
(
aTKa
) 1
2 , (18)
e a norma da solução exata por
||u||A =
(∫ L
0
EA(u′(x))2dx
) 1
2
. (19)
O procedimento de MEF descrito anteriormente é aplicado em malhas de 2, 4, 8, 16 elementos lineares e
quadráticos com nós igualmente espaçados. Esta escolha foi feita para poder realizar a comparação com o que
foi realizado no Trabalho 03. A figura 1 ilustra o erro relativo versus o tamanho da malha em escala Log-Log.
Não houve diferença significativa dos valores obtidos numericamente com aqueles avaliados de forma anaĺıtica
e impostas no código (a diferença de ≈ 10−8%). Desta forma, tanto a taxa de convergência quanto os demais
resultados apresentados anteriormente no Trabalho 03 são idênticos, e portanto omitidos na sequência.
3
Figure 1: Taxas de convergência.
Item D: O Algoritmo utilizado encontra-se descrito na sequência:
% Codigo Alterado Trabalho 06
clear
Lb = 1; % bar length
EA = 1e+4; % Young modulus Versus Cross sectional area
% FEM setup: all element with the same length
Nel = 16; % number of elements
intorder = 2; % interpolation order
Elnodes = intorder + 1; % number of nodes per element
Nvert = intorder*Nel+1; % number of nodes
he = Lb/Nel; % element lenght
Ndofs = intorder*Nel+1; % number of d.o.f. (1 d.o.f. per node)
% Allocation of the stifness matrix
Kg = zeros(Ndofs,Ndofs);
% Allocation of the load vector
Fg = zeros(Ndofs,1);
% Allocation of the solution vector
ug = zeros(Ndofs,1);
% Incidence matrix
Incid = zeros(Nel,Elnodes);
i = 1;
for e = 1:Nel
for j = 2:Elnodes
Incid(e,1) = (e-1)*(Elnodes-1)+1;
Incid(e,j) = (i+1);
i = i+1;
end
end
% Nodes coordinates
X = linspace(0,Lb,Nvert);
% Jacobi root and weight
[r,w] = GetPointsWeight1D(0,0,10,’GJ’,’NL’);
if intorder==1 % Linear interpolation
tipo = ’lineares’;
% dN = [-0.5+0*r 0.5+0*r];
4
% N = [0.5*(1-r) 0.5*(1+r)];
elseif intorder==2 % Quadratic interpolation
tipo = ’Quadraticos’;
% dN = [-0.5+r -2*r 0.5+r];
% N = [-0.5*r.*(1-r) 1-r.ˆ2 0.5*r.*(1+r)];
end
N = GetPolynomialValue(intorder,r,linspace(-1,1,Elnodes))’;
dN = Get1stDerivative(intorder,r,linspace(-1,1,Elnodes))’;
% Assembly of global matrix and global load vector without b.c.’s
for e = 1 : Nel
xl = X(Incid(e,:))’; % Element coordinates (local)
he = xl(end)-xl(1); % Element lenght
dJac = he/2; % Determinant of the Jacobian
dJinv = 1/dJac; % Jacobian inverse
%Stiffness matrix
Ke = zeros(Elnodes,Elnodes);
for i = 1:Elnodes
for j = 1:Elnodes
for k = 1:length(w)
Ke(i,j) = Ke(i,j) + EA*w(k)*dN(k,i)*dN(k,j)*dJinv;
end
end
end
Kg(Incid(e,:),Incid(e,:)) = Kg(Incid(e,:),Incid(e,:)) + Ke;
% Load vector
Fe = zeros(Elnodes,1);
for i = 1:Elnodes
for k = 1:length(w)
Fe(i) = Fe(i) + (pi/Lb)*sin(pi*(N(k,:)*xl)/Lb)*N(k,i)*w(k)*dJac;
end
end
Fg(Incid(e,:)) = Fg(Incid(e,:)) + Fe;
end
% Application of b.c. disp. condition at x = 0 and x = L
% elimination of the first/last line and column of the global system
ncc = 2:(Ndofs-1);
ug(ncc,1) = Kg(ncc,ncc)\Fg(ncc);
eh = sqrt(ug’*Kg*ug);
syms x real
ua = (Lb/(pi*EA))*sin(pi*x/Lb);
dua = (1/EA)*cos(pi*x/Lb);
ea = sqrt(eval(int(EA*duaˆ2,x,0,Lb)));
% Relative error
erro = (ea - eh)/ea;
fprintf(’O erro relativo utilizando %i elementos %s é: %e\n’,Nel,tipo,erro)
5