Uma Matriz de Rigidez Espectral
Marco L. Bittencourt,
Thais G. Vazquez,
Depto de Projeto Mecânico, FEM, UNICAMP,
13083-970, Campinas, SP
E-mail: [email protected], [email protected].
Neste trabalho novas funções de forma 1D
são propostas para a versão p do método dos
elementos finitos (MEF) com o objetivo de encontrar uma matriz de rigidez espectral. A
função de forma proposta fornece uma matriz
de rigidez unidimensional quase diagonal e que
apresenta melhor condicionamento do que a
matriz de rigidez quando usada as bases de Lagrange ou Jacobi.
1
Introdução
A escolha de um conjunto de funções de forma
para a solução de Problemas de Valor de Contorno (PVC) usando MEF deve ser analisada
quanto à sua eficiência numérica, condicionamento e independência linear, tal que forme
uma base para o espaço das funções propostas.
Um problema de aproximação pode se tornar
equivalente à resolver uma equação matricial da forma M û = f , a qual é resolvida
invertendo-se a matriz M , tal que û = M −1 f ,
onde M é uma matriz, û o vetor solução e f
um carregamento [4].
Porém, o custo computacional para inverter
sistemas matriciais é muito alto. Portanto,
quanto mais próxima de uma matriz diagonal
for possı́vel representar o problema, mais fácil
será a solução do problema inicial. Também,
o condicionamento numérico da matriz M está
relacionado com a independência linear da expansão.
Quando invertemos numericamente um sistema matricial existe um erro associado com
a representação exata da matriz, devido ao
arredondamento. Se a matriz for mal condicionada, o erro de arredondamento pode conduzir a erros muito grandes na solução [4].
Um conjunto de funções de forma espectrais
bastante citados na literatura são baseados nos
polinômios de Lagrange de grau P colocados
em (P + 1) pontos nodais [4]. Esses polinômios
obedecem à propriedade de colocação, em que
os coeficientes representam a aproximação da
solução nos pontos de colocação. Os pontos de
colocação são chamados de nós e portanto a expansão usando base de Lagrange é referenciada
como uma expansão nodal.
Existem também expansões modais, como
por exemplo um conjunto de funções usando
polinômios ortogonais de Jacobi. Com a escolha ótima das ponderação desses polinômios é
possı́vel melhorar a esparsidade e o condicionamento das matrizes de massa e rigidez [4].
2
Integração
Numérica
Quadratura Gaussiana
e
A integral
da função u(ξ), ξ ∈ [−1, 1], dada
R1
por −1
u(ξ)dξ, pode ser resolvida usando
integração numérica baseadas nas regras de
quadraturas dadas por [4]
Z
1
u(ξ)dξ =
−1
Q−1
X
wi u(ξi ),
i=0
onde wi são constantes especı́ficas ou pesos e ξi
são Q pontos distintos no intervalo [−1, 1].
Na quadratura gaussiana, o integrando u(ξ)
é aproximado por um polinômio de Lagrange
usando Q pontos ξi , da seguinte maneira [4]
u(ξ) =
Q−1
X
u(ξi )hi (ξ) + ǫ(u),
i=0
onde
hi (ξ)
=
Q−1
Y
(ξ − ξq )
(ξ
p − ξq )
q=0,q6=i
são
os
polinômios de Lagrange e ǫ(u) é o erro da
aproximação.
Assim
Z
1
u(ξ)dξ =
−1
Q−1
X
i=0
wi u(ξi ) + R(u),
onde
A matriz de massa é cheia se avaliarmos o
produto
interno exato. Se, contudo, usarmos a
1
1
wi =
hi (ξ)dξ and R(u) =
ǫ(u)dξ.
regra de quadratura de Gauss-Lobatto-Jacobi
−1
−1
nos pontos de integração ξi iguais aos pontos
Como a equação para wi define a ex- de colocação no qual a expansão foi definida,
pressão para o cálculo dos pesos em termos a matriz de massa torna-se diagonal, ou seja,
da integração dos polinômios de Lagrange, é obedece a propriedade de colocação [4]
necessário conhecer os pontos ξi . Como u(ξ) é
P
P
X
X
representado por um polinômio de grau Q−1, a
e
M
[p][q]
=
w
N
(ξ
)N
(ξ
)
=
wi δpi δqi
i p i
q i
integração será exata se u(ξ) for um polinômio
i=0
i=0
de grau menor ou igual a Q − 1. Assim, ao es= wp δpq .
(3)
colher os pontos ξi , podemos escolher os pontos igualmente espaçados no intervalo [−1, 1] ou
Pode ser observado que R(u) = 0 em (2) se
outros pontos que forneçam a integração exata
u(ξ)
∈ P2P −3 ([−1, 1]). Se wi são os pesos da
para polinômios de ordem maior que Q − 1.
regra
de quadratura para (P + 2) pontos, então
Existem diferentes tipos de regras de
quadratura de Gauss.
Aqui, usaremos a a integração é exata, mas a ortogonalidade não
quadratura de Gauss-Lobatto-Jacobi, a qual é preservada. O erro cometido quando usada
usa os polinômios de Lagrange de grau P colo- essa redução na ordem da regra de quadratura
é consistente com o erro da aproximação da excados nas seguintes (P + 1) raı́zes [4]
pansão, isto é, o erro na aproximação da função

usando a expansão polinomial é da mesma ori=0

 −1
α+1,β+1
(1) dem que o erro da integração Gaussiana com
ξi =
ξi−1,P −1 i = 1, 2, . . . , P − 1.


um ponto a menos do que o necessário para a
1
i=P
integração exata [4]. Os elementos da diagonal
α+1,β+1
Os pontos ξi−1,P −1 são os zeros dos polinômios da matriz de massa espectral usando a regra de
de Jacobi cim pesos α + 1 e β + 1, denotado por quadratura são iguais a soma dos elementos das
linhas respectivas da matriz de massa usando
PPα+1,β+1
(ξ).
−1
A regra de quadratura de Gauss-Lobatto- integração exata. A expressão (3) representa
uma ortogonalidade discreta das funções de inJacobi é, então, dada por [4]
terpolação nodais.
Z
Z
1
−1
Z
α
β
(1−ξ) (1+ξ) u(ξ)dξ =
P
X
wiα,β u(ξi )+R(u),
i=0
(2) 4 Matriz de Rigidez Espectral
onde u(ξ) é uma função polinomial, R(u) é o
Unidimensional
erro cometido na aproximação e os pesos corOs polinômios de Jacobi unidimensionais de
respondentes são
grau n, Pnα,β (ξ), são uma famı́lia de soluções

α,β

(β
+
1)C
i
=
0

polinomiais para o problema singular de Sturm0,P −1

α,β
wiα,β =
,
Ci,P
i
=
1,
2,
.
.
.
,
P
−
1
Liouville [4]. Estes polinômios são ortogonais
−1


 (α + 1)C α,β
[2] no intervalo [−1, 1] com respeito à função
P,P −1 i = P
peso (1 − ξ)α (1 + ξ)β , (α, β > −1; α, β ∈ ℜ), da
seguinte maneira
2α+β+1 (α + P )!(β + P )!
α,β
Ci,P
.
−1 =
P P !(α + β + P + 1)![PPα,β (ξi )]2
Z 1
α,β
(1 − ξ)α (1 + ξ)β Pnα,β (ξ)Pm
(ξ)dξ = Cδnm ,
3
Matriz de Massa Espectral
Unidimensional
−1
(4)
ξ ∈ [−1, 1] e
Conforme descrito em [3, 4], uma classe de ele2α+β+1
(n + α)!(n + β)!
C
=
mentos nodais, conhecidos como elementos es2n + α + β + 1 n!(n + α + β)!
pectrais usam os polinômios de Lagrange de
grau P colocados nos pontos ξi dados por (1). uma constante.
As funções de interpolação modais podem ser
definidas por [4, 1]
Np (ξ)
=
φp (ξ)
 1
 2 (1 − ξ)
1
(1 + ξ)
 21
α,β
4 (1 − ξ)(1 + ξ)Pp−1 (ξ)
=
p=0
p=P
0<p<P
Observe que as funções de vértice são as mesmas funções do elemento linear de Lagrange.
Os coeficientes das matrizes de rigidez dos
elementos unidimensionais são dados por
Z
e
K [p][q] =
1
Np′ (ξ)Nq′ (ξ)dξ,
−1
com 0 ≤ p, q ≤ P . Os coeficientes correspondentes aos modos internos podem ser escritos,
após a integração por partes, como [4]
K e [p][q] =
−
1
4
Z
1
α,β
(1 − ξ)α (1 + ξ)β Pp−1
(ξ)
−1
d2
φq (ξ)dξ.
d2 ξ
Usando a relação de ortogonalidade dos
polinômios de Jacobi (4), com α = β = 1 e
d2
observando que 2 φq (ξ) é um polinômio de
d ξ
grau q − 1, os coeficientes internos da matriz
de rigidez zeram para p > q, e analogamente
p < q. Portanto, a matriz relacionada aos modos internos é diagonal [4].
Para construir uma matriz de rigidez spectral nodal similar à matriz de massa, o seguinte
conjunto de funções é definido


1
2 (1
− ξ),
(1 − ξ)φp (ξ),
Np (ξ) =
 1
2 (1 + ξ),
p=0
1 ≤ p ≤ P + 1 , (5)
p=P +2
com
φ1 (ξ)/8 =
Essa escolha garante que as funções de
vértice são lineares e as funções internas zeram
para ξ0 = −1 e ξP = 1. Em adição, a C 0 . continuidade entre os elementos é preservada
[5].
Observe que, por exemplo, para P = 4, serão
necessárias 7 funções de forma para completar
o espaço, já que o grau das funções é 6. Assim,
teremos sempre P + 3 funções de interpolação.
A constante 8 foi usada apenas para melhorar o condicionamento numérico da matriz de
rigidez.
As funções φp (ξ) são tais que as derivadas
′
φp (ξ) das funções internas tenham a propriedade de colocação. Para exemplificar, considere P = 4 e os pontos de colocação ξ0 = −1,
ξ1 , ξ2 , ξ3 e ξ4 = 1. As expressões das derivadas
φ′p (ξ) são
φ′1 (ξ) = 2(ξ − ξ4 )(ξ − ξ1 )(ξ − ξ2 )(ξ − ξ3 ),
φ′2 (ξ) = 2(ξ − ξ0 )(ξ − ξ2 )(ξ − ξ3 )(ξ − ξ4 ),
φ′3 (ξ) = 2(ξ − ξ0 )(ξ − ξ1 )(ξ − ξ3 )(ξ − ξ4 ),
φ′4 (ξ) = 2(ξ − ξ0 )(ξ − ξ1 )(ξ − ξ2 )(ξ − ξ4 ),
φ′5 (ξ) = 2(ξ − ξ0 )(ξ − ξ1 )(ξ − ξ2 )(ξ − ξ3 ).
As expressões anteriores correspondem à
polinômios de Lagrange, a menos de uma constante.
Para o bloco das funções internas na matriz
de rigidez, a expressão dos coeficientes será
e
K [p][q] =
= 64
Z
Z
1
Np′ (ξ)Nq′ (ξ)dξ
−1
1
[(1 − ξ)φ′p (ξ) − φp (ξ)][(1 − ξ)φ′q (ξ) − φq (ξ)]dξ.
−1
Após a integração por partes concluı́mos que
2
(ξ − ξP ) (ξ − ξ1 )(ξ − ξ2 ) . . . (ξ − ξP −1 )
P +1
+
X
(−1)
i 2
i!
(ξ − ξP )
id
e
K [p][q] = 64
i−2
dξ
((ξ − ξ1 ) . . . (ξ − ξP −1 ))
i=3
2
−(−1 − ξP ) (−1 − ξ1 )(−1 − ξ2 ) . . . (−1 − ξP −1 )
P +1
−
X
(−1)
i 2
i!
(−1 − ξP )
id
i−2
dξ
((ξ − ξ1 ) . . . (ξ − ξP −1 ))ξ=−1
i=3
e, para 2 < p < P + 1,
φp (ξ)/8 =
(ξ − ξ0 )
+
P +1
X
i=3
1
−1
(1 − ξ)2 φ′p (ξ)φ′q (ξ)dξ.(6)
Se escolhermos como pontos de colocação e
integração os P + 1 zeros de Gauss-LobattoJacobi para α = 2 e β = 0, a integral anterior
pode ser aproximada por
K e [p][q] = 64
P
X
φ′p (ξi2,0 )φ′q (ξi2,0 )wi2,0 = 64wi2,0 ,
i=0
2 (ξ − ξ1 )(ξ − ξ2 ) . . . (ξ − ξP )
(−1)
Z
(ξ − ξp−1 )
i 2
i!
(ξ − ξ1 )
id
i−2
dξ
(ξ − ξ1 ) . . . (ξ − ξP )
(ξ − ξp−1 )
.
onde a propriedade de colocação das derivadas
φ′p (ξ) foram usadas. Baseado nisso, o bloco interno da matriz de rigidez será diagonal. O
erro cometido na integração será da mesma ordem que o erro da aproximação da função polinomial. Observe que o termo (1 − ξ)2 corresponde à função de ponderação da quadratura
(2) e não precisa ser integrado explicitamente
tomando-se α = 2 e β = 0.
1
2
3
4
5
6
7
8
9
0
1
2
3
4
5
6
7
8
9
nz = 10
Os coeficientes do bloco das funções de
1
vértice da matriz de rigidez são iguais a ± . Figura 1: Perfil de esparsidade para P = 5 (nz
2 é o número de coeficientes diferentes de zero).
Já os coeficientes do bloco de acoplamento
vértice-interno zeram naturalmente quando inde colocação são, respectivamente,
tegramos por partes, pois
Z
1
−1
1
±
2
5






K=





Np′ (ξ)Nq′ (ξ)dξ =
Z
1
−1
[(1 − ξ)φ′q (ξ) − φq (ξ)]dξ = 0.
Considere o problema de Poisson unidimenKap
sional
d2 φ(ξ)
= ξ P −2 ,
dξ 2
ξ ∈ [−1, 1],
P = 2, 3, . . . , 8,
(7)
com condições de contorno homogêneas
i. φ(−1) = φ(1) = 0;
ii.
dφ
dξ (−1)
= φ(1) = 0;
iii. φ(−1) =
dφ
dξ (1)
= 0.
O problema anterior foi resolvido usando
polinômios de Lagrange integrados analiticamente e as funções propostas integrados analiticamente e por colocação com (P + 1) pontos de integração. As respectivas matrizes de
rigidez são denotadas por Klag , K e Kap , respectivamente.
Os perfis de esparsidade da matriz de rigidez
para as funções de forma propostas, integrada
nos pontos de colocação, tem a mesma caracterı́stica (para qualquer grau P ), como
ilustrado na Figura 1 para P = 5. Esse é
o mesmo perfil de esparsidade da matriz de
rigidez unidimensional empregando polinômios
de Jacobi [4].
As matrizes de rigidez para P = 4, integradas analiticamente e nos (P + 1) pontos
0.0000
0.0000
0.0000
0.0000
0.0000
−0.5
0.0
18.9137
−1.8914
−1.8914
−1.8914
−1.8914
0.0
0.0
−1.8914
18.9137
−1.8914
−1.8914
−1.8914
0.0
0.0
−1.8914
−1.8914
18.9137
−1.8914
−1.8914
0.0
0.0
−1.8914
−1.8914
−1.8914
18.9137
−1.8914
0.0
0.0
−1.8914
−1.8914
−1.8914
−1.8914
5.0437
0.0
−0.5
0.0000
0.0000
0.0000
0.0000
0.0000
0.5

Resultados
0.5





=





0.5
0.0000
0.0000
0.0000
0.0000
0.0000
−0.5
0.0
20.8051
0.0000
0.0000
0.0000
0.0000
0.0
0.0
0.0000
20.8051
0.0000
0.0000
0.0000
0.0
0.0
0.0000
0.0000
20.8051
0.0000
0.0000
0.0
0.0
0.0000
0.0000
0.0000
20.8051
0.0000
0.0
0.0
0.0000
0.0000
0.0000
0.0000
6.9350
0.0
−0.5
0.0000
0.0000
0.0000
0.0000
0.0000
0.5
Observe que os coeficientes fora da diagonal
principal relacionados às funções internas da
matriz K são iguais. Nesse caso, esse valor é
k̄ = −1.8914. Pode ser verificado que os coeficientes do bloco interno da matriz Kap são
dados pela seguinte expressão
Kapii = Kii + |k̄|.
O erro na norma de energia entre a solução
exata u e aproximada uap calculado para cada
grau P = 2 até P = 8 está mostrado na Figura
2. Observe que para P = 4 o erro será exatamente igual a zero.
Os números de condição das matrizes de
µ
rigidez são calculados por k = max
min µ , sendo µ
o conjunto de valores singulares das matrizes.
Como a matriz de rigidez é semi-definida positiva, min µ é o primeiro valor singular diferente
de zero.
Os condicionamentos numéricos das matrizes
de rigidez obtidas com as funções propostas,
integradas analiticamente e por colocação,
polinômios de Lagrange e Jacobi com α = β =
1, são mostrados na Figura 3. Observa-se que






,











.





(9)
1.4
1.2
A propriedade de colocação será preservada
somente quando p 6= r e q 6= s simultaneamente. Isso acontece somente para alguns coeficientes do bloco interno da matriz de rigidez,
como ilustrado na Figura 4(a) para P = 4. Esse
perfil pode ser comparado com a esparsidade
do bloco interno da matriz de rigidez obtida
Figura 2: Erro na norma de Energia para 2 ≤ usando polinômios de Jacobi com α = β = 1
P ≤ 8.
na Figura 4(b).
Error
1
0.8
0.6
0.4
0.2
0
2
3
4
5
6
7
8
Degree
0
as matrizes obtidas com as funções propostas
apresentam um melhor condicionamento com
relaçao às outras bases.
5
10
One−dimensional Stiffness
4
15
10
KJacobi
KLagrange
K
Kap
20
3
10
Condition Number
25
0
5
10
15
20
25
nz = 481
(a) Funções propostas.
2
10
0
1
10
5
10
0
10
2
3
4
5
6
7
8
Polynomial degree
15
Figura 3: Condicionamento das matrizes de
rigidez usando polinômios de Jacobi e Lagrange, calculadas analiticamebte e aproximadas para 2 ≤ P ≤ 8.
O resultado anterior não se extendeu para o
bloco das funções internas da matriz de rigidez
do caso bidimensional (quadrado).
As funções 2D, Np (ξ, η) e Nq (ξ, η) são obtidas pelo produto tensorial usual de funções
unidimensionais, i.e. Np (ξ, η) = Ni (ξ)Nj (η)
e Nq (ξ, η) = Nk (ξ)Nl (η). Portanto, os coeficientes do bloco interno da matriz de rigidez
para o problema de Poisson nos elementos
quadrangulares são dados por
20
25
0
5
10
15
20
25
nz = 189
(b) Jacobi com α = 1, β =
1.
Figura 4: Perfis de esparsidade da matriz de
rigidez para quadrados e grau P = 4 (nz é o
número de coeficientes diferentes de zero).
6
Conclusões
O objetivo inicial deste trabalho era encontrar
funções de forma que fornecessem, ao mesmo
Z 1Z 1
e
K [p][q] =
∇Np (ξ, η) · ∇Nq (ξ, η)dξdη tempo, as matrizes de massa e rigidez espec−1 −1
trais. A primeira idéia foi dobrar o grau de toZ 1Z 1
dNi (ξ) dNk (ξ)
dos os monômios do polinômio de Lagrange e,
=
Nj (η)Nl (η)
dξ
dξ
−1 −1
desta forma, ter o polinômio e sua derivada nas
dNj (η) dNl (η)
mesmas raı́zes. Contudo, o problema deste pro+Ni (ξ)Nk (ξ)
dξdη
dη
dη
cedimento foi aumentar muito o grau da função
Z 1
Z 1
dNi (ξ) dNk (ξ)
de forma e não ter pontos suficientes para a in=
Nj (η)Nl (η)dη
dξ
dξ
dξ
tegração numérica por colocação.
−1
−1
Z 1
Z 1
Além disso, ainda precisávamos garantir que
dNj (η) dNl (η)
+
Ni (ξ)Nk (ξ)dξ
dη.
(8) as funções internas fossem zero nos dois exdη
dη
−1
−1
tremos para podermos eliminar o termo do contorno quando se aplica a forma fraca na matriz de rigidez, e os vértices respectivos fossem unitário no seu ponto de colocação (o
primeiro vértice igual a 1 no primeiro ponto
de integração e o último vértice igual a 1 no
último ponto de integração), para garantir a
continuidade C 0 dos elementos.
De acordo com a equação (8), como a matriz de massa unidimensional não é diagonal
e a matriz de rigidez unidimensional tem dois
coeficientes diferentes de zero no bloco vérticeinterna o processo de tensorização usado para
obter a matriz de rigidez bidimensional fornece
uma matriz de rigidez não-diagonal.
Portanto, a diagonalização da matriz de
rigidez espectral para elementos 2D e 3D requer
a propriedade de colocação das funções unidimensionais e as suas respectivas derivadas.
Referências
[1] M. L. Bittencourt, Fully tensorial nodal
and modal shape functions for triangles
and tetrahedra, Int. J. for Num. Meth. in
Eng.., 63 (2005) 1530-1558.
[2] Theodore S. Chihara, An Introduction
to Orthogonal Polynomials, Mathematics
and its Applications Series, Gordon and
Breach - Science Publishers., (1978), New
York.
[3] A.T. Patera, A spectral method for fluid
dynamics: Laminar flow in a channel expasion, Journal of Computational
Physics., 54 (1984) 468.
[4] G. E. Karniadakis e S. J. Sherwin, Spectral/hp Element Methods for CFD, Oxford
University Press., (1999).
[5] O. C. Zienkiewicz and R. L. Taylor, The
Finite Element Method, McGraw-Hill International Editions., (1989).