Método direto da rigidez

O método direto da rigidez, também denominado método matricial da rigidez e método dos deslocamentos (em inglês direct stiffness method (DSM)) é um método de cálculo aplicável a estruturas hiperestáticas de barras que se comportam de forma elástica linear. Este método é adequado para realizar análise computacional de qualquer estrutura, incluindo estruturas estaticamente indeterminadas. O método matricial se baseia em estimar as componentes das relações de rigidez para resolver as forças ou os deslocamentos com uso de um computador. O método direto da rigidez é a implementação mais comum do Método dos elementos finitos (MEF). As propriedades de rigidez do material são computadas em uma única equação matricial que governa o comportamento interno da estrutura idealizada. Os dados que se desconhecem na estrutura são as forças e deslocamentos, que podem ser determinados resolvendo esta equação. O método direto da rigidez é o mais comum nos programas de cálculo de estruturas (tanto comerciais quanto livres).

O método foi originado na área da aeronáutica. Os investigadores conseguiram aproximar o comportamento estrutural das partes de um avião mediante equações simples que requeriam porém grandes tempos de computação. Com o advento dos computadores estas equações começaram a ser resolvidas de forma rápida e simples.

Introdução

O método consiste em atribuir à estrutura de barras um objeto matemático chamado matriz de rigidez, que relaciona os deslocamentos de um conjunto de pontos da estrutura, chamados nós, com as forças exteriores que é necessário aplicar para obter estes deslocamentos (as componentes desta matriz são forças generalizadas associadas a deslocamentos generalizados). A matriz de rigidez relaciona as forças nodais equivalentes e deslocamentos sobre os nó da estrutura, mediante a equação

${\displaystyle {\begin{Bmatrix}F_{1}+R_{1}\\F_{2}+R_{2}\\...\\F_{n}+R_{n}\end{Bmatrix}}_{G}={\begin{bmatrix}k_{11}&k_{12}&\cdots &k_{1n}\\k_{21}&k_{22}&\cdots &k_{2n}\\\cdots &\cdots &\ddots &\cdots \\k_{n1}&k_{n2}&\cdots &k_{nn}\end{bmatrix}}_{G}{\begin{Bmatrix}\delta _{1}\\\delta _{2}\\...\\\delta _{n}\end{Bmatrix}}_{G}~~.}$

Nesta equação ${\displaystyle F_{i}\,}$ são as forças nodais equivalentes associadas às forças exteriores aplicadas sobre a estrutura; ${\displaystyle R_{i}\,}$ são as reações hiperestáticas inicialmente desconhecidas sobre a estrutura; ${\displaystyle \delta _{i}\,}$ os deslocamentos nodais incógnitos da estrutura e ${\displaystyle n\,}$ é o número de graus de liberdade da estrutura.

A energia de deformação elástica também pode ser expressa em termos da matriz de rigidez mediante a relação

${\displaystyle E_{def}={\frac {1}{2}}{\boldsymbol {\delta }}\cdot [\mathbf {K} ({\boldsymbol {\delta }})]={\frac {1}{2}}\sum _{i,j}k_{ij}\delta _{i}\delta _{j}~~.}$

Do teorema de Maxwell-Betti deduz-se que a matriz de rigidez deve ser simétrica, e portanto

${\displaystyle k_{ij}=k_{ji}~~.}$

Fundamentação teórica

Em geral, um sólido deformável real, como qualquer meio contínuo, é um sistema físico com um número infinito de graus de liberdade. Assim, para descrever a deformação de um sólido é necessário estabelecer um campo vetorial de deslocamentos sobre cada um de seus pontos. Este campo de deslocamentos em geral não é redutível a um número finito de parâmetros, e portanto um sólido deformável de forma totalmente geral não tem um número finito de graus de liberdade.

No entanto, para barras elásticas longas (prismas mecânicos) de grande comprimento comparado com a área de sua seção transversal, o campo de deslocamentos é descrito pela curva elástica, cuja deformação pode ser reduzida a um conjunto finito de parâmetros. De forma concreta, fixados os deslocamentos e rotações das seções extremas de uma viga elástica, é completamente determinada sua forma. Assim, para uma estrutura formada por barras longas elásticas, fixados os deslocamentos dos nós, é completamente determinada sua forma deformada. Resulta assim que as estruturas de barras longas podem ser tratadas aproximadamente mediante um número finito de graus de liberdade e que podem ser calculadas resolvendo um número finito de equações algébricas. O método matricial proporciona estas equações em forma de um sistema matricial que relaciona os deslocamentos das extremidades de barras com variáveis dependentes das forças exteriores.

Isto contrasta com a situação geral dos sólidos elásticos, onde o cálculo de suas tensões internas e deformações envolve a resolução de complexos sistemas de equações diferenciais em derivadas parciais.

Descrição do método

O método matricial requer atribuir a cada barra elástica da estrutura uma matriz de rigidez, chamada matriz de rigidez do elemento. A partir do conjunto de matrizes de cada elemento, mediante uma matriz de conectividade, que estabelece a forma como cada barra é conectada com outra através dos nós, é obtida uma matriz de rigidez global, que relaciona os deslocamentos dos nós com as forças equivalentes sobre os mesmos.

Igualmente, a partir das forças aplicadas sobre cada barra se constrói o chamado vetor de forças nodais equivalentes que depende das ações externas sobre a estrutura. Junto com estas forças devem ser consideradas as possíveis reações sobre a estrutura em seus apoios (cujos valores são incógnitos).

Finalmente se constrói um sistema linear de equações, para os deslocamentos e as incógnitas. O número de reações incógnitas e deslocamentos incógnitos depende do número de nós: é igual a 3N para problemas bidimensionais, e igual a 6N para um problema tridimensional. Este sistema pode ser sempre dividido em dois subsistemas de equações desacopladas tal que

• Subsistema 1. Que agrupa todas as equações lineares do sistema original que contém somente deslocamentos incógnitos.
• Subsistema 2. Que agrupa o resto das equações, e que uma vez resolvido o subsistema 1 e substituídos seus valores no subsistema 2 permite encontrar os valores das reações incógnitas.

Resolvido o subsistema 1 que dá os deslocamento, é substituído o valor destes no subsistema 2, que é trivial de resolver. Finalmente, a partir das reações, forças nodais equivalentes e deslocamentos são calculados os esforço nos nós ou uniões das barras, a partir dos quais podem ser determinados os esforços em qualquer ponto da estrutura e portanto suas tensões máximas, que permitem dimensionar adequadamente todas as seções da estrutura.

Matrizes de rigidez elementares

Para construir a matriz de rigidez da estrutura é necessário atribuir previamente a cada barra individual (elemento) uma matriz de rigidez elementar. Esta matriz depende exclusivamente de

1. as condições de ligação em seus dois nós extremos
2. as características da seção transversal da barra: área, momento de inércia de área e características geométricas gerais como comprimento da barra, curvatura, etc.
3. número de graus de liberdade por nó, que depende de tratar-se de problemas bidimensionais ou tridimensionais.

Barra reta bidimensional com transmissão de momentos

Um nó onde se unem duas barras é denominado rígido ou engastado se o ângulo formado pelas duas barras após a deformação não muda em relação ao ângulo que formavam antes da deformação. Apesar de estarem impossibilitadas de mudar o ângulo entre as barras, as duas barras em conjunto podem rotacionar em relação ao nó, mantendo porém o ângulo formado. Para barras unidas rigidamente em seus dois extremos a matriz de rigidez elementar que representa adequadamente seu comportamento é expressa por

${\displaystyle \left[K^{(e)}\right]={\begin{bmatrix}{\frac {EA}{L}}&0&0&-{\frac {EA}{L}}&0&0\\0&{\frac {12EI}{L^{3}}}&{\frac {6EI}{L^{2}}}&0&-{\frac {12EI}{L^{3}}}&{\frac {6EI}{L^{2}}}\\0&{\frac {6EI}{L^{2}}}&{\frac {4EI}{L}}&0&-{\frac {6EI}{L^{2}}}&{\frac {2EI}{L}}\\-{\frac {EA}{L}}&0&0&{\frac {EA}{L}}&0&0\\0&-{\frac {12EI}{L^{3}}}&-{\frac {6EI}{L^{2}}}&0&{\frac {12EI}{L^{3}}}&-{\frac {6EI}{L^{2}}}\\0&{\frac {6EI}{L^{2}}}&{\frac {2EI}{L}}&0&-{\frac {6EI}{L^{2}}}&{\frac {4EI}{L}}\end{bmatrix}}~~.}$

Sendo aqui ${\displaystyle L,A,I,\,}$ as magnitudes geométricas (comprimento, área e momento de inércia) e ${\displaystyle E\,}$ o módulo de de elasticidade longitudinal (módulo de Young).

A matriz de rigidez de uma barra reta pode ser expressa abreviadamente na forma

${\displaystyle \left[K^{(e)}\right]={\frac {EI}{L^{3}}}{\begin{bmatrix}\lambda _{k}^{2}&0&0&-\lambda _{k}^{2}&0&0\\0&12&6L&0&-12&6L\\0&6L&4L^{2}&0&-6L&2L^{2}\\-\lambda _{k}^{2}&0&0&\lambda _{k}^{2}&0&0\\0&-12&-6L&0&12&-6L\\0&6L&2L^{2}&0&-6L&4L^{2}\end{bmatrix}}~~,}$

sendo ${\displaystyle \lambda _{k}:={\sqrt {\frac {AL^{2}}{I}}}={\frac {L}{i_{giro}}}}$ o índice de esbeltez.

Barra reta bidimensional com um nó articulado e outro rígido

Neste caso quando se impõem rotações no nó articulado não são transmitidos esforços sobre o nó não articulado. Neste caso a matriz de rigidez, usando a mesma notação da seção anterior, é expressa por

${\displaystyle \left[K^{(e)}\right]={\frac {EI}{L^{3}}}{\begin{bmatrix}\lambda _{k}^{2}&0&0&-\lambda _{k}^{2}&0&0\\0&3&3L&0&-3&0\\0&3L&3L^{2}&0&-3L&0\\-\lambda _{k}^{2}&0&0&\lambda _{k}^{2}&0&0\\0&-3&-3L&0&3&0\\0&0&0&0&0&0\end{bmatrix}}}$

Onde foi suposto que o nó articulado é o segundo. Se fosse o primeiro, seria necessário permutar os elementos da matriz anterior para obter

${\displaystyle \left[K^{(e)}\right]={\frac {EI}{L^{3}}}{\begin{bmatrix}\lambda _{k}^{2}&0&0&-\lambda _{k}^{2}&0&0\\0&3&0&0&-3&3L\\0&0&0&0&0&0\\-\lambda _{k}^{2}&0&0&\lambda _{k}^{2}&0&0\\0&-3&0&0&3&-3L\\0&3L&0&0&-3L&3L^{2}\end{bmatrix}}}$

Barra reta bidimensional com dois nós articulados

Posto que uma barra reta com nós articulados pode transmitir esforços somente ao longo de seu eixo, a correspondente matriz de rigidez desta barra tem componentes diferentes de zero apenas para os graus de liberdade longitudinais. Neste caso a matriz de rigidez é expressa por

${\displaystyle \left[K^{(e)}\right]={\frac {EA}{L}}{\begin{bmatrix}+1&0&0&-1&0&0\\0&0&0&0&0&0\\0&0&0&0&0&0\\-1&0&0&+1&0&0\\0&0&0&0&0&0\\0&0&0&0&0&0\end{bmatrix}}~~.}$

Arco circular bidimensional de nós engastados

Barra reta tridimensional de nós rígidos

Uma barra reta tridimensional tem 6 graus de liberdade por nó (3 de traslação e 3 de rotação), e como a barra tem dois nós, a matriz de rigidez é uma matriz 12 x 12. Ademais, uma barra tridimensional pode transmitir torção e também flexão e esforço cortante em duas direções diferentes, esta maior complexidade de comportamento estrutural faz com que uma barra tridimensional requeira mais graus de liberdade e uma matriz de rigidez mais complexa para descrever seu comportamento. Esta matriz é composta de 3 submatrizes

${\displaystyle \left[K^{(e)}\right]={\begin{bmatrix}\mathbf {A} _{1}&\mathbf {B} ^{T}\\\mathbf {B} &\mathbf {A} _{2}\end{bmatrix}}~~.}$

As submatrizes são

${\displaystyle \left[\mathbf {A} _{i}\right]={\begin{bmatrix}{\frac {EA}{L}}&0&0&0&0&0\\0&{\frac {12EI_{z}}{L^{3}}}&0&0&0&\epsilon _{i}{\frac {6EI_{z}}{L^{2}}}\\0&0&{\frac {12EI_{y}}{L^{3}}}&0&-\epsilon _{i}{\frac {6EI_{y}}{L^{2}}}&0\\0&0&0&{\frac {GJ}{L}}&0&0\\0&0&-\epsilon _{i}{\frac {6EI_{y}}{L^{2}}}&0&{\frac {4EI_{y}}{L}}&0\\0&\epsilon _{i}{\frac {6EI_{z}}{L^{2}}}&0&0&0&{\frac {4EI_{z}}{L}}\\\end{bmatrix}}}$

${\displaystyle \left[\mathbf {B} \right]={\begin{bmatrix}-{\frac {EA}{L}}&0&0&0&0&0\\0&-{\frac {12EI_{z}}{L^{3}}}&0&0&0&-{\frac {6EI_{z}}{L^{2}}}\\0&0&-{\frac {12EI_{y}}{L^{3}}}&0&+{\frac {6EI_{y}}{L^{2}}}&0\\0&0&0&-{\frac {GJ}{L}}&0&0\\0&0&-{\frac {6EI_{y}}{L^{2}}}&0&{\frac {2EI_{y}}{L}}&0\\0&{\frac {6EI_{z}}{L^{2}}}&0&0&0&{\frac {2EI_{z}}{L}}\\\end{bmatrix}}}$

As magnitudes geométricas e mecânicas associadas à barra são

${\displaystyle L,A;I_{y},I_{z};J\,}$ são as magnitudes geométricas: comprimento da barra e sua área transversal, momentos de área nas direções y e z e módulo de torção, respectivamente.

${\displaystyle E,G\,}$ o módulo de Young e o módulo de cisalhamento.

${\displaystyle \epsilon _{1}=+1,\epsilon _{2}=-1\,}$ são sinais relativos.

Forças nodais

Para cada barra se define um vetor elementar de forças nodais generalizadas, que seja estaticamente equivalente às forças aplicadas sobre a barra. A dimensão do vetor de forças nodais depende da dimensionalidade da barra:

${\displaystyle \{\mathbf {F} ^{(e)}\}\in {\begin{cases}\mathbb {R} ^{6}&{\mbox{bidimensional}}\\\mathbb {R} ^{12}&{\mbox{tridimensional}}\end{cases}}}$

Exemplo

Para as cargas mostradas na figura o vetor de forças nodais consiste em duas forças verticais (F_Vd, F_Vi) aplicadas em cada um dos dois extremos, duas forças horizontais (F_Hd, F_Hi) aplicadas em cada um dos extremos e dois momentos (M_d, M_i) aplicados em cada um dos extremos. Estas seis componentes formam o vetor de forças nodais. A força e o momento resultante destas seis componentes são estaticamente equivalentes ao sistema de forças original formado por P e q com

${\displaystyle \{\mathbf {F} ^{(e)}\}={\begin{Bmatrix}F_{Hd}\\F_{Vd}\\M_{d}\\F_{Hi}\\F_{Vi}\\M_{i}\end{Bmatrix}}={\frac {1}{\sqrt {2}}}{\begin{Bmatrix}{\frac {2}{3}}P\\-{\frac {20}{27}}P\\-{\frac {4}{27}}PL\\{\frac {1}{3}}P\\-{\frac {7}{27}}P\\+{\frac {2}{27}}PL\end{Bmatrix}}+{\begin{Bmatrix}0\\+{\frac {13}{54}}qL\\-{\frac {1}{108}}qL^{2}\\0\\-{\frac {31}{54}}qL\\{\frac {1}{324}}qL^{2}\end{Bmatrix}}}$

Cálculo de deslocamentos

Uma vez encontrada a matriz de rigidez global e o vetor de forças nodais global se constrói um sistema de equações. Este sistema tem a propriedade de poder ser decomposto em dois subsistemas de equações:

1. El primero de estos sistemas relaciona únicamente los desplazamientos incógnita con algunas de las componentes del vector de fuerzas nodales global y constituye siempre un sistema compatible determinado
2. El segundo subsistema contiene también las reacciones incógnita y una vez resuelto el primer subsistema es de resolución trivial.

Resolvendo o primeiro subsistema são determinados os deslocamentos incógnitos de todos os nós da estrutura. Inserindo a solução do primeiro subsistema no segundo resultam as reaçõnes.

Ilustrando o cálculo dos deslocamentos com um exemplo, consideremos a flexão no plano XY da viga reta da seção anterior, considerando uma viga biarticulada unida em suas extremidade a duas rótulas fxas

${\displaystyle {\frac {1}{\sqrt {2}}}{\begin{Bmatrix}{\frac {2}{3}}P\\-{\frac {20}{27}}P\\-{\frac {4}{27}}PL\\{\frac {1}{3}}P\\-{\frac {7}{27}}P\\+{\frac {2}{27}}PL\end{Bmatrix}}+{\begin{Bmatrix}0\\+{\frac {13}{54}}qL\\-{\frac {1}{108}}qL^{2}\\0\\-{\frac {31}{54}}qL\\{\frac {1}{324}}qL^{2}\end{Bmatrix}}+{\begin{Bmatrix}R_{H1}\\R_{V1}\\0\\R_{H2}\\R_{V2}\\0\end{Bmatrix}}={\begin{bmatrix}{\frac {EA}{L}}&0&0&-{\frac {EA}{L}}&0&0\\0&{\frac {12EI}{L^{3}}}&{\frac {6EI}{L^{2}}}&0&-{\frac {12EI}{L^{3}}}&{\frac {6EI}{L^{2}}}\\0&{\frac {6EI}{L^{2}}}&{\frac {4EI}{L}}&0&-{\frac {6EI}{L^{2}}}&{\frac {2EI}{L}}\\-{\frac {EA}{L}}&0&0&{\frac {EA}{L}}&0&0\\0&-{\frac {12EI}{L^{3}}}&-{\frac {6EI}{L^{2}}}&0&{\frac {12EI}{L^{3}}}&-{\frac {6EI}{L^{2}}}\\0&{\frac {6EI}{L^{2}}}&{\frac {2EI}{L}}&0&-{\frac {6EI}{L^{2}}}&{\frac {4EI}{L}}\end{bmatrix}}{\begin{Bmatrix}0\\0\\\theta _{1}\\0\\0\\\theta _{2}\end{Bmatrix}}}$

As linhas 3 e 6 contém as rotações incógnitas dos extremos da viga e consideradas em conjunto conformam o primeriro subsistema para os deslocamentos. Ignorando os termos nulos e reescrevendo em forma matricial o subsistema de equações para os deslocamentos resulta

${\displaystyle {\frac {1}{\sqrt {2}}}{\begin{Bmatrix}-{\frac {4}{27}}PL\\+{\frac {2}{27}}PL\end{Bmatrix}}+{\begin{Bmatrix}-{\frac {1}{108}}qL^{2}\\{\frac {1}{324}}qL^{2}\end{Bmatrix}}={\begin{bmatrix}{\frac {4EI}{L}}&{\frac {2EI}{L}}\\{\frac {2EI}{L}}&{\frac {4EI}{L}}\end{bmatrix}}{\begin{Bmatrix}\theta _{1}\\\theta _{2}\end{Bmatrix}}~~.}$

A solução fornece os valores das rotações nas extremidades direita e esquerda sob estas cargas

${\displaystyle \theta _{1}=-{\frac {5}{81{\sqrt {2}}}}{\frac {PL^{2}}{EI}}-{\frac {7}{1944}}{\frac {qL^{3}}{EI}}\qquad \theta _{2}=+{\frac {4}{81{\sqrt {2}}}}{\frac {PL^{2}}{EI}}+{\frac {5}{1944}}{\frac {qL^{3}}{EI}}}$

Conhecidos estes valores e inseridos na matriz nas linhas 1, 2, 4 e 5, são determinados os valores das quatro reações hiperestáticas previamente desconhecidas.

Cálculo das reações

Calculados os deslocamentos e resolvendo um sistema de equações, o cálculo das reações é simples e direto

${\displaystyle {\begin{Bmatrix}R_{1}\\R_{2}\\...\\R_{n}\end{Bmatrix}}_{G}={\begin{bmatrix}k_{11}&k_{12}&\cdots &k_{1n}\\k_{21}&k_{22}&\cdots &k_{2n}\\\cdots &\cdots &\ddots &\cdots \\k_{n1}&k_{n2}&\cdots &k_{nn}\end{bmatrix}}_{G}{\begin{Bmatrix}\delta _{1}\\\delta _{2}\\...\\\delta _{n}\end{Bmatrix}}_{G}-{\begin{Bmatrix}F_{1}\\F_{2}\\...\\F_{n}\end{Bmatrix}}_{G}}$

Considerando o mesmo exemplo da última seção o cálculo de reações na viga biarticulada com carga P e q resulta

${\displaystyle {\begin{Bmatrix}R_{H1}\\R_{V1}\\0\\R_{H2}\\R_{V2}\\0\end{Bmatrix}}={\begin{bmatrix}{\frac {EA}{L}}&0&0&-{\frac {EA}{L}}&0&0\\0&{\frac {12EI}{L^{3}}}&{\frac {6EI}{L^{2}}}&0&-{\frac {12EI}{L^{3}}}&{\frac {6EI}{L^{2}}}\\0&{\frac {6EI}{L^{2}}}&{\frac {4EI}{L}}&0&-{\frac {6EI}{L^{2}}}&{\frac {2EI}{L}}\\-{\frac {EA}{L}}&0&0&{\frac {EA}{L}}&0&0\\0&-{\frac {12EI}{L^{3}}}&-{\frac {6EI}{L^{2}}}&0&{\frac {12EI}{L^{3}}}&-{\frac {6EI}{L^{2}}}\\0&{\frac {6EI}{L^{2}}}&{\frac {2EI}{L}}&0&-{\frac {6EI}{L^{2}}}&{\frac {4EI}{L}}\end{bmatrix}}{\begin{Bmatrix}0\\0\\\theta _{1}\\0\\0\\\theta _{2}\end{Bmatrix}}-{\frac {1}{\sqrt {2}}}{\begin{Bmatrix}{\frac {2}{3}}P\\-{\frac {20}{27}}P\\-{\frac {4}{27}}PL\\{\frac {1}{3}}P\\-{\frac {7}{27}}P\\+{\frac {2}{27}}PL\end{Bmatrix}}-{\begin{Bmatrix}0\\+{\frac {13}{54}}qL\\-{\frac {1}{108}}qL^{2}\\0\\-{\frac {31}{54}}qL\\{\frac {1}{324}}qL^{2}\end{Bmatrix}}}$

Introduzindo os valores das rotações nas extremidades e multiplicando a matriz de rigidez pelo vetor de deslocamentos resulta

${\displaystyle {\begin{Bmatrix}R_{H1}\\R_{V1}\\R_{H2}\\R_{V2}\end{Bmatrix}}={\begin{Bmatrix}0\\{\frac {6EI}{L^{2}}}(\theta _{1}+\theta _{2})\\0\\{\frac {6EI}{L^{2}}}(\theta _{1}+\theta _{2})\end{Bmatrix}}-{\frac {1}{\sqrt {2}}}{\begin{Bmatrix}{\frac {2}{3}}P\\-{\frac {20}{27}}P\\{\frac {1}{3}}P\\-{\frac {7}{27}}P\end{Bmatrix}}-{\begin{Bmatrix}0\\+{\frac {13}{54}}qL\\0\\-{\frac {31}{54}}qL\end{Bmatrix}}={\begin{Bmatrix}{\frac {2}{3{\sqrt {2}}}}P\\{\frac {18}{27{\sqrt {2}}}}P-{\frac {20}{81}}qL\\{\frac {1}{3{\sqrt {2}}}}P\\{\frac {5}{27{\sqrt {2}}}}P+{\frac {46}{81}}qL\end{Bmatrix}}}$

Isto completa o cálculo das reações.

Cálculo de esforços

O cálculo de esforços se realiza examinando em coordenadas locais das barras o esforço normal, esforço cortante, momento fletor e momento torsor em cada uma das barras, conhecidos os deslocamentos de todos os nós da estrutura. Isto pode ser realizado usando as matrizes de rigidez expressas em coordenadas locais e os deslocamentos nodais expressos também em coordenadas locais.

Referências

Bibliografia

• Timoshenko, Stephen; Godier J.N. (1951). McGraw-Hill, ed. Theory of elasticity. [S.l.: s.n.] 
• Ortiz Berrocal, Luis (1991). McGraw-Hill, ed. Resistencia de Materiales. [S.l.: s.n.] ISBN 84-7651-512-3 Verifique |isbn= (ajuda) 
• Monleón Cremades, S., Análisis de vigas, arcos, placas y láminas, Ed. UPV, 1999, ISBN 84-7721-769-6.
• Felippa, Carlos A. (2001), «A historical outline of matrix structural analysis: a play in three acts» (PDF), Computers & Structures, ISSN 0045-7949, 79 (14): 1313–1324, doi:10.1016/S0045-7949(01)00025-6, consultado em 5 de junho de 2017, cópia arquivada (PDF) em |arquivourl= requer |arquivodata= (ajuda) 🔗  |título= e |titulo= redundantes (ajuda)
• Felippa, Carlos A. Introduction to Finite Element Method. Fall 2001. University of Colorado. 18 Sept. 2005 <http://www.devdept.com/fem/books.php>
• Robinson, John. Structural Matrix Analysis for the Engineer. New York: John Wiley & Sons, 1966
• Rubinstein, Moshe F. Matrix Computer Analysis of Structures. New Jersey: Prentice-Hall, 1966
• McGuire, W., Gallagher, R. H., and Ziemian, R. D. Matrix Structural Analysis, 2nd Ed. New York: John Wiley & Sons, 2000.

Programas