Coeficientes uj ao longo do domínio do problema

O Método dos Elementos Finitos: uma introdução

O ano era 1901, iniciava-se o século XX e com ele as bases para o que viria a ser conhecido posteriormente como o Método dos Elementos Finitos. De um lado, a indústria aeroespacial em busca de resolver problemas estruturais complexos. Do outro, mas talvez em sintonia, acadêmicos de engenharia estrutural debruçados sobre a possibilidade de expandir os horizontes da análise matricial de estruturas. 

Assim, em torno da década de 50 o Método dos Elementos Finitos começou a ser utilizado na indústria, e em 1960 R.W. Clough de Berkley usou pela primeira vez o termo Finite Element em uma publicação. Não deixando o mérito apenas para Clough, a literatura, e particularmente os alemães, afirma que o professor J.H. Argyris (alemão) atuou ativamente nesse período no desenvolvimento do método.

Como se não fosse o bastante, matemáticos também buscam seu lugar ao sol nessa criação. Há evidências de que na década de 40 o grande matemático Courant já teria publicado um artigo onde ele resolveu um problema de vibrações utilizando elementos Finitos triangulares (sem nomeação do método). E não para por aí. Se formos mais fundo veremos os nomes Jhon William Strutt, Walther Ritz, Boris Galerkin, Hrennikoff, entre outros que contribuíram de alguma forma para o método desde o início do século XX.

Após essa introdução histórica com alta densidade de protagonistas, vamos caminhar cuidadosamente em direção ao entendimento do método sob a luz da matemática e da engenharia.

O que é o Método dos Elementos Finitos e qual a sua utilidade

À princípio, vou ser evasivo no que diz respeito à definição do Método dos Elementos Finitos e de suas aplicações. Dada a grandeza e generalidade do método, não vou me comprometer a limitá-lo a uma gama de problemas de engenharia atuais. As pesquisas podem revelar futuramente que ele é mais poderoso do que já temos considerado. Dessa forma, recorrendo à matemática, me sinto confortável apenas em defini-lo como um método aproximado de solução de equações diferenciais parciais através da minimização de funcionais.

Desculpe-me pela pretensão expressa no título deste tópico, mas vou ficar por aqui com a definição. Acredito fortemente que seja mais didático que você entenda por si só do que se trata o método e qual a sua utilidade através do desenvolvimento que será apresentado abaixo. Um pequeno alerta que faço, porém, é que aqui nos concentraremos em apenas apresentar uma introdução ao método. Portanto, ao final da leitura, saiba que o Método dos Elementos Finitos é/pode isso e muito mais.

Primeiros passos no Método dos elementos Finitos

Vamos iniciar com algo simples. Proponhamos solucionar uma equação diferencial unidirecional. Esse problema terá, além da equação, um domínio no qual ela é válida e valores para a função ou sua derivada no contorno do domínio, caracterizando, dessa forma, um problema clássico de valor de contorno unidimensional. Para aplicação do MEF (Método dos Elementos Finitos) ou FEM em inglês (Finite Element Method), não é necessário que o problema tenha solução analítica. Porém, para efeito de verificação dos resultados obtidos através do método, escolheremos um caso que possua solução analítica. Dito isso, abaixo segue a expressão que, em linguagem matemática, resume todo o texto desenvolvido até aqui neste tópico.

(1)   \begin{gather*} \frac{d^{2}u}{dx^{2}}+f(x)=0 \hspace{10}   em\hspace{2} \Gamma \\ \\ u(x_{i})= \={u}\hspace{5}(condição\hspace{2}de\hspace{2}contorno\hspace{2}essencial) \\ \\ \frac{du}{dx}(x_{f})=\alpha \hspace{5}(condição\hspace{2}de\hspace{2}contorno\hspace{2}natural) \\ \\ \Gamma = [x_{i}, x_{f}] \end{gather*}

Obs.: A ideia por trás dos termos essencial e natural das condições de contorno será revelada no decorrer da leitura.

O Método dos Elementos Finitos nos permitirá obter uma função û(x) que aproxima a função u(x).  Abaixo serão apresentadas as etapas necessárias para isso.

Obtenção da formulação fraca do problema

A equação original (1) representa a formulação forte do problema. Essa formulação é caracterizada por equações de governo, normalmente diferenciais, e condições de contorno que exigem uma solução contínua (suave).  Em contraste, a formulação fraca, com o objetivo de simplificar o problema, faz concessões no que diz respeito à condição de continuidade da solução, possuindo exigências de continuidade mais fracas. Em suma, a formulação fraca de um problema é a forma integral equivalente de sua formulação forte. Obter a formulação fraca da equação (1) será o nosso desafio dessa etapa.

Utilizaremos o método dos resíduos ponderados para obter a forma fraca da nossa equação diferencial. Há outras possibilidades de se obter a forma fraca de equações diferenciais, porém, esse método é o mais genérico entre os seus pares.

Sem delongas, então, integremos a nossa equação “forte” ao longo de todo o seu domínio.

(2)   \begin{equation*} \int_{\Gamma}\left(\frac{d^{2}u}{dx^{2}}+f(x)\right)dx=0 \end{equation*}

Bom, temos agora o problema inicial na forma de uma integral. Mas é equivalente? Vamos analisar: a equação 1 exige que todo o termo dentro da integral acima valha zero ao longo de todo o domínio. Ao integrarmos esse termo e forçarmos o resultado zero não necessariamente garantimos que ele valerá zero ao longo de todo o domínio. Basta pensar numa equação hipotética h(x) com o seguinte gráfico:

Gráfico de função no intervalo [0, 1] com equilíbrio entre área positiva (acima do gráfico) e negativa (abaixo do gráfico)

Pela antissimetria do gráfico no intervalo [x_{i}, x_{f}] fica óbvio que o fato de ser verdade que

(3)   \begin{equation*} \int_{x_{i}}^{x_{f}}h(x)dx=0 \end{equation*}

Não implica em h(x) = 0 em [x_{i}, x_{f}].

Então, para não abrirmos espaço em nossa formulação fraca para uma solução que apresente valores diferentes da formulação forte ao longo do domínio, temos que garantir que o integrante da equação 2 seja nulo ao longo de todo o domínio, assim como na equação 1. Fazer isso não é tão difícil. Basta multiplicar o integrando por uma função qualquer. Sendo uma função qualquer, para que a equação seja satisfeita, obrigatoriamente a outra parte do produto deve ser zero. Pronto: a condição da equação 1 foi preservada. Matematicamente, temos:

(4)   \begin{equation*} \int_{\Gamma}\left(\frac{d^{2}u}{dx^{2}}+f(x)\right)w(x)dx=0\hspace{10}\forallw \end{equation*}

Desvendando parcialmente o nome do método que estamos aplicando (resíduos ponderados), chamaremos w(x) de função de ponderação. Vamos, por conveniência, considerar essa função de ponderação pertencente à um espaço vetorial equivalente ao espaço vetorial da solução do problema. O efeito prático que usufruiremos dessa consideração é que as condições de contorno de w(x) serão as mesmas de u(x).  Aqui vem uma importante observação para que o entendimento não se perca. A condição de contorno essencial só será introduzida na resolução do sistema de equações final, portanto, independente do seu valor, devemos considerá-la igual a zero (para este problema: u(x_{i})=0). Dito isto, vamos manipular a equação 4.

Primeiro, separando a soma em duas integrais temos:

(5)   \begin{equation*} \int_{\Gamma}\frac{d^{2}u}{dx^{2}}w(x)dx+\int_{\Gamma}f(x)w(x)dx=0 \end{equation*}

Segundo, aplicando a integração por partes na primeira integral, temos:

(6)   \begin{equation*} \frac{du}{dx}w(x)\bigg|_{\Gamma}+\int_{\Gamma}\frac{du}{dx}\frac{dw}{dx}dx+\int_{\Gamma}f(x)w(x)dx=0 \end{equation*}

Sabendo que \Gamma=[x_{i}, x_{f}], podemos abrir o termo da esquerda, o que resulta em:

(7)   \begin{equation*} \frac{du}{dx}(x_{f})w(x_{f})-\frac{du}{dx}(x_{i})w(x_{i})-\int_{\Gamma}\frac{du}{dx}\frac{dw}{dx}dx+\int_{\Gamma}f(x)w(x)dx=0 \end{equation*}

Repare que aparece aqui a condição de contorno \frac{du}{dx}(x_{f}) da formulação forte do problema. O fato dela surgir naturalmente no desenvolvimento da forma fraca do problema concede a ela o nome que já conhecemos lá no início – condição de contorno natural. Podemos substituí-la na equação 7 por \alpha, bem como podemos substituir w(x_{i}) por u(x_{i}), que (atenção!) será zero, pois o valor de fato (\={u}) será considerado na equação final conforme já discutimos. Com isso, e reorganizando a equação, apresento a forma fraca do nosso problema.

(8)   \begin{equation*} \int_{\Gamma}\frac{du}{dx}\frac{dw}{dx}dx=\int_{\Gamma}f(x)w(x)dx+\alpha w(x_{f}) \hspace{10}\forallw \end{equation*}

Aproximação por elementos finitos

Vamos agora subdividir o domínio do problema em n-1 elementos finitos e, consequentemente, n pontos nodais.

Domínio do problema discretizado em n elementos

Então, a solução do problema por elementos finitos consiste em se obter coeficiente constantes (u_{j}) para cada um dos n nós resultantes da discretização do domínio e, posteriormente, interpolá-los ao longo de todo o domínio utilizando funções de interpolação (N_{j}(x)).

Coeficientes uj ao longo do domínio do problema

Matematicamente, a aproximação final para u(x) fica da seguinte forma:

(9)   \begin{equation*} û(x) = \sum_{j=1}^{n}N_{j}(x)u_{j} \end{equation*}

A função de interpolação deverá ser escolhida, e os coeficientes constantes serão obtidos como solução do problema em sua forma fraca (8). Veremos a parte de obtenção dos coeficientes mais adiante. Por agora, nos concentremos em entender as tão importantes funções de interpolação.

Definição das funções de interpolação

Expandindo a equação 9, temos a seguinte expressão:

(10)   \begin{equation*} û(x)=N_{1}u_{1}+N_{2}u_{2}+N_{3}u_{3}+...+N_{n-2}u_{n-2}+N_{n-1}u_{n-1}+N_{n}u_{n} \end{equation*}

Essa expressão “aberta” nos ajudará na análise que se seguirá. Observando o gráfico da figura acima vemos que no nó a função de aproximação û(x) é igual à constante u associada àquele nó.

(11)   \begin{gather*} û(x_{1})=u_{1} \\ \\ û(x_{2})=u_{2} \\ \\ û(x_{3})=u_{3} \\ \\ \vdots \\ \\ û(x_{n-2})=u_{n-2} \\ \\ û(x_{n-1})=u_{n-1} \\ \\ û(x_{n})=u_{n} \end{gather*}

Note que para isso ocorrer é necessário que

(12)   \begin{gather*} N_{1} = 1, \hspace{2} N_{2}=0, \hspace{2} N_{3}=0, \hspace{2}\hdots, \hspace{2} N_{n-2}=0, \hspace{2} N_{n-1}=0 \hspace{5} e \hspace{5} N_{n}=0 \hspace{10} para \hspace{2}x=x_{1} \\ \\ N_{1} = 0, \hspace{2} N_{2}=1, \hspace{2} N_{3}=0, \hspace{2}\hdots, \hspace{2} N_{n-2}=0, \hspace{2} N_{n-1}=0 \hspace{5} e \hspace{5} N_{n}=0 \hspace{10} para \hspace{2}x=x_{2} \\ \\ \vdots \\ \\ N_{1} = 0, \hspace{2} N_{2}=0, \hspace{2} N_{3}=0, \hspace{2}\hdots, \hspace{2} N_{n-2}=0, \hspace{2} N_{n-1}=1 \hspace{5} e \hspace{5} N_{n}=0 \hspace{10} para \hspace{2}x=x_{n-1} \\ \\ N_{1} = 0, \hspace{2} N_{2}=0, \hspace{2} N_{3}=0, \hspace{2}\hdots, \hspace{2} N_{n-2}=0, \hspace{2} N_{n-1}=0 \hspace{5} e \hspace{5} N_{n}=1 \hspace{10} para \hspace{2}x=x_{n} \\ \\ \end{gather*}

Em palavras significa dizer que a função de interpolação do nó j (N_{j}) valerá 1 no nó j e zero nos demais nós. Simples, não? Vejamos um exemplo: a função de interpolação linear para o nó 2.

Exemplo de uma função de interpolação linear para o nó 2

Em diversas aplicações do Método dos Elementos Finitos são utilizadas funções de interpolação lineares. Portanto, nos limitaremos a elas neste artigo. As equações que reproduzem o que está exemplificado na imagem acima estão apresentadas abaixo.

(13)   \begin{equation*} N_{j}(x)= \begin{cases} 1, \hspace{20} para \hspace{5} x=x_{j}\\ 0, \hspace{20} para \hspace{5} x=x_{i}(i\neq j) \end{cases}\ \end{equation*}

(14)   \begin{equation*} N_{j}(x)= \begin{cases} -\frac{x-x_{j}}{(x_{j+1}-x_{j}}+1, \hspace{20}x_{j}\leq x \leq x_{j+1}\\ \\ \frac{x-x_{j-1}}{(x_{j}-x_{j-1}}, \hspace{47}x_{j-1}\leq x \leq x_{j} \end{cases}\ \end{equation*}

Aplicando a ideia central do Método dos Elementos Finitos

Vamos agora desvendar por completo o método dos resíduos ponderados introduzido anteriormente. Já entendemos a ponderação, basta agora compreender a ideia de resíduo. Temos discorrido bastante sobre a função û(x) que se propõe a aproximar a solução do nosso problema. Vamos agora usá-la na formulação forte do problema (Equação 1). Considerando que ela não é exatamente a função u(x) (digo considerando, pois em alguns casos a solução aproximada pode ser igual  à solução exata), no lado direito da equação surgirá um valor diferente de zero. Esse valor, representado por R, é chamado de resíduo. Temos a seguinte equação então:

(15)   \begin{equation*} \frac{d^{2}û}{dx^{2}}+f(x)=R(x) \hspace{10}   em\hspace{2} \Gamma \\ \\ \end{equation*}

Vamos agora exigir que esse erro seja zero. Mas não na formulação forte (senão caímos no mesmo problema inicial), e sim na formulação fraca. Então, relembrando o que desenvolvemos no tópico onde discorremos sobre a obtenção da formulação fraca, teríamos:

(16)   \begin{equation*} \int_{\Gamma}Rwdx=0 \end{equation*}

Que pode ser escrita em termos da função de aproximação simplesmente trocando R(x) pelo lado esquerdo da equação 15.

(17)   \begin{equation*} \int_{\Gamma}\left(\frac{d^{2}û}{dx^{2}}+f(x)\right)wdx=0 \end{equation*}

E, por analogia com o que já foi visto, sabemos que se desenvolvermos essa equação chagaremos na forma fraca apresentada abaixo.

(18)   \begin{equation*} \int_{\Gamma}\frac{dû}{dx}\frac{dw}{dx}dx=\int_{\Gamma}f(x)w(x)dx+\alpha w(x_{f}) \hspace{10}\forallw \end{equation*}

Já definimos também que û(x) pode ser aproximado por:

(19)   \begin{equation*} û(x) = \sum_{j=1}^{n}N_{j}(x)u_{j} \end{equation*}

Falta-nos agora definir a função qualquer w(x) para prosseguirmos com o desenvolvimento do problema. Seguindo o mesmo princípio utilizado para a função objetivo (u(x)), vamos utilizar uma aproximação para a função peso. Além disso, nos apropriando de um método largamente conhecido criado por Galerkin, e nomeado com seu próprio nome (Método de Galerkin), vamos fazer a aproximação da função peso com as mesmas funções de interpolação utilizadas para aproximar a função u(x). Dessa forma, temos:

(20)   \begin{equation*} \^w(x) = \sum_{i=1}^{n}N_{i}(x)w_{i} \end{equation*}

Enfim, podemos aplicar essas aproximações para û(x) e \^w(x) na formulação fraca e prosseguir o desenvolvimento da solução aproximada com elementos finitos.

Formulação discreta do problema

Com o esforço aplicado até aqui, temos a nossa proposta de solução do problema no seguinte estágio:

(21)   \begin{equation*} \int_{\Gamma}\frac{dû}{dx}\frac{d\^w}{dx}dx=\int_{\Gamma}f(x)\^w(x)dx+\alpha \^w(x_{f}) \hspace{10}\forallw \end{equation*}

Onde

(22)   \begin{gather*} û(x) = \sum_{j=1}^{n}N_{j}(x)u_{j} \\\\ \^w(x) = \sum_{i=1}^{n}N_{i}(x)w_{i} \\\\ n = Número\hspace{2} de\hspace{2} nós\hspace{2} resultantes\hspace{2} da\hspace{2} discretização\hspace{2} do\hspace{2} domínio\hspace{2} do\hspace{2} problema. \end{gather*}

Injetando as equações das aproximações dentro da formulação fraca (21), vem:

(23)   \begin{equation*} \int_{\Gamma}\frac{d}{dx}\left(\sum_{j=1}^{n}N_{j}u_{j}\right)\frac{d}{dx}\left(\sum_{i=1}^{n}N_{i}w_{i}\right)dx=\int_{\Gamma}f(x)\sum_{i=1}^{n}N_{i}w_{i}dx+\alpha \sum_{i=1}^{n}N_{i}(x_{f})w_{i} \end{equation*}

Focando primeiro no lado esquerdo da equação acima, aplicaremos a seguinte sequência de manipulação:

  • Lançando mão da propriedade da derivada que diz que a derivada da soma é igual a soma das derivadas, o operador de derivada pode ser transportado para dentro do somatório.

(24)   \begin{equation*} \int_{\Gamma}\sum_{j=1}^{n}\frac{d}{dx}\left(N_{j}u_{j}\right)\sum_{i=1}^{n}\frac{d}{dx}\left(N_{i}w_{i}\right)dx=\hdots \end{equation*}

  • Recordando que os coeficientes u_{j} e w_{i} são constantes, eles podem sair da derivada.

(25)   \begin{equation*} \int_{\Gamma}\sum_{j=1}^{n}u_{j}\frac{dN_{j}}{dx}\sum_{i=1}^{n}w_{i}\frac{dN_{i}}{dx}dx=\hdots \end{equation*}

  • Pelo mesmo motivo anterior, os coeficientes podem ainda ser enviados para fora da integral. Nesse caso, os somatórios devem ir juntos. E isso pode ser realizado sem problemas, pois há a propriedade das integrais que garante que a integral da soma é igual a soma das integrais.

(26)   \begin{equation*} \sum_{i=1}^{n}\sum_{j=1}^{n}\left(\int_{\Gamma}\frac{dN_{j}}{dx}\frac{dN_{i}}{dx}dx\right)w_{i}u_{j}=\hdots \end{equation*}

Vamos agora ao lado direito da equação. Além de ser mais simples, com a experiência obtida na manipulação do termo do lado esquerdo da equação, sinto adequado apresentar diretamente o resultado final, evitando ser prolixo.

(27)   \begin{equation*} \dots=\sum_{i=1}^{n}\left(\int_{\Gamma}f(x)N_{i}dx+\alpha N_{i}(x_{f}) \right) w_{i} \end{equation*}

No geral, então, temos até o momento:

(28)   \begin{equation*} \sum_{i=1}^{n}\sum_{j=1}^{n}\left(\int_{\Gamma}\frac{dN_{j}}{dx}\frac{dN_{i}}{dx}dx\right)w_{i}u_{j}=\sum_{i=1}^{n}\left(\int_{\Gamma}f(x)N_{i}dx+\alpha N_{i}(x_{f})\right)w_{i} \end{equation*}

Sinto necessidade nesse momento de recordá-lo de que estamos buscando pelos valores dos coeficientes u_{j} (que são as aproximações para u(x) nos nós) para, então, obtermos a aproximação por elementos finitos û(x) fazendo \sum_{j=1}^{n}N_{j}u_{j}. Devo lembrá-lo também de que N_{j}(x) e N_{i}(x) são as funções de interpolação que nós devemos escolher, logo não são incógnitas. f(x) e g são dados do problema e w_{i} são coeficientes de ponderação arbitrários. Em resumo, as incógnitas do problema acima são os n coeficientes u_{j}. Com isso esclarecido, simplificaremos a representação dos termos que não são incógnitas que se encontram dentro dos parênteses substituindo-os pelas letras k e f.

(29)   \begin{equation*} \sum_{i=1}^{n}\sum_{j=1}^{n}k_{ij}w_{i}u_{j}=\sum_{i=1}^{n}f_{i}w_{i} \end{equation*}

Bem melhor, não acha? Repare que o k recebe tanto o índice i quanto o índice j dado que o termo que ele representa depende dos dois índices. Vamos agora “abrir” esses somatórios. Pelo motivo que será desvendado mais adiante, farei isso em um formato especial como segue.

(30)   \begin{gather*} i=1 \hspace{20} k_{11}w_{1}u_{1}+k_{12}w_{1}u_{2}+\hdots+k_{1n}w_{1}u_{n} \hspace{35} f_{1}w_{1} \\\\ i=2 \hspace{20} k_{21}w_{2}u_{1}+k_{22}w_{2}u_{2}+\hdots+k_{2n}w_{2}u_{n} \hspace{10}=\hspace{10} f_{2}w_{2} \\\\ +\\ \vdots \\ +\\\\ i=1 \hspace{20} k_{n1}w_{n}u_{1}+k_{n2}w_{n}u_{2}+\hdots+k_{nn}w_{n}u_{n} \hspace{35} f_{n}w_{n} \end{gather*}

Subtraindo os termos que se encontram no lado direito em toda a equação e colocando os coeficientes w_{i} em evidência, temos:

(31)   \begin{gather*} i=1 \hspace{20} w_{1}(k_{11}u_{1}+k_{12}u_{2}+\hdots+k_{1n}u_{n} -f_{1})\hspace{33}  \\\\ i=2 \hspace{20} w_{2}(k_{21}u_{1}+k_{22}u_{2}+\hdots+k_{2n}u_{n}-f_{2}) \hspace{10}=\hspace{10} 0 \\\\ +\\ \vdots \\ +\\\\ i=1 \hspace{20} w_{n}(k_{n1}u_{1}+k_{n2}u_{2}+\hdots+k_{nn}u_{n}-f_{n}) \hspace{33} \end{gather*}

A solução trivial onde todos os coeficientes são nulos não nos interessa, pois permitiria qualquer solução para o problema. Então, necessariamente os termos entre os parênteses devem ser nulos. Com isso, temos o seguinte sistema de n equações.

(32)   \begin{equation*} \left\{ \begin{array}{ll} k_{11}u_{1}+k_{12}u_{2}+\hdots+k_{1n}u_{n} -f_{1}=0\\ k_{21}u_{1}+k_{22}u_{2}+\hdots+k_{2n}u_{n}-f_{2}=0\\ \hspace{50}\vdots\\ k_{n1}u_{1}+k_{n2}u_{2}+\hdots+k_{nn}u_{n}-f_{n}=0 \end{array}\right \end{equation*}

Em resumo, o sistema final de equações do Método de Elementos Finitos para o problema será:

(33)   \begin{equation*} \sum_{j=1}^{n}k_{ij}u_{j}=f_{i} \hspace{20} com \hspace{4}i=1,...,n \end{equation*}

Onde

(34)   \begin{gather*} k_{ij}=\int_{\Gamma}\frac{dN_{j}}{dx}\frac{dN_{i}}{dx}dx \\\\ f_{i}=\int_{\Gamma}f(x)N_{i}dx-\alpha N_{i}(x_{f}) \end{gather*}

Temos n equações para n incógnitas. Mas essas equações são linearmente independentes? A resposta é: nem todas. No caso desse problema duas das equações são LD entre si. Mas não há motivo para preocupação, pois na verdade não temos n incógnitas e sim n-1. O valor do coeficiente u_{j} em x_{i} é igual a \={u}. Essa é a condição de contorno essencial do nosso problema (aliás, por isso que ela recebe esse nome). Com o número de incógnitas reduzido em 1, é aceitável que duas das equações sejam LD entre si sem prejuízo à resolução do sistema.

Vamos agora pôr em prática toda essa teoria aprendida.

Resolução de uma equação diferencial utilizando o Método dos Elementos Finitos

Vamos resolver a seguinte equação diferencial:

(35)   \begin{gather*} \frac{d^{2}u}{dx^{2}}+f(x)=0 \hspace{10}   em\hspace{2} \Gamma \\ \\ u(0)= 1\\ \\ \frac{du}{dx}(1)=\alpha = 0 \\ \\ \Gamma = [0, 1] \\ \\ f(x) = 2 \end{gather*}

Obtenção da forma fraca do problema

Como essa equação é similar à equação geral que foi utilizada nos tópicos anteriores, podemos obter a forma fraca do problema simplesmente substituindo f(x), \alpha e \Gamma da equação 8 pelos valores deste problema.

(36)   \begin{equation*} \int_{0}^{1}\frac{du}{dx}\frac{dw}{dx}dx=\int_{0}^{1}2w(x)dx \hspace{10}\forallw \end{equation*}

Discretização do domínio em elementos finitos

O domínio, nesse caso linear, será dividido em 3 elementos finitos. Isso resultará em 4 nós, enumerados conforme a figura a seguir.

Discretização do domínio

 

Definição das funções de interpolação

Para cada um dos nós teremos uma função de interpolação. Elas serão, por escolha, lineares. Abaixo estão apresentadas as equações de interpolação de cada nó e as respectivas representações gráficas delas.

(37)   \begin{equation*} N_{1}= \left\{ \begin{array}{ll} 1-3x \hspace{40} para \hspace{10} 0\leq x\leq 1/3\\ 0 \hspace{40} para \hspace{10} 1/3<x\leq 1 \end{array}\right \end{equation*}

Função de interpolação

 

(38)   \begin{equation*} N_{2}= \left\{ \begin{array}{ll} 3x \hspace{40} para \hspace{10} 0\leq x\leq 1/3\\ 2-3x \hspace{40} para \hspace{10} 1/3<x\leq 2/3 \\ 0 \hspace{40} para \hspace{10} 2/3<x\leq 1 \end{array}\right \end{equation*}

Função de interpolação 2

 

(39)   \begin{equation*} N_{3}= \left\{ \begin{array}{ll} 0 \hspace{40} para \hspace{10} 0\leq x\leq 1/3\\ 3x-1 \hspace{40} para \hspace{10} 1/3<x\leq 2/3 \\ 3-3x \hspace{40} para \hspace{10} 2/3<x\leq 1 \end{array}\right \end{equation*}

Função de interpolação 3

 

(40)   \begin{equation*} N_{4}= \left\{ \begin{array}{ll} 0 \hspace{40} para \hspace{10} 0\leq x\leq 2/3\\ 3x-2 \hspace{40} para \hspace{10} 2/3<x\leq 1 \end{array}\right \end{equation*}

Função de interpolação 4

 

Repare que seguimos fielmente a regra das funções de interpolação: ter valor 1 no nó ao qual ela diz respeito e valer 0 nos demais nós.

Montagem do sistema algébrico de elementos finitos

Construiremos o sistema algébrico para este problema conforme a Equação 33 já apresentada. Para isso, é necessário calcular os valores de k_{ij} e f_{i}. Como o problema que estamos resolvendo é similar ao problema considerado no desenvolvimento teórico, esses valores podem ser obtidos pelas expressões da Equação 34, substituindo apenas os valores de f(x) e \alpha. Temos, então, a seguinte sequência de cálculos:

(41)   \begin{equation*} \frac{dN_{1}}{dx}= \left\{ \begin{array}{ll} -3 \hspace{40} para \hspace{10} 0\leq x\leq 1/3\\ 0 \hspace{40} para \hspace{10} 1/3<x\leq 1 \end{array}\right \end{equation*}

 

(42)   \begin{equation*} \frac{dN_{2}}{dx}= \left\{ \begin{array}{ll} 3 \hspace{40} para \hspace{10} 0\leq x\leq 1/3\\ -3 \hspace{40} para \hspace{10} 1/3<x\leq 2/3 \\ 0 \hspace{40} para \hspace{10} 2/3<x\leq 1 \end{array}\right \end{equation*}

 

(43)   \begin{equation*} \frac{dN_{3}}{dx}= \left\{ \begin{array}{ll} 0 \hspace{40} para \hspace{10} 0\leq x\leq 1/3\\ 3 \hspace{40} para \hspace{10} 1/3<x\leq 2/3 \\ -3 \hspace{40} para \hspace{10} 2/3<x\leq 1 \end{array}\right \end{equation*}

 

(44)   \begin{equation*} \frac{dN_{4}}{dx}= \left\{ \begin{array}{ll} 0 \hspace{40} para \hspace{10} 0\leq x\leq 2/3\\ 3 \hspace{40} para \hspace{10} 2/3<x\leq 1 \end{array}\right \end{equation*}

 

(45)   \begin{equation*} k_{11}=\int_{0}^{1}\frac{dN_{1}}{dx}\frac{dN_{1}}{dx}dx=\int_{0}^{1/3}(-3)^2dx=3 \end{equation*}

(46)   \begin{equation*} k_{12}=\int_{0}^{1}\frac{dN_{1}}{dx}\frac{dN_{2}}{dx}dx=\int_{0}^{1/3}(-3)\times 3dx=-3 \end{equation*}

(47)   \begin{equation*} k_{13}=\int_{0}^{1}\frac{dN_{1}}{dx}\frac{dN_{3}}{dx}dx=\int_{0}^{1}0dx=0 \end{equation*}

(48)   \begin{equation*} k_{14}=\int_{0}^{1}\frac{dN_{1}}{dx}\frac{dN_{4}}{dx}dx=\int_{0}^{1}0dx=0 \end{equation*}

(49)   \begin{equation*} k_{21}=\int_{0}^{1}\frac{dN_{2}}{dx}\frac{dN_{1}}{dx}dx=\int_{0}^{1/3}3\times (-3)dx=-3 \end{equation*}

(50)   \begin{equation*} k_{22}=\int_{0}^{1}\frac{dN_{2}}{dx}\frac{dN_{2}}{dx}dx=\int_{0}^{1/3}3\times 3dx+\int_{1/3}^{2/3}(-3)\times (-3)dx=6 \end{equation*}

(51)   \begin{equation*} k_{23}=\int_{0}^{1}\frac{dN_{2}}{dx}\frac{dN_{3}}{dx}dx=\int_{1/3}^{2/3}(-3)\times 3dx=-3 \end{equation*}

(52)   \begin{equation*} k_{24}=\int_{0}^{1}\frac{dN_{2}}{dx}\frac{dN_{4}}{dx}dx=\int_{0}^{1}0dx=0 \end{equation*}

(53)   \begin{equation*} k_{31}=\int_{0}^{1}\frac{dN_{3}}{dx}\frac{dN_{1}}{dx}dx=\int_{0}^{1}0dx=0 \end{equation*}

(54)   \begin{equation*} k_{32}=\int_{0}^{1}\frac{dN_{3}}{dx}\frac{dN_{2}}{dx}dx=\int_{1/3}^{2/3}3\times (-3)dx=-3 \end{equation*}

(55)   \begin{equation*} k_{33}=\int_{0}^{1}\frac{dN_{3}}{dx}\frac{dN_{3}}{dx}dx=\int_{1/3}^{2/3}3^2dx+\int_{2/3}^{1}(-3)^2dx=6 \end{equation*}

(56)   \begin{equation*} k_{34}=\int_{0}^{1}\frac{dN_{3}}{dx}\frac{dN_{4}}{dx}dx=\int_{2/3}^{1}(-3)\times 3dx=-3 \end{equation*}

(57)   \begin{equation*} k_{41}=\int_{0}^{1}\frac{dN_{4}}{dx}\frac{dN_{1}}{dx}dx=\int_{0}^{1}0dx=0 \end{equation*}

(58)   \begin{equation*} k_{42}=\int_{0}^{1}\frac{dN_{4}}{dx}\frac{dN_{2}}{dx}dx=\int_{0}^{1}0dx=0 \end{equation*}

(59)   \begin{equation*} k_{43}=\int_{0}^{1}\frac{dN_{4}}{dx}\frac{dN_{3}}{dx}dx=\int_{2/3}^{1}3\times (-3)dx=-3 \end{equation*}

(60)   \begin{equation*} k_{44}=\int_{0}^{1}\frac{dN_{4}}{dx}\frac{dN_{4}}{dx}dx=\int_{2/3}^{1}3^2dx=3 \end{equation*}

 

Nesse momento preciso destacar duas características dos coeficientes k_{ij} que podem ser extraídas dessa sequência árdua de cálculo de integrais presenciada acima.

  1. O valor de k_{ij} sempre será zero se os nós i e j não estiverem conectados por um elemento finito comum; e
  2. Os valores k_{ij} são simétricos, isto é, k_{ij} = k_{ji}.

Saber isso anteriormente nos teria poupado algumas integrações. Mas assim é melhor, pois o leitor pôde visualizar essas duas afirmações sendo validadas.

Vamos agora calcular os valores de f_{i}.

(61)   \begin{equation*} f_{1}=\int_{0}^{1}2\times N_{1}}dx - 0 \times  N_{1}(1)=\int_{0}^{1/3}2 \times (1-3x) dx=1/3 \end{equation*}

(62)   \begin{equation*} f_{2}=\int_{0}^{1}2\times N_{2}}dx - 0 \times  N_{2}(1)=\int_{0}^{1/3}2 \times 3x dx + \int_{1/3}^{2/3}2 \times (2-3x) dx=2/3 \end{equation*}

(63)   \begin{equation*} f_{3}=\int_{0}^{1}2\times N_{3}}dx - 0 \times  N_{3}(1)=\int_{1/3}^{2/3}2 \times (3x-1) dx + \int_{2/3}^{1}2 \times (3-3x) dx=2/3 \end{equation*}

(64)   \begin{equation*} f_{4}=\int_{0}^{1}2\times N_{4}}dx - 0 \times  N_{4}(1)=\int_{2/3}^{1}2 \times (3x-2) dx=1/3 \end{equation*}

 

Assim, já podemos montar o sistema de equações de elementos finitos.

(65)   \begin{equation*} \left\{ \begin{array}{ll} 3u_{1}+-3u_{2}=1/3\\ -3u_{1}+6u_{2}-3u_{3}=2/3\\ -3u_{2}+6u_{3}-3u_{4}=2/3\\ -3u_{3}+3u_{4}=1/3\\ \end{array}\right \end{equation*}

Que, em formato matricial, se apresenta da seguinte forma:

(66)   \begin{equation*} \begin{bmatrix} 3 & -3 & 0 & 0 \\ -3 & 6 & -3 & 0 \\ 0 & -3 & 6 & -3 \\ 0 & 0 & -3 & 3 \end{bmatrix} \begin{bmatrix} u_{1}\\ u_{2} \\ u_{3} \\ u_{4} \end{bmatrix} = \begin{bmatrix} 1/3 \\ 2/3 \\ 2/3 \\ 1/3 \end{bmatrix} \end{equation*}

 

Resolução do sistema de equações de elementos finitos

Se tentarmos resolver o sistema de equações final que obtivemos para o nosso exemplo não teremos êxito, pois o sistema obtido não tem solução. Precisamos utilizar a condição de contorno essencial (u(0)=u_{1}=1) para reduzir o número de incógnitas. O sistema então passa a ser:

(67)   \begin{equation*} \begin{bmatrix} 3 & -3 & 0 & 0 \\ -3 & 6 & -3 & 0 \\ 0 & -3 & 6 & -3 \\ 0 & 0 & -3 & 3 \end{bmatrix} \begin{bmatrix} 1\\ u_{2} \\ u_{3} \\ u_{4} \end{bmatrix} = \begin{bmatrix} 1/3 \\ 2/3 \\ 2/3 \\ 1/3 \end{bmatrix} \end{equation*}

Que pode ser simplificado para:

(68)   \begin{equation*} \begin{bmatrix} 6 & -3 & 0 \\ -3 & 6 & -3 \\ 0 & -3 & 3 \end{bmatrix} \begin{bmatrix} u_{2} \\ u_{3} \\ u_{4} \end{bmatrix} = \begin{bmatrix} 2/3 + 3 \times 1\\ 2/3 + 0 \times 1\\ 1/3 + 0 \times 1 \end{bmatrix}= \begin{bmatrix} 11/3\\ 2/3\\ 1/3 \end{bmatrix} \end{equation*}

Resolvendo-o, temos o seguinte resultado:

(69)   \begin{equation*} u_{2}=14/9; \hspace{20} u_{3}=17/9;\hspace{20} u_{4}=2; \end{equation*}

Pronto! Agora sabemos os valores da função u(x) nos nós resultantes da discretização em elementos finitos que fizemos no domínio do problema. Vamos agora gerar a função de aproximação û(x).

Função de aproximação de elementos finitos

Lembrando o que já foi visto, a função de aproximação se apresenta da seguinte forma:

(70)   \begin{equation*} û(x) = \sum_{j=1}^{n}N_{j}(x)u_{j} \end{equation*}

No problema em questão será, então:

(71)   \begin{equation*} û(x) = 1N_{1}(x) + 14/9N_{2}(x) + 17/9N_{3}(x) + 2N_{4}(x) \end{equation*}

Substituindo as funções de aproximação que foram definidas para o problema, tem-se por fim:

(72)   \begin{equation*} û(x) = \left\{ \begin{array}{ll} 5/3x+1 \hspace{40} para \hspace{10} 0\leq x\leq 1/3\\ x+11/9 \hspace{40} para \hspace{10} 1/3<x\leq 2/3 \\ x/3+5/3 \hspace{40} para \hspace{10} 2/3<x\leq 1 \end{array}\right \end{equation*}

Sabendo que a solução exata para esse problema é a equação

(73)   \begin{equation*} u(x) = -2x^2+2x+1 \end{equation*}

Podemos plotar o gráfico da curva da aproximação de elementos finitos sobre o gráfico da curva exata. Fazendo isso obtemos:

Comparação - solução exata e solução aproximada - elementos finitos
Solução exata em azul e solução aproximada em vermelho

Acho que, observando o gráfico acima, o leitor há de concordar que conseguimos aproximar a solução do problema utilizando o Método dos Elementos Finitos. E se usássemos funções de interpolação quadráticas, como seria o resultado? Como a solução exata desse problema é quadrática, a resposta é que encontraremos exatamente a solução do problema. Com essa observação final, fecho esse primeiro exemplo de aplicação do método de elementos finitos.

Resolução do problema da barra carregada axialmente utilizando o Método dos Elementos Finitos

Vamos agora a um exemplo específico de engenharia: obter os deslocamentos ao longo de uma barra carregada axialmente utilizando elementos finitos. Seguindo a notação que vem sendo utilizada, a incógnita do problema (campo de deslocamento), será chamada de u(x).

Entrando em mais detalhes, o problema consiste em uma barra de comprimento L com geometria de seção constante ao longo do comprimento. Ela está engastada em uma extremidade e livre na outra e possui um carregamento distribuído f(x) ao longo do comprimento e outro pontual R na extremidade livre, ambos na direção axial da barra. O esquema abaixo ilustra o problema.

Barra carregada axialmente

Vamos às etapas do método de elementos finitos. Diferentemente do caso anterior, aqui iniciaremos um passo atrás. Desenvolveremos a formulação forte do problema.

Desenvolvimento da formulação forte do problema

Consideremos um elemento infinitesimal de comprimento \Delta x da barra. Atuará nele a força externa f(x+\Delta x/2) e as forças internas (N(x)) em cada uma de suas extremidades.

Elemento infinitesimal - barra carregada axialmente

Para que o elemento infinitesimal esteja em equilíbrio, é necessário que a soma de todas as forças atuando sobre ele seja igual a zero. Assim,

(74)   \begin{equation*} N(x+\Delta x) - N(x) + f(x+\frac{\Delta x}{2}) \Delta x = 0 \end{equation*}

Dividindo todos os termos por \Delta x, temos:

(75)   \begin{equation*} \frac{N(x+\Delta x) - N(x)}{\Delta x} + f(x+\frac{\Delta x}{2}) = 0 \end{equation*}

Fazendo o limite quando \Delta x tende a zero, teremos o primeiro termo do lado esquerdo como sendo a própria definição de derivada e no segundo termo a parcela \Delta x/2 poderá ser desprezada. Disso resulta:

(76)   \begin{equation*} \frac{dN(x)}{dx} + f(x) = 0 \end{equation*}

Agora vamos definir o N(x), a força interna. Sabe-se pela Lei de Hooke que

(77)   \begin{equation*} N(x) = EA\epsilon \end{equation*}

E, por definição, para o elemento infinitesimal

(78)   \begin{equation*} \epsilon = \frac{u(x+\Delta x)-u(x)}{\Delta x} \end{equation*}

Analogamente ao caso anterior, fazendo o limite quando \Delta x tende a zero, teremos

(79)   \begin{equation*} \epsilon = \frac{du}{dx} \end{equation*}

Substituindo então na equação 77, obtemos:

(80)   \begin{equation*} N(x) = EA\frac{du}{dx} \end{equation*}

Levando esse resultado para a equação 76, obtemos, por fim, a equação “forte” que modela o nosso problema.

(81)   \begin{equation*} \frac{d^2u}{dx^2} + f(x) = 0 \hspace{20} em \hspace{5}[0, L] \end{equation*}

Para finalizar esse tópico ainda é necessário definir as condições de contorno do problema. Como em x=0 a barra está engastada, é fácil concluir que o deslocamento nesse ponto será zero. Além disso, também é possível notar que em x=L a força interna na barra (N(x)) será igual a força R aplicada na extremidade da barra. Assim, podemos intuir diretamente as seguintes condições de contorno:

(82)   \begin{gather*} u(0)=0 \\ \\ EA\frac{du}{dx}(L)=R \end{gather*}

Feito! Vamos agora obter a forma fraca.

Obtenção da forma fraca do problema

Se o leitor é atento, vai perceber que a forma forte do problema da barra é análoga à forma forte do problema anterior. Com isso, por semelhança, já posso escrever direto que a forma fraca desse problema será:

(83)   \begin{equation*} EA\int_{0}^{L}\frac{du}{dx}\frac{dw}{dx}dx=\int_{0}^{L}f(x)wdx+Rw(L) \end{equation*}

Discretização do domínio em elementos finitos

O domínio será dividido em 3 elementos finitos, conforme no problema anterior. Para agilizar o processo vamos considerar o comprimento L da barra como sendo 1m. Dessa forma, a discretização do domínio recairá exatamente na discretização do exemplo anterior. Como vamos considerar funções de interpolação lineares nesse exemplo, poderemos reutilizar as funções de interpolação calculadas no exemplo anterior. Como o domínio ([0,1]) e a discretização (3 elementos finitos de tamanhos iguais) para os dois exemplos são semelhantes, as funções de interpolação lineares também serão iguais. Dessa forma, vamos direto à montagem do sistema de equações de elementos finitos.

Montagem do sistema algébrico de elementos finitos

Por similaridade com o exemplo anterior, as equações de elementos finitos para se definir os valores de k_{ij} e de f_{i} são as apresentadas abaixo.

(84)   \begin{gather*} k_{ij}=AE\int_{\Gamma}\frac{dN_{j}}{dx}\frac{dN_{i}}{dx}dx \\\\ f_{i}=\int_{\Gamma}f(x)N_{i}dx-\alpha N_{i}(x_{f}) \end{gather*}

Para obter os valores de k_{ij} vamos adotar E igual a 2\times 10^{11}Pa e A igual a 0,005 m^2. Com isso, o produto AE será igual a 1\times 10^9 N. Como esse fator AE é o único elemento que difere do exemplo anterior, para obter a matriz K deste problema, basta então aplicar esse fator à matriz do exemplo anterior. Dessa forma, temos que:

(85)   \begin{equation*} K=1\times 10^9 \times \begin{bmatrix} 3 & -3 & 0 & 0 \\ -3 & 6 & -3 & 0 \\ 0 & -3 & 6 & -3 \\ 0 & 0 & -3 & 3 \end{bmatrix} \end{equation*}

Para calcular os valores de f_{i} vamos assumir que f(x) vale 10N/m e que R=50N. Dessa forma, seguem os cálculos:

(86)   \begin{equation*} f_{1}=\int_{0}^{1}10\times N_{1}}dx +50 N_{1}(1)=\int_{0}^{1/3}10 \times (1-3x) dx+50\times 0=5/3 \end{equation*}

(87)   \begin{equation*} f_{2}=\int_{0}^{1}10\times N_{2}}dx + 50 N_{2}(1)=\int_{0}^{1/3}10 \times 3x dx + \int_{1/3}^{2/3}10 \times (2-3x) dx+50\times 0=10/3 \end{equation*}

(88)   \begin{equation*} f_{3}=\int_{0}^{1}10\times N_{3}}dx + 50 N_{3}(1)=\int_{1/3}^{2/3}10 \times (3x-1) dx + \int_{2/3}^{1}10 \times (3-3x) dx+50\times 0=10/3 \end{equation*}

(89)   \begin{equation*} f_{4}=\int_{0}^{1}10\times N_{4}}dx + 50   N_{4}(1)=\int_{2/3}^{1}10 \times (3x-2) dx+50\times 1=155/3 \end{equation*}

Assim, juntando as partes, temos o sistema de equação para o problema de elementos finitos:

(90)   \begin{equation*} 1\times 10^9 \times \begin{bmatrix} 3 & -3 & 0 & 0 \\ -3 & 6 & -3 & 0 \\ 0 & -3 & 6 & -3 \\ 0 & 0 & -3 & 3 \end{bmatrix} \begin{bmatrix} u_{1}\\ u_{2} \\ u_{3} \\ u_{4} \end{bmatrix} = \begin{bmatrix} 5/3 \\ 10/3 \\ 10/3 \\ 155/3 \end{bmatrix} \end{equation*}

Resolução do sistema de equações de elementos finitos

O primeiro passo aqui consiste em reduzir o sistema de equações utilizando a condição de contorno essencial – u(0)=0. Fazendo isso, temos:

(91)   \begin{equation*} 1\times 10^9 \times \begin{bmatrix} 3 & -3 & 0 & 0 \\ -3 & 6 & -3 & 0 \\ 0 & -3 & 6 & -3 \\ 0 & 0 & -3 & 3 \end{bmatrix} \begin{bmatrix} 0\\ u_{2} \\ u_{3} \\ u_{4} \end{bmatrix} = \begin{bmatrix} 5/3 \\ 10/3 \\ 10/3 \\ 155/3 \end{bmatrix} \end{equation*}

Que pode ser simplificado para:

(92)   \begin{equation*} 1\times 10^9 \times \begin{bmatrix} -6 & -3 & 0 \\ -3 & 6 & -3 \\ 0 & -3 & 3 \end{bmatrix} \begin{bmatrix} u_{2} \\ u_{3} \\ u_{4} \end{bmatrix} = \begin{bmatrix} 10/3 \\ 10/3 \\ 155/3 \end{bmatrix} \end{equation*}

Resolvendo-o, temos o seguinte resultado:

(93)   \begin{equation*} u_{2}=1,94\hspace{5} [10^{-8}m]; \hspace{20} u_{3}=3,78\hspace{5} [10^{-8}m];\hspace{20} u_{4}=5,50\hspace{5} [10^{-8}m] \end{equation*}

Vamos agora gerar a função de aproximação por elementos finitos û(x).

Função de aproximação de elementos finitos

A função de aproximação por elementos finitos desse problema será da seguinte forma:

(94)   \begin{equation*} û(x) = 0N_{1}(x) + 1,94N_{2}(x) + 3,78N_{3}(x) + 5,50N_{4}(x) \hspace{5} [10^{-8}m] \end{equation*}

Substituindo as funções de aproximação que foram definidas para o problema, tem-se por fim:

(95)   \begin{equation*} û(x) = \left\{ \begin{array}{ll} 5,82x \hspace{5} [10^{-8}m] \hspace{40} para \hspace{10} 0\leq x\leq 1/3\\ 0,1+5,52x \hspace{5} [10^{-8}m]\hspace{40} para \hspace{10} 1/3<x\leq 2/3 \\ 0,34+5,16x\hspace{5} [10^{-8}m] \hspace{40} para \hspace{10} 2/3<x\leq 1 \end{array}\right \end{equation*}

A solução exata para esse problema é a equação apresentada abaixo:

(96)   \begin{equation*} u(x) = -\frac{x^2}{2}+6x\hspace{5} [10^{-8}m] \end{equation*}

Se plotássemos os gráficos da função u(x) e de sua aproximação teríamos um resultado semelhante ao do exemplo anterior.

Pronto, mais um exemplo de aplicação do Método dos Elementos Finitos concluído.

Livro sobre elementos finitos

Há muitas referências de qualidade para se estudar o método dos Elementos Finitos. Dentre elas, indico o livro em língua portuguesa do Dr. Fernando Luiz B. Ribeiro, professor da disciplina de Elementos Finitos da UFRJ por décadas. Esse livro traz teoria e prática numa dosagem didaticamente confortável. Ele pode ser obtido na loja virtual da Amazon através do seguinte endereço: introdução ao Método dos Elementos Finitos.

Considerações finais sobre o Método dos Elementos Finitos

Apesar de termos atualmente uma gama de softwares especializados em análises numéricas envolvendo o Método dos Elementos Finitos, é imprescindível ter pelo menos os conhecimentos básicos do método para tomar decisões assertivas na etapa de modelagem e discretização do modelo e na etapa final de crítica aos resultados encontrados. Finalizo enfatizando dramaticamente: software de elementos finitos na mão de desconhecedores do método, apesar de toda a robustez que possa ter, será uma caixa preta que receberá números como input e devolverá números como output.

Deixe um comentário

O seu endereço de e-mail não será publicado. Campos obrigatórios são marcados com *