Código Computacional Avançado para Análise de Fenômenos Físicos de Campo Escalar da Engenharia

RESUMO

As equações diferenciais governam inúmeros fenômenos físicos presentes em diversos campos de conhecimento da física-matemática. Uma vez que os modelos são compreendidos pela teoria de campo escalar, possibilitam a descrição de uma gama de aplicações práticas da engenharia, transitando de problemas na área da acústica, térmica, mecânica dos sólidos e fluidos, eletromagnetismo, dentre outras. Neste cenário, o vigente artigo concentra-se em apresentar o módulo específico do programa computacional, desenvolvido em Matlab, denominado NASEN, destinado à solução da equação de campo escalar por meio do método dos elementos finitos (MEF). Para tanto, a discretização espacial é realizada com base nos procedimentos matemáticos do método de Galerkin e a discretização temporal é realizado com base no método de integração de Newmark e o método-θ. Estuda-se problemas de diferentes naturezas, como problemas parabólicos não lineares, hiperbólicos e de autovalor. A metodologia empregada para solução dos problemas não lineares tem como base as estratégias incrementais e iterativas de Newton-Raphson. A avaliação de desempenho do programa computacional é realizada com as soluções analíticas, numéricas ou experimentais disponíveis na literatura. Em síntese, os resultados obtidos mostram-se satisfatórios, indicando que o código computacional é capaz de prever adequadamente os comportamentos físicos dos problemas propostos.

Palavras-chave:
elementos finitos; teoria de campo escalar; programa computacional; física-matemática

ABSTRACT

Differential equations govern innumerable physical phenomena present in different fields of knowledge of physics-mathematical. Once the models are understood by the scalar field theory, they enable the description of a range of practical engineering applications, moving through problems in the area of ??acoustics, thermal, solids and fluids mechanics, electromagnetism, among others. In this scenario, the current article focuses on presenting the specific module of the computer program, developed in Matlab, called NASEN, aimed at solving the scalar field equation using the finite element method (FEM). Therefore, the spatial discretization is performed based on the mathematical procedures of the Galerkin method and the temporal discretization is performed based on the Newmark integration method and the A-method. Problems of different nature are studied, such as nonlinear parabolic, hyperbolic and eigenvalue problems. The methodology used to solve nonlinear problems is based on the incremental and iterative strategies of Newton-Raphson. The performance evaluation of the computer program is carried out with the analytical, numerical or experimental solutions available in the literature. In summary, the results obtained are satisfactory, indicating that the computational code is capable of adequately predicting the physical behavior of the proposed problems.

Keywords:
finite elements; scalar field theory; computer program; physics-mathematics

INTRODUÇÃO

A modelagem matemática-numérica de problemas físicos é um pilar que movimenta um vasto volume de pesquisa científica destinada à análise do comportamento de diferentes fenômenos práticos da engenharia, fornecendo conceitos, parâmetros e uma sensibilidade física acerca do fenômeno. Dentre os problemas estudados pela engenharia e áreas afins, uma parcela desses se enquadram na teoria de campo escalar. Essa teoria foi desenvolvida por volta do inicio do século XIX, visando integrar diversos fenômenos físicos presentes na natureza compreendidos pelo mesmo modelo matemático. O interesse visível no desenvolvimento dessa teoria era construir um modelo baseado em ideias de linha de fluxos, campo e potencial, sendo adicional aos modelos clássico de partículas ou mecanicistas ^{2), (28}. A aceitação da teoria de campo escalar direciona uma grande unificação entre diferentes áreas de conhecimento, possibilitando, por meio do mesmo modelo matemático, a solução de problemas de condução de calor ou elétrica, difusão de massa, propagação de ondas planas e cisalhantes, escoamento potencial, torção de barras prismáticas, deflexão de membranas, problemas de autovalor e autovetor, dentre outras aplicações da física-matemática ^{42), (44}.

As equações diferenciais parciais ou ordinárias representativas dos fenômenos físicos e com base na teoria de campo escalar, podem apresentar diversas complexidades físicas, geométricas e matemáticas, impossibilitando a solução por meio de técnicas analíticas. Sendo assim, a solução é encaminhada com base no emprego de métodos numéricos, sendo eles recorrentemente utilizados em aplicações de engenharia. Usualmente, destaca-se quatro grandes técnicas aproximativas: o método das diferenças finitas (MDF), volumes finitos (MVF), elementos finitos (MEF) e elementos de contorno (MEC). Cada método apresenta particularidades na sua implementação e formulação matemática ²⁷. Todavia, em modo geral, ao fim dos procedimentos e desenvolvidos numéricos de cada técnica, deve-se resolver um sistema algébrico efetivo, fornecendo, por exemplo, o campo de temperatura, deslocamento ou pressão. Nesse estudo, aplica-se os conceitos e formulações com base nos procedimentos numéricos de elementos finitos de Galekin ⁴².

Inúmeras pesquisas de natureza numérica, desenvolvidas com base no MDF, MEF, MVF e MEC, são elaboradas visando solucionar os diferentes fenômenos físicos governados pela teoria de campo escalar, por exemplo, os problemas térmicos não lineares da área de estruturas em condição de incêndio ^{9), (31), (39), (40), (46}, os fenômenos descritos pela equação hiperbólica ^{16), (24), (25}, os problemas de autovalor ^{3), (11), (22}, os problemas estruturais de torção de barras prismáticas ^{15), (35}, e os problemas clássicos puramente difusivos ^{2), (13}.

Nesse contexto, o presente artigo visa apresentar o módulo computacional específico do programa desenvolvido NASEN (Numerical Analysis System for Engineering), que contempla a solução de problemas governados pela teoria de campo escalar para diferentes contextos físicos. Essa ferramenta de solução desenvolvida é estruturada em módulos computacionais, capazes de realizar análises numéricas avançadas de inúmeros problemas clássicos de engenharia. O sistema computacional enquadra-se na solução de problemas em diferentes áreas de conhecimento da engenharia, como na análise térmica e termomecânica de estruturas aço, concreto e mistas em condição de incêndio, análise estrutural linear e não linear geométrica, e análises preliminares de problemas de acústica/vibração ^{29), (31), (32), (33}. Os módulos computacionais foram implementados em ambiente Matlab ²⁶, pelo fato ser uma linguagem de programação de alto nível que permite utilizar funções e recursos avançados predefinidos no sistema, como solução de sistemas lineares, cálculos simbólicos e recursos gráficos.

Em particular, são analisadas três classes de problemas neste trabalho, os problemas descritos pelas equações de natureza parabólica e hiperbólica, e os problemas de autovalor. A avaliação de desempenho do módulo computacional desenvolvido para as soluções dos problemas físicos de engenharia, são realizadas com base nos resultados obtidos por meio das soluções analíticas, simulações numéricas e previsões experimentais encontradas na literatura.

2 PROCEDIMENTOS FÍSICOS-MATEMÁTICOS E NUMÉRICOS GERAIS

Os problemas físicos de engenharia, direcionados pela teoria de campo escalar, são governados por uma equação diferencial parcial (EDP) caracterizada pela presença dos operadores diferenciais espaciais de 2ª-ordem e temporais de 1ª e 2ª-ordem, conforme posto na Eq. (2.1).

$-{\nabla }^T\mathbf{D}\nabla u+c_{0}u+c_1\dot{u}+c_2\ddot{u}=f(x,y,t)$ (2.1)

Onde u é um escalar associado a um campo físico, f é o termo fonte, $\dot{u}$ e ü são as derivadas temporais de 1ª e 2ª-ordem, respectivamente. Além disso, c ₀, c ₁ e c ₂ são constantes físicas e D é conhecida como tensor de difusão do material, uma vez considerado o caso bidimensional e meio isotrópico, vale a Eq. (2.2).

$\mathbf{D}=\left[\begin{array}{c}\begin{array}{cc}k& 0\\ 0& k\end{array}\end{array}\right]=k\mathbf{I}$ (2.2)

Na Eq. (2.1), caso as propriedades físicas do meio dependam do valor do potencial u, a equação diferencial, nesta condição, é dita sendo não linear. Na área de estruturas em condição de incêndio, as propriedades físicas dos materiais, como o aço, concreto e madeira, variam com o aumento de temperatura ^{7), (8}. Usualmente, problemas dessa natureza são resolvidos com base em uma estratégia incremental-iterativa.

Paralelamente, para a solução da equação diferencial de campo escalar, apresentada na Eq. (2.1), deve-se impor ao sistema condições de contorno e iniciais apropriadas ao problema. Considerando um domínio genérico Ω associado à fronteira Γ, simplificadamente, o contorno pode ser composto por parcelas representativas das condições essenciais e naturais atuantes nas bordas do domínio, conforme mostra as Eqs. (2.3), (2.4) e (2.5).

$u=\overline{u}\text{em}{\Gamma }_u$ (2.3)

$q=\overline{q}\text{em}{\Gamma }_q$ (2.4)

$q_{\iota }=\alpha (u-u_{\infty })\text{em}{\Gamma }_{\iota }$ (2.5)

Onde ū representa o potencial prescrito, $\overline{q}$ é o fluxo prescrito, q _ι é um fluxo linearizado, u _∞ é um potencial conhecido no problema, como uma temperatura externa de um fluido, e α é o coeficiente efetivo, que pode variar com o tempo e com campo potencial primal. Sendo assim, o fluxo normal q _n é computado, com base na lei de Fourier, pela sentença matemática a seguir:

$-\mathbf{D}\nabla u=q_n=\overline{q}+\alpha \left(u-u_{\infty }\right)$ (2.6)

Pelo fato do modelo diferencial geral de campo escalar considerar os operadores diferenciais temporais de 1ª e 2ª-ordem, o problema necessita de duas condições iniciais, ou seja, impõem-se ao sistema $u\left(0\right)=u_0$ e $\dot{u}\left(0\right)=v_0$ no instante inicial.

A solução numérica da equação diferencial do problema físico de interesse, dada pela Eq. (2.1), parte das premissas do método de resíduos ponderados, que busca encontrar uma solução aproximada de modo que o resíduo da equação diferencial seja mínimo ⁴². Desta forma, para facilitar a compreensão conceitual da aplicação das ideias do método dos resíduos ponderados e dos processos numéricos de elementos finitos, a Fig. 1 apresenta uma estrutura global utilizada nas etapas de cálculo.

Figura 1:
Estrutura geral dos procedimentos numéricos de elementos finitos.

Como é possível observar no fluxograma da Fig. 1, o desenvolvimento do modelo de elementos finitos para problemas dependentes do tempo envolve dois estágios principais. A primeira, chamada de semidiscretização, é desenvolver a forma fraca das equações sobre um elemento e buscar aproximações espaciais das variáveis dependentes do problema. O resultado final desta etapa é um conjunto de equações diferenciais ordinárias no tempo entre os valores nodais das variáveis dependentes. Para problemas transientes, o segundo estágio consiste na aproximação temporal da equações diferenciais ordinárias.

Desta maneira, primeiramente, realiza-se a discretização do espaço, com base nas sentenças de resíduos ponderados e nos procedimentos do cálculo diferencial-integral, apresenta-se, na (2.7), a formulação variacional fraca do problema físico.

$\underset{\Omega }{\int }{\nabla }^{T}w(\mathbf{D}\nabla u)d\Omega +\underset{\Omega }{\int }{c}_{0}uwd\Omega +\underset{\Omega }{\int }{c}_1\dot{u}wd\Omega +\underset{\Omega }{\int }{c}_2\ddot{u}wd\Omega -\underset{\Omega }{\int }fwd\Omega -\underset{\Gamma }{\int }w{q}_{n}d\Gamma =0$ (2.7)

Onde w é denominada usualmente como função auxiliar ou peso, e q _n é o fluxo normal no contorno (Eq. (2.6)). As sentenças integrais, mostradas na Eq. (2.7), são válidas em um domínio contínuo. Realizando as aproximações de elementos finitos com base na funções de forma N _i e os valores nodais $u_i^e$, conforme posto na Eq. (2.8).

$u(x,y)=\sum _{i=1}^n{N}_i{u}_i^e=\mathbf{N}{\left\{\mathbf{u}\right\}}^e$ (2.8)

Definindo-se o vetor de funções de forma N e o vetor de potenciais nodais u por meio das expressões apresentadas na Eq. (2.9).

$\mathbf{N}=\left[\begin{array}{cccc}{N}_1& N_2& \cdots & N_n\end{array}\right]\mathbf{u}={\left[\begin{array}{cccc}{u}_1& u_2& \cdots & u_n\end{array}\right]}^T$ (2.9)

Satisfazendo as características do método dos elementos finitos, as funções de forma são definidas de modo a assumirem o valor unitário nó i e zero nos demais nós do elemento, garantindo, por exemplo, que as matrizes sejam esparsas ⁴².

O sistema de equações resultante local é proveniente da substituição da aproximação de elementos finitos, conforme apresentado na Eq. (2.8), na formulação variacional fraca do problema físico, resultando na Eq. (2.18), escrita em notação matricial fechada.

${\left[\mathbf{M}\right]}^e{\left\{\ddot{\mathbf{u}}\right\}}^e+{\left[\mathbf{C}\right]}^e{\left\{\dot{\mathbf{u}}\right\}}^e+{\left[\mathbf{K}\right]}^e{\left\{\mathbf{u}\right\}}^e={\left\{{\mathbf{f}}_{\Gamma }\right\}}^e+{\left\{{\mathbf{f}}_{\Omega }\right\}}^e={\left\{\mathbf{F}\right\}}^e$ (2.10)

Onde:

${\left[\mathbf{M}\right]}^e=\underset{\Omega }{\int }{\mathbf{N}}^T\left[c_2\right]\mathbf{N}d\Omega$ (2.11)

${\left[\mathbf{C}\right]}^e=\underset{\Omega }{\int }{\mathbf{N}}^T\left[c_1\right]\mathbf{N}d\Omega$ (2.12)

${\left[\mathbf{K}\right]}^e={\left[{\mathbf{K}}_d\right]}^e+{\left[{\mathbf{K}}_r\right]}^e=\underset{\Omega }{\int }\left[{\mathbf{B}}^T\mathbf{DB}+{\mathbf{N}}^T\left[c_0\right]\mathbf{N}\right]d\Omega$ (2.13)

${\left\{{\mathbf{f}}_{\Gamma }\right\}}^e=\underset{\Gamma }{\int }{\mathbf{N}}^{\mathrm{T}}{q}_{n}d\Gamma$ (2.14)

${\left\{{\mathbf{f}}_{\Omega }\right\}}^e=\underset{\Omega }{\int }{\mathbf{N}}^{\mathrm{T}}fd\Omega$ (2.15)

A Eq. (2.10) apresenta um sistema geral de solução da equação diferencial de campo escalar com operadores temporais de 1ª e 2ª-ordem. Fisicamente, as Eqs.(2.11) e (2.12) representam os efeitos inerciais e de dissipação de energia, enquanto, na Eq. (2.13), a matriz K _d e K _r correspondem aos efeitos físicos difusivos e reativos. A Eq. (2.14) representa as contribuições relativas das ações dos fluxos normais ao contorno e a Eq. (2.15) representa a geração de energia ou a fonte de energia interna no domínio. Em aplicações de mecânica estrutural, a matriz M é denominada de matriz de massa, C é a matriz de amortecimento do sistema e K é usualmente denotada de matriz de rigidez. Além disso, a Eq. (2.10) pode ser interpretada dizendo que a forças externas do sistema são equilibradas pela combinação das forças inerciais, de amortecimento e elásticas.

O vetor gradiente pode ser escrito em conjunto com a função de interpolação de elementos finitos $\nabla u=\nabla \mathbf{N}\tilde{u}=\mathbf{B}\tilde{u}$, conforme expresso na Eq. (2.16), na forma expandida.

$\mathbf{B}=\left[\begin{array}{cccc}\frac{\partial N_1}{\partial x}& \frac{\partial N_2}{\partial x}& & \frac{\partial N_n}{\partial x}\\ \frac{\partial N_1}{\partial y}& \frac{\partial N_2}{\partial y}& \stackrel{...}{}& \frac{\partial N_n}{\partial y}\end{array}\right]$ (2.16)

As matrizes globais do sistema são computadas por meio do assemblamento das matrizes locais de cada elemento contido na malha bidimensional de elementos finitos, conforme apresenta a Eq. (2.17).

$\mathbf{C}=\sum _{e=1}^{ne}{\left[\mathbf{C}\right]}^e\mathbf{M}=\sum _{e=1}^{ne}{\left[\mathbf{M}\right]}^e\mathbf{K}=\sum _{e=1}^{ne}{\left[\mathbf{K}\right]}^e\mathbf{F}=\sum _{e=1}^{ne}{\left[\mathbf{F}\right]}^e$ (2.17)

A partir da Eq. (2.17), o sistema de equações resultante global é posto a seguir na Eq. (2.18).

$\left[\mathbf{M}\right]\left\{\ddot{\mathbf{u}}\right\}+\left[\mathbf{C}\right]\left\{\dot{\mathbf{u}}\right\}+\left[\mathbf{K}\right]\left\{\mathbf{u}\right\}=\left\{{\mathbf{f}}_{\Gamma }\right\}+\left\{{\mathbf{f}}_{\Omega }\right\}=\left\{\mathbf{F}\right\}$ (2.18)

Com base nos processos apresentados na Fig. 1, após a discretização dos termos espaciais presentes na equação de governo do problema físico (Eq. (2.1)), deve-se realizar o tratamento numérico das parcelas temporais de 1ª e 2ª-ordem.

2.1 Equação Parabólica

A partir do sistema resultante, apresentado na Eq. (2.18), considere o caso com $\mathbf{M}=0$. O sistema é reduzido ao modelo representativo de fenômenos físicos de condução de calor transiente, em algum caso específico da dinâmica dos fluidos ou problemas similares definidos em uma região bidimensional, conforme posta na Eq. (2.19).

$\left[\mathbf{C}\right]\left\{\dot{\mathbf{u}}\right\}+\left[\mathbf{K}\right]\left\{\mathbf{u}\right\}=\left\{\mathbf{F}\right\}$ (2.19)

A condição inicial requerida ao sistema é definida pela imposição do potencial no instante inicial, ${\left\{\mathbf{u}\right\}}_{t=0}=\left\{{\mathbf{u}}_0\right\}$. O operador diferencial temporal de 1ª-ordem $\left(\dot{\mathbf{u}}\right)$ é aproximado por meio do método-θ^{29), (42}, em que uma média ponderada da derivada de tempo de uma variável dependente é aproximada em duas etapas de tempo consecutivas por interpolação linear dos valores da variável nas duas etapas, conforme mostra a Eq. (2.20),

$\left(1-\theta \right){\left\{\dot{\mathbf{u}}\right\}}_t+\alpha {\left\{\dot{\mathbf{u}}\right\}}_{t+1}=\frac{{\left\{\mathbf{u}\right\}}_{t+1}-{\left\{\mathbf{u}\right\}}_t}{\Delta t}$ (2.20)

Onde ∆t é o passo de tempo e o parâmetro θ varia entre 0 e 1. Com base na Eq. (2.20), pode-se escrever as Eqs. (2.21) e (2.22), no passo de tempo $t+\theta$.

${\left\{\mathbf{u}\right\}}_{t+1}={\left\{\mathbf{u}\right\}}_t+\Delta t{\left\{\dot{\mathbf{u}}\right\}}_{t+\theta }$ (2.21)

${\left\{\dot{\mathbf{u}}\right\}}_{t+\theta }=\left(1-\theta \right){\left\{\dot{\mathbf{u}}\right\}}_t+\theta {\left\{\dot{\mathbf{u}}\right\}}_{t+\theta }$ (2.22)

O parâmetro θ pode assumir variações distintas, acarretando em esquemas de integração no tempo diferentes. Caso $\theta =1$, 1/2 ou 2/3, têm-se o esquema Euler implícito, Crank-Nicolson ou Galerkin, respectivamente, sendo esses incondicionalmente estáveis. Por sua vez, para $\theta =0$, tem-se o esquema de Euler explícito, um método condicionalmente estável ^{17), (42), (43}. No presente trabalho, utiliza-se o esquema de Galerkin ($\theta =2/3$) para todas as simulações dos problemas dessa natureza.

Aplicando as aproximações, apresentadas nas Eqs. (2.20), (2.21) e (2.22), na Eq. (2.19), após reorganização dos termos, chega-se no sistema algébrico efetivo para cada passo de tempo, conforme apresenta a Eq. (2.23).

$\left[\mathbf{C}+\theta \Delta t\mathbf{K}\right]{\left\{\mathbf{u}\right\}}_{t+1}=\left[\mathbf{C}-\left(1-\theta \right)\Delta t\mathbf{K}\right]{\left\{\mathbf{u}\right\}}_t+\left(1-\theta \right)\Delta t{\mathbf{F}}_t+\gamma \Delta t{\mathbf{F}}_{t+1}$ (2.23)

Assume-se, na Eq. (2.23), que as matrizes são independentes do tempo, contudo, essas podem depender do campo potencial, tornando o problema não linear, devendo-se aplicar um método incremental-iterativo para se obter a solução do problema físico. Neste trabalho aplica-se o método de Newton-Raphson ^{1), (42}. De forma geral, a Fig. 2 apresenta a estrutural geral de programação de problemas físicos descritos pela equação parabólica não linear.

Figura 2:
Estrutura geral do código computacional para problemas parabólicos não lineares.

Como pode-se notar, na linha 13 da Fig. 2, a solução numérica deve satisfazer a um critério de convergência em cada passo de tempo. Para a análise térmica de estruturas em condição de incêndio (solução não linear), adota-se uma tolerância de 0,1 ^{31), (39}. Além disso, as propriedades físicas do meio c ₁ e k dependem do campo primal u, e devem ser atualizadas em cada iteração.

2.2 Equação Hiperbólica

Os fenômenos físicos descritos pela equação hiperbólica são definidos, a partir da Eq. (2.18), considerando $\mathbf{C}=0$, resultando na equação apresentada a seguir:

$\left[\mathbf{M}\right]\left\{\ddot{\mathbf{u}}\right\}+\left[\mathbf{K}\right]\left\{\mathbf{u}\right\}=\left\{\mathbf{F}\right\}$ (2.24)

A natureza hiperbólica é ditada por $c_2\ne 0$, não por $c_1=0$; de fato, em sistemas estruturais dinâmicos (problema vetorial), c ₁ e c ₁ podem ser diferentes de zero ³⁰. Todavia, os problemas tratados no presente trabalho são governados pela teoria de campo escalar, adotando $c_1=0$.

A equação hiperbólica, apresentada na Eq. (2.24), pode ser reduzida a um sistema de equações algébricas por meio da aproximação do operador diferencial temporal de 2ª-ordem. Para tanto, aplica-se o esquema de integração de tempo de Newmark ³⁶, conforme posto nas Eqs. (2.25) e (2.26).

${\left\{\dot{\mathbf{u}}\right\}}_{t+1}={\left\{\dot{\mathbf{u}}\right\}}_t+\left[\left(1-\delta \right){\left\{\ddot{\mathbf{u}}\right\}}_t+\delta {\left\{\ddot{\mathbf{u}}\right\}}_t\right]\Delta t$ (2.25)

${\left\{\mathbf{u}\right\}}_{t+1}={\left\{\mathbf{u}\right\}}_t+\Delta t{\left\{\dot{\mathbf{u}}\right\}}_t+\left[\left(1/2-\delta \right){\left\{\dot{\mathbf{u}}\right\}}_t+\alpha {\left\{\dot{\mathbf{u}}\right\}}_{t+1}\right]\Delta t^2$ (2.26)

Onde α e δ são parâmetros que podem ser determinados para obter precisão e estabilidade de integração. Quando $\alpha =1/6$ e $\delta =1/2$, tem-se o método da aceleração linear, condicionalmente estável. Por sua vez, para $\alpha =1/4$ e $\delta =1/2$, define-se o método da aceleração média constante, um esquema incondicionalmente estável ¹. Dessa forma, com base nas Eqs. (2.24), (2.25) e (2.26), em cada passo de tempo, deve-ser resolver o seguinte sistema de equações:

$\left[a_0\mathbf{M}+\mathbf{K}\right]{\left\{\mathbf{u}\right\}}_{t+1}={\mathbf{F}}_{t+1}+\mathbf{M}\left[a_0{\left\{\mathbf{u}\right\}}_t+a_2{\left\{\dot{\mathbf{u}}\right\}}_t+a_3{\left\{\ddot{\mathbf{u}}\right\}}_t\right]$ (2.27)

Onde $a_0=1/\alpha \Delta t^2,a_2=1/\alpha \Delta t$ e $a_3=1/2\alpha -1$. A lógica utilizada na implementação dos procedimentos de solução utilizando o método de Newmark é mostrado na Fig. 3.

Figura 3:
Solução step-by-step usando método de integração de Newmark.

Como pode-se notar, na Fig. 3, após determinar o vetor u no instante $t+1$, deve-se atualizar os vetores $\dot{\mathbf{u}}$ e ü, representativos, por exemplo, do vetor de velocidade e aceleração do sistema, respectivamente. Além disso, essa rotina de cálculo é restringida somente ao caso do sistema com $\mathbf{C}=0$.

2.3 Análise de Autovalor

O problema de autovalor, associado a Eq. (2.24), busca encontrar $\mathbf{u}=\overline{\mathbf{u}}{e}^{-\lambda t}$, de modo que as condições iniciais e de contorno sejam homogêneas e $f=0$. As matrizes K e M são independentes de u. Sendo assim, a expressão característica do problema de autovalor é mostrada na Eq. (2.28).

$\left(\left[\mathbf{K}\right]-\lambda \left[\mathbf{M}\right]\right)\left\{\overline{\mathbf{u}}\right\}=\left\{{\mathbf{f}}_{\Gamma }\right\}$ (2.28)

Onde $\lambda ={\omega }^2$ é chamado de autovalor e ū é denominado de autovetor. Para um problema escalar de membrana, ω representa as frequências naturais de vibração do sistema. O número de autovalores do sistema discreto (Eq. (2.28)) do problema é igual ao número de valores nodais desconhecidos de u na malha.

2.4 Elemento triangular linear

O módulo de análise bidimensional de problemas físicos do programa NASEN, governados pela teoria de campo escalar, apresenta duas classes de elementos finitos planos para discretização do domínio bidimensional ²⁹, primeiramente, (i) a família de elementos triangulares: T3, T6 - Elemento triangular de 3 e 6 nós, respectivamente, e (ii) a família de elementos quadriláteros: Q4, Q8 e Q9 Elemento quadrilateral de 4, 8 e 9 nós, respectivamente.

No presente trabalho, as simulações são realizadas somente com elementos do tipo triangular de 3 nós. Esse tipo de elemento é caracterizado pela aproximação linear entre nós, conforme esquematizado na Fig. 4. Sendo assim, as funções de forma do elemento triangular, escritas em termos de variáveis locais (ξ e η), são expressas na Eq. (2.29).

$N_1(\xi ,\eta )=\xi N_2(\xi ,\eta )=\eta N_3(\xi ,\eta )=1-\xi -\eta$ (2.29)

Figura 4:
Características do elemento bidimensional do tipo triangular linear.

Essas funções de formas (Eq. (2.29)) são associadas com as variáveis globais por meio do processo de mapeamento ^{1), (29), (42}. O elemento triangular linear apresenta uma formulação matemática simples e versátil para a discretização de diferentes tipos de geometrias utilizadas na prática da engenharia, possibilitando computar, analiticamente, sem esforços matemáticos, as matrizes elementares (ou locais) do problema, conforme posto nas Eqs. (2.30), (2.31), (2.32) e (2.33). Todavia, para os elementos de alta-ordem (T6, Q4, Q8 e Q9), emprega-se a integração numérica de Gauss ¹.

${\left[\mathbf{C}\right]}^e=c_1\frac{A_e}{12}\left[\begin{array}{ccc}2& 1& 1\\ 1& 2& 1\\ 1& 1& 2\end{array}\right]{\left[\mathbf{M}\right]}^e=\frac{c_2}{c_1}{\left[\mathbf{C}\right]}^e{\left\{\mathbf{Q}\right\}}^e=f\frac{A_e}{3}\left\{\begin{array}{c}1\\ 1\\ 1\end{array}\right\}$ (2.30)

${\left[\mathbf{K}\right]}^e={\left[{\mathbf{K}}_d\right]}^e+{\left[{\mathbf{K}}_r\right]}^e=\frac{k}{4A_{\mathbf{e}}}\left[\begin{array}{ccc}{k}_{11}& k_{12}& k_{13}\\ k_{21}& k_{22}& k_{23}\\ k_{31}& k_{32}& k_{33}\end{array}\right]+\frac{c_0}{c_1}{\left[\mathbf{C}\right]}^e$ (2.31)

Onde A _e é a área de um elemento finito triangular e as componentes da matriz de difusão $(k_{ij},i,j=1,2,3)$, em função das coordenadas globais de cada elemento, são apresentadas na Eq. (2.32).

$$\begin{gathered} k_{11}={\left(y_2-y_3\right)}^2+{\left(x_3-x_2\right)}^2 \\ {k}_{22}={\left(y_3-y_1\right)}^2+{\left(x_1-x_3\right)}^2 \\ {k}_{33}={\left(y_1-y_2\right)}^2+{\left(x_2-x_1\right)}^2 \\ {k}_{12}=k_{21}=\left(y_2-y_3\right)\left(y_3-y_1\right)+\left(x_3-x_2\right)\left(x_1-x_3\right) \\ {k}_{13}=k_{31}=\left(y_2-y_3\right)\left(y_1-y_2\right)+\left(x_3-x_2\right)\left(x_2-x_1\right) \\ {k}_{23}=k_{32}=\left(y_3-y_1\right)\left(y_1-y_2\right)+\left(x_1-x_3\right)\left(x_2-x_1\right) \end{gathered}$$ (2.32)

Para computar as matrizes/vetores correspondentes as condições de contorno do problema, utiliza-se um elemento unidimensional linear de comprimento L, representativo de uma das faces do elemento triangular sob ação da condição de contorno do tipo natural (Eq. (2.6)), conforme esquematizado na Fig. 4. Matematicamente, as matrizes e os vetores locais, característicos dos efeitos do fluxo linearizado e prescrito, são mostrados na Eq. (2.33).

$K_{\iota }=\frac{\alpha L}{6}\left[\begin{array}{cc}2& 1\\ 1& 2\end{array}\right]\left\{{\mathbf{F}}_{\iota }\right\}=\frac{\alpha L{u}_{\infty }}{2}\left\{\begin{array}{c}1\\ 1\end{array}\right\}\left\{{\mathbf{F}}_q\right\}=\overline{q}\frac{L}{2}\left\{\begin{array}{c}1\\ 1\end{array}\right\}$ (2.33)

Para garantir a correta aderência entre o elemento unidimensional de dois nós e o elemento triangular bidimensional de 3 nós, deve-se assegurar a conectividade dos nós coincidentes entre os elementos.

3 EXEMPLOS NUMÉRICOS DE APLICAÇÃO

Nesta seção busca-se realizar uma investigação numérica acerca dos problemas físicos regidos pela equação de campo escalar. No presente estudo, as aplicações tratam da análise do problema de condução de calor, do comportamento térmico não linear de estruturas em condição de incêndio, da solução do problema de autovalor e da avaliação do comportamento de membranas elásticas. Em todos os testes numéricos são utilizados elementos triangulares lineares.

Dessa forma, são analisados problemas físicos bidimensionais governados pela equação parabólica e hiperbólica, bem como as análises dos problemas de autovalor. Os problemas estudados apresentam características físicas e geométricas bem definidas, sendo, então, empregado neste trabalho uma malha de elementos finitos única para cada caso-teste estudado, uma vez que em todos os casos as avaliações de desempenho são aferidas com soluções de referências de natureza analítica, numérica e/ou experimental. Mesmo que em alguns casos-testes analisados apresentem efeitos não lineares ou condições não usuais, os problemas não sofrem variações significativas com sucessivos refinamentos das malhas numéricas, como pode ser visto nos inúmeros testes numéricos e trabalhos desenvolvidos, utilizando o programa NASEN, para análise de problemas físicos similares em domínios bidimensionais ^{29), (31), (32), (33), (34}.

Para quantificar a precisão das respostas aproximadas (x _aprox) obtidas com o programa NASEN. Para tanto, define-se a norma || e || e o erro relativo |∆|, conforme expressas nas equações a seguir:

$||\mathbf{e}||=\sqrt{\sum _{i=1}^{\text{n}°\text{pontos}}{\left({\mathbf{x}}_{\text{ref}}^i-{\mathbf{x}}_{\text{aprox}}^i\right)}^2}$ (3.1)

$\left|\Delta \right|=100\times \left|{\mathbf{x}}_{\text{ref}}-{\mathbf{x}}_{\text{aprox}}\right|/\left|{\mathbf{x}}_{\text{ref}}\right|$ (3.2)

Neste trabalho, a Eq. (3.1) é aplicada para verificar a precisão de toda amostra de dados, enquanto, a Eq. (3.2) é utilizada para avaliar a precisão entre cada elemento da amostra de dados. Ressalta-se que, em alguns casos não é possível aplicar essas equações, uma vez que os dados numéricos das soluções de referências (x _ref) dos casos estudados não estão disponíveis para realizar os tratamentos adequados da amostra de dados, sendo então realizada somente uma avaliação gráfica qualitativa entre as curvas.

3.1 Problemas de Natureza Parabólica

3.1.1 Placa infinita mantida a temperatura prescrita na fronteira

No primeiro caso-teste analisado, considera-se a equação diferencial de natureza parabólica, sendo que o potencial genérico u representa o campo de temperatura T . Sendo assim, esse exemplo trata-se de uma condução de calor em uma placa infinita com a condutividade térmica do material variando linearmente com a temperatura, respeitando a função $k=2+0,1T$, em W/m^◦ C. A capacidade calorífica é mantida constante e igual a $c_1=\rho c=8,0$, em Ws/m^3◦ C, onde c é o calor específico e ρ é a massa específica do material. O passo de tempo utilizado na simulação é igual 0,05 s e aplica-se o esquema de Galerkin ($\theta =2/3$).

Esse problema foi estudado em inúmeras pesquisas a fim de verificar eficiência de códigos computacionais, vide ^{4), (5), (38), (41}. O domínio é simplificadamente representado por um retângulo associado as condições de contorno do tipo essencial aplicadas nas fronteiras, conforme ilustrado na Fig. 5(a). A malha numérica utilizada na simulação computacional é tipo triangular linear (ver Fig. 5(b)).

Figura 5:
(a) Dimensões, características térmicas e condições de contorno da placa infinita e (b) malha de elementos finitos.

A face direita do domínio é mantida à temperatura constante e igual a 100^◦ C, enquanto a temperatura prescrita na face esquerda varia com tempo, assumindo um valor constante de 200^◦ C até 10 s e instantes seguintes, subitamente muda para uma temperatura constante de 100^◦ C. Em relação à condição inicial do problema, considera-se que a placa é mantida a uma temperatura de 100^◦ C.

Cabe ressaltar que esse problema apresenta duas características importantes e significava na reposta numérica. Primeiramente, a condutividade variar com a temperatura, tornando o problema não linear. Em segundo momento, a dificuldade no tratamento numérico devido à presença de uma descontinuidade na condição de contorno. Os resultados são apresentados graficamente por meio da distribuição de temperatura, ao longo do eixo x, da placa para os tempos de 10 e 13 s, conforme mostram as Figs. 6(a) e 6(b), respectivamente.

Figura 6:
Perfil de temperatura da placa infinita ao fim de (a) 10 e (b) 13 s.

Para avaliar quantitativamente os resultados, a Tabela 1 apresenta os erros relativos (Eq. (3.2)) medidos em relação à solução de Orivuori ³⁷. É possível observar que o programa desenvolvido NASEN apresenta baixos níveis de erros relação à solução de Orivuori ³⁷, enquanto, a solução de Vila Real ⁴⁵ apresenta uma ligeira diferença no comportamento, uma vez que o autor utiliza na sua simulação um esquema de discretização diferente ($\theta =1/2$) e um passo de tempo igual a 1,0 s, podendo acarretar essas possíveis diferenças.

Tabela 1:
Comparação entre os níveis de erro obtidos em cada solução.

	10 s		13 s
x (cm)	|∆0−1|(%)	|∆0−2|(%)	|∆0−1| (%)	|∆0−2| (%)
0	0,000	0,000	0,000	0,000
1	0,027	2,326	0,741	4,675
2	0,043	4,639	0,127	3,494
3	0,042	6,533	0,548	0,451
4	0,001	6,987	0,650	5,349
5	0,017	5,861	0,099	5,865
6	0,020	3,897	0,136	5,075
7	0,005	2,302	0,023	3,657
8	0,020	1,169	0,007	2,279
9	0,018	0,522	0,007	1,257
10	0,011	0,190	0,011	0,609

• índices: 0 Orivuori (1979); 1 NASEN; 2 Vila Real (1988)

3.1.2 Viga de concreto armado em situação de incêndio

O próximo caso-teste visa analisar um problema no contexto da engenharia de segurança contra incêndio. Desta forma, sendo um elemento frequentemente utilizado na construção civil, considera-se uma viga de concreto com armaduras de aço submetida à ação do incêndio, conforme mostra a Fig. 7(a). A viga quadrada de 200x200 mm com oito armaduras de diâmetro de 16 mm é discretizada com elementos triangulares lineares (ver Fig. 7(b)) e é exposta ao incêndio-padrão descrito pela curva ISO 834 ¹⁴, conforme posto na Eq. (3.3).

$T_g=T_{amb}+345\log (8t+1)$ (3.3)

Figura 7:
(a) Detalhamento e dimensões da viga de concreto armado, (b) malha numérica de elementos lineares triangulares e (c) distribuição de temperatura na seção transversal para 30 e 120 min de exposição ao fogo.

Onde t é o tempo de exposição ao incêndio, em minutos, T _amb é a temperatura ambiente, em graus Celsius, que usualmente é adotada 20^◦ C, e $T_g=u_{\infty }$ é a temperatura dos gases do ambiente aquecido, em graus Celsius.

No contexto de estruturas em situação de incêndio, o coeficiente α, contido na expressão matemática do fluxo linearizado de calor, conforme mostrado na Eq. (2.5), é caracterizado pelos efeitos combinados de convecção e radiação provocados pela ação do incêndio no elemento de concreto ³², conforme mostra a Eq. 3.4.

$\alpha ={\alpha }_c+{\alpha }_r={\alpha }_c+ϵ\sigma (T+T_g)(T^2+{T_g}^2)$ (3.4)

Onde α _c é o coeficiente de transferência de calor por convecção, adotado igual a 25 W/m^2oC, conforme as prescrições do EN 1992 ⁷, e α _r é o coeficiente de radiação. Em que σ é a constante de Stefan-Boltzmann igual a $5,6697·{10}^{-8}$ W/m²K⁴, ε é a emissividade, adotada igual a 0,7. Observa-se, na Eq. (3.4), que o coeficiente de radiação apresenta maior influência da temperatura em relação ao coeficiente de convecção, ou seja, após um certo tempo de exposição ao fogo, o efeito de radiação torna-se predominante.

Esse problema é caracterizado por uma natureza não linear por conta da variação das propriedades dos materiais com a temperatura e pela presença do fluxo de calor do tipo radiação. O comportamento das armaduras de aço e do concreto, em função do aumento de temperatura, seguem as recomendações do EN 1993 ⁸ e EN 1992 ⁷, conforme mostram as Figs. 8, 9 e 10.

Figura 8:
Massa específica do (a) concreto e do (b) aço em função da temperatura.

Figura 9:
Calor específico do (a) concreto e do (b) aço em função da temperatura.

Figura 10:
Condutividade do (a) concreto e do (b) aço em função da temperatura.

As análises são realizadas ao longo da linha de simetria vertical, onde são definidos os pontos 16, localizados, respectivamente, na profundidade 0 até 100 mm em relação à face inferior da viga, conforme esquematizado na Fig. 7(a). Os resultados numéricos obtidos com programa NASEN são comparados com dados experimentais fornecidos em Haksever e Anderberg ¹² e com as simulações numéricas realizadas no trabalho de Di Capua e Mari ⁶. Neste problema, analisase duas configurações físicas: (i) concreto sem umidade (ou concreto seco) e (ii) concreto com teor de umidade de U _m = 6%. A umidade do concreto é considerada, no modelo de elementos finitos, por meio da função que descreve o calor específico do concreto, conforme apresentado na Fig. 9(a). Para valores de umidade do concreto diferentes de 0, 1,5 e 3%, são realizadas interpolações lineares, conforme as recomendações do EN 1992 ⁷.

Para as duas configurações analisadas neste problema, os resultados são apresentados nas Figs. 11(a) e 11(b), onde é possível observar que os resultados obtidos com o código computacional NASEN são próximos as simulações computacionais realizadas por Di Capua e Mari ⁶. Além disso, os níveis de temperatura são maiores no ponto 1, uma vez que este está localizado na superfície em contato direto com o fogo, em contrapartida, o ponto 6, localizado no centro da viga, apresenta os menores níveis de temperatura, evidenciando a baixa condutividade do concreto, conforme pode ser visualizado graficamente na distribuição de temperatura apresentada na Fig. 7(c).

Figura 11:
Comparação entre os resultados obtidos com programa NASEN e os dados experimentais/numéricos da literatura, admitindo o concreto (a) sem e (b) com umidade, para o caso-teste da viga de concreto armado em situação de incêndio.

Em relação aos resultados experimentais ¹², o modelo com concreto seco fornece resultados preliminares satisfatórios, conforme mostra a Fig. 11(a). Quando se considera, no modelo numérico, o teor de umidade do concreto, a precisão dos resultados melhora em relação aos pontos 1, 2, 3 e 4, conforme apresentado na Fig. 11(b). Todavia, é possível observar que os pontos 5 e 6 apresentaram resultados com certa divergência na solução, esse comportamento pode ser por conta de desprezar a ocorrência da migração da umidade ao longo da seção transversal de concreto, uma vez que a umidade, na realidade, apresenta níveis mais altas nos pontos centrais da viga de concreto (como os pontos 5 e 6), conforme relatado na literatura ^{6), (19), (20), (21}.

3.2 Problemas de Natureza Hiperbólica

3.2.1 Membrana quadrada sob deslocamento inicial prescrito em todo o domínio

Neste exemplo de aplicação, analisa-se um problema físico descrito pela equação de natureza hiperbólica (Eq. (2.24)), caracterizado por uma membrana elástica retangular de largura a e comprimento b, e com os todos os lados da membrana engastados, conforme mostra a Fig. 12(a).

Figura 12:
(a) Características da membrana retangular e (b) malha de elementos finitos triangular.

Assume-se que a tração na membrana é igual a $k=12,5$ N/m e a densidade é igual $c_2=2,5$ kg/m². Além disso, a condição inicial acerca do campo de deslocamento da membrana, no instante $t=0$, é dada pela expressão matemática $u(x,y,0)=u_0=0,1(4x-x^2)(2y-y^2)$. Além disso, impõem-se ao sistema a condição $\dot{u}(x,y,0)=v_0=0$.

O objetivo dessa análise numérica é determinar o campo de deslocamento u(x,y,t) da membrana retangular, medido no ponto central (a/2, b/2), em função do tempo. A membrana é discretizada com elementos planos triangulares de 3 nós, conforme apresenta a Fig. 12(b). Para avaliação do desempenho do programa desenvolvido, utiliza-se a solução analítica, que pode ser encontrada na literatura clássica ^{10), (18}.

A Fig. 13(a) mostra a comparação dos resultados entre a solução aproximada e a exata do problema. O programa NASEN apresenta um comportamento semelhante a solução exata. Quanto à norma de diferença em relação à solução exata (Eq. (3.1)), foi obtido o valor 0,015, indicando o bom desempenho do programa na predição do deslocamento da membrana.

Figura 13:
(a) Comparação da deflexão central obtida usando o programa NASEN e a solução analítica de um membrana retangular com deflexão inicial, $∆t=0.025$; (b) características de estabilidade dos esquemas de aceleração média constante e aceleração linear, $∆t=0.25$.

Adicionalmente, na Fig. 13(b), apresenta uma investigação da estabilidade da solução numérica conforme o esquema de integração utilizado. Para tanto, utiliza-se o método da aceleração linear ($\delta =1/2$ e $\alpha =1/6$) e o método da aceleração média constante ($\delta =1/2$ e $\alpha =1/4$), conforme discutido na Sec. 2.2. Observa-se que a solução numérica, utilizando o esquema de aceleração média constante, apresenta uma resposta estabilizada, uma vez que este é classificado como incondicionalmente estável. Por sua vez, quando considera-se o método da aceleração linear, observa-se que para $t>1$ s, a solução apresenta pequenas oscilações, e em instantes seguintes, a perda total de estabilidade da solução. Desta forma, para métodos condicionalmente estáveis, deve-se obrigatoriamente atender a um critério de estabilidade, conforme pode ser visto em Reddy ⁴².

3.2.2 Membrana quadrada sob velocidade inicial prescrita sobre parte do domínio

A membrana quadrada, representada na Fig. 14(a), com o campo de velocidade inicial $v_0=V$ prescrito sobre o subdomínio Ω₀ e os deslocamentos nulos prescritos em toda a fronteira, é analisada neste exemplo. A solução analítica para este problema é dada por ¹⁸.

Figura 14:
(a) Características da membrana sob velocidade inicial prescrita sobre parte do domínio e (b) malha de elementos finitos.

Neste exemplo, por conta da descontinuidade da condição inicial, a malha de elementos finitos triangulares lineares foi refinada em torno do subdomínio Ω₀, conforme mostra a Fig. 14(b). A avaliação dos resultados numéricos obtidos com programa NASEN é realizada com a evolução do deslocamento, no ponto central da membrana, em função do tempo. Utiliza-se o método de Newmark com esquema de aceleração média constante ^36),(42, o passo de tempo é igual a 0,01 s.

Aplicando a Eq. (3.1), a norma de diferença em relação aos resultados analíticos, fornece o valor 0,0686. Além disso, na Fig. 15, os resultados obtidos mostram que a curva numérica apresenta um comportamento físico coerente com a solução analítica, indicando a boa precisão do programa desenvolvido na solução do problema estudado.

Figura 15:
Deslocamento no ponto central (a/2, a/2) da membrana quadrada sob campo de velocidade inicial sobre parte do domínio.

3.3 Problema de Autovalor

O último problema estudado visa realizar uma avaliação do desempenho do programa computacional desenvolvido, NASEN, para determinar os autovalores, com base na Eq. (2.28). Sendo assim, estuda-se, numericamente, dois casos-testes com condições de contorno diferentes, conforme representado nas Figs. 16(a) e 16(b), respectivamente. As dimensões e propriedades físicas dos casos são consideradas unitárias.

Figura 16:
Características e condições de contorno do (a) caso 1 e do (b) caso 2, e (c) malha triangular de elementos finitos.

A malha numérica utilizada para simulação, em ambos os casos, é formada por elementos triangulares lineares, totalizando 80 nós e 98 elementos, conforme mostra a Fig. 16(c). Além disso, para a verificação do desempenho do programa NASEN, são utilizados soluções de natureza numérica encontradas no trabalho de Loeffler et al.²³, com base no método dos elementos de contorno com integração direta (MECID) e no método de elementos finitos (MEF), além da solução analítica conhecida do caso 2.

A Tabela 2 apresenta os valores numéricos obtidos pelo programa NASEN e os respectivos valores absolutos de erro relativo (Eq. (3.2)) obtidos pela comparação entre cada método. Ao fim, pode-se verificar que o código desenvolvido no presente trabalho atingiu valores satisfatórios com erros percentuais diminutos. É possível observar os resultados entre NASEN e MEF são próximos, enquanto, quando comparados ao MECID, os resultados apresentam um aumento na ordem de grandeza do erro a medida que as frequências aumentam. Essa diferença entre os técnicas numéricas pode ser devido a diferentes fatores, por exemplo, a formulação e tratamentos matemáticos de cada método, a solução do sistema característico do problema de autovalor, a escolha das funções de interpolação, dentre outros.

Tabela 2:
Caso 1 Frequências naturais da membrana.

ω n	MEF1	MECID2	NASEN3	|∆1−3| (%)	|∆2−3| (%
1	3,143	3,144	3,142	0,0191	0,0509
2	4,448	4,449	4,446	0,0381	0,0605
3	6,292	6,288	6,290	0,0375	0,0261
4	7,036	7,035	7,033	0,0361	0,0219
5	7,043	7,042	7,040	0,0428	0,0287
6	8,917	8,919	8,913	0,0444	0,0668
7	9,451	9,454	9,447	0,0472	0,0789
8	9,969	9,945	9,965	0,0437	0,1975
9	9,969	9,984	9,965	0,0421	0,1923
10	11,371	11,335	11,366	0,0427	0,2747
11	11,401	11,423	11,396	0,0460	0,2385
12	12,624	12,647	12,618	0,0488	0,2306
13	13,020	12,970	13,014	0,0440	0,3413
14	13,023	13,071	13,018	0,0416	0,4086
15	13,425	13,420	13,419	0,0428	0,0056
16	14,151	14,038	14,144	0,0470	0,7576
17	14,153	14,255	14,146	0,0470	0,7622
18	15,810	15,833	15,803	0,0456	0,1908
19	15,826	15,886	15,818	0,0487	0,4262
20	15,902	15,913	15,895	0,0453	0,1143

A Tabela 3 apresenta os resultados obtidos para o segundo caso, determinando as frequências naturais da membrana. A solução exata do problema é igual ${\omega }_w=\pi \sqrt{n^2+m^2}$. É possível verificar que os resultados numéricos do programa NASEN são satisfatórios, com baixo níveis de erro, evidenciando o bom desempenho do programa na resolução do problema.

Tabela 3:
Caso 2 Frequências naturais da membrana.

ω mn		Exata0	MECID1	|∆0−1| (%)	NASEN2	|∆0−2| (%)
		m	n
1	1	4,4429	4,4460	0,0702	4,4463	0,0771
1	2	7,0248	7,0240	0,0116	7,0341	0,1324
2	1	7,0248	7,0240	0,0116	7,0393	0,2066
2	2	8,8858	8,8640	0,2450	8,9131	0,3071
1	3	9,9346	9,9450	0,1048	9,9646	0,3026
3	1	9,9346	9,9670	0,3263	9,9648	0,3037
2	3	11,3272	11,3310	0,0338	11,3673	0,3541
3	2	11,3272	11,3310	0,0338	11,3946	0,5952
1	4	12,9531	13,0170	0,4932	13,0151	0,4785
4	1	12,9531	13,3980	3,4346	13,0166	0,4904
3	3	13,3286	13,0170	2,3382	13,4193	0,6802
2	4	14,0496	14,0900	0,2873	14,1443	0,6739
4	2	14,0496	14,0910	0,2945	14,1461	0,6866
3	4	15,7080	15,8590	0,9615	15,8149	0,6806
4	3	15,7080	15,8590	0,9615	15,8924	1,1741
1	5	16,0190	16,1110	0,5741	16,1336	0,7150
5	1	16,0190	16,2110	1,1983	16,1336	0,7154
2	5	16,9180	16,9880	0,4138	17,0707	0,9024
5	2	16,9180	16,9880	0,4138	17,0772	0,9412
4	4	17,7715	18,0400	1,5107	17,9776	1,1593

4 CONCLUSÃO

No presente estudo são apresentados testes numéricos acerca da verificação de um módulo computacional específico desenvolvido do programa NASEN para a solução de problemas físicos descritos pela teoria de campo escalar. Os problemas tratados são divididos em três classes, os problemas parabólicos, hiperbólicos e de autovalor. O desenvolvido do modelo de elementos finitos para problemas dependentes do tempo e de autovalor, são divididos em dois estágios. A discretização do espaço, baseada nos procedimentos numéricos clássicos do método dos resíduos ponderados e as propostas de elementos finitos, e a discretização do tempo, fundamentada nas estratégias do método de Newmark e do método-θ.

As análises dos casos-testes realizados neste trabalho, buscam estudar diferentes configurações físicas e geométricas. Foi tratado, numericamente, problemas considerando uma condição essencial descontinua variando com tempo, propriedades dos materiais em função do campo potencial, condições de contorno do tipo natural não lineares, além das condições iniciais não usuais. Em linhas gerais, em todos os casos-testes estudados, os resultados obtidos pelo programa NASEN mostram-se satisfatório, direcionando o bom desempenho do código nas aplicações propostas.

Em pesquisas futuras, visa-se ampliar as aplicabilidades do programa computacional NASEN para outras áreas de conhecimento e adicionar novas funcionalidades ao programa, por exemplo, implementação de novos esquemas de integração do tempo, tratamentos sobre os efeitos provocados pelo ar enclausurado em cavidades nas estruturas sob altas temperaturas e aprimoramentos gráficos dos módulos.

REFERÊNCIAS

• ¹
  K.J. Bathe. “Finite element procedures”. Klaus-Jurgen Bathe (2006).
• ²
  A. Bulcão. “Formulação do Método dos Elementos de Contorno com Dupla Reciprocidade Usando Elementos de Ordem Superior Aplicada a Problemas de Campo Escalar Generalizado”. Dissertação de mestrado, Universidade Federal do Espírito Santo, Espírito Santo (1999).
• ³
  D.H. Choi & W.J. Hoefer. The finite-difference-time-domain method and its application to eigenvalue problems. IEEE Transactions on Microwave Theory and Techniques, 34(12) (1986), 1464-1470.
• ⁴
  F. Damjanić & D.R.J. Owen. Practical considerations for thermal transient finite element analysis using isoparametric elements. Nuclear Engineering and Design, 69(1) (1982), 109-126.
• ⁵
  F.B. Damjanic. “Reinforced concrete failure prediction under both static and transient conditions”. Ph.D. thesis, University College of Swansea (1983).
• ⁶
  D. Di Capua & A.R. Mari. Nonlinear analysis of reinforced concrete cross-sections exposed to fire. Fire Safety Journal, 42(2) (2007), 139-149.
• ⁷
  EN 1992-1-2. Eurocode 2: Design of concrete structures-Part 1-2: General rules-structural fire design.European Committee for Standardization, (2004).
• ⁸
  EN 1993-1-2. Eurocode 3: Design of steel structures part 1-2: General rules structural fire design. European Committee for Standardization, (2005).
• ⁹
  A. Farcas & D. Lesnic. The boundary-element method for the determination of a heat source dependent on one variable. Journal of Engineering Mathematics, 54(4) (2006), 375-388.
• ¹⁰
  A.M. Ferri & W. Pihua. “Boundary Element Methods In Engineering And Sciences”, volume 4. World Scientific (2010).
• ¹¹
  J. Gary. Computing eigenvalues of ordinary differential equations by finite differences. Mathematics of Computation, 19(91) (1965), 365-379.
• ¹²
  A. Haksever, Y. Anderberg, A. Anderberg & Y. Anderberg. Comparison between measured and computed structural response of some reinforced concrete columns in fire. Fire Safety Journal , 4(4) (1982), 293-297.
• ¹³
  F. Hermeline. A finite volume method for the approximation of diffusion operators on distorted meshes. Journal of computational Physics, 160(2) (2000), 481-499.
• ¹⁴
  ISO 834. Fire resistance tests-elements of building construction. International Organization for Standardization, Geneva, Switzerland, (1999).
• ¹⁵
  C. Jog & I.S. Mokashi. A finite element method for the Saint-Venant torsion and bending problems for prismatic beams. Computers & Structures, 135 (2014), 62-72.
• ¹⁶
  M.K. Kadalbajoo, A. Kumar & L.P. Tripathi. A radial basis functions based finite differences method for wave equation with an integral condition. Applied Mathematics and Computation, 253 (2015), 8-16.
• ¹⁷
  W. Ko¨hler & J. Pittr. Calculation of transient temperature fields with finite elements in space and time dimensions. International journal for numerical methods in engineering, 8(3) (1974), 625-631.
• ¹⁸
  E. Kreyszig. “Advanced Mathematics for Engineering”, volume 6. John Wiley, New York (1988).
• ¹⁹
  T.T. Lie. Fire resistance of circular steel columns filled with bar-reinforced concrete. Journal of structural engineering, 120(5) (1994), 1489-1509.
• ²⁰
  T.T. Lie & B. Celikkol. Method to calculate the fire resistance of circular reinforced concrete columns.Materials Journal, 88(1) (1991), 84-91.
• ²¹
  T.T. Lie & R.J. Erwin. Method to calculate the fire resistance of reinforced concrete columns with rectangular cross section. ACI structural Journal, 90(1) (1993), 52-60.
• ²²
  C. Loeffler, A. Frossard & L. Lara. Testing complete and compact radial basis functions for solution of eigenvalue problems using the boundary element method with direct integration. International Journal for Computational Methods in Engineering Science and Mechanics, 19(2) (2018), 117-128.
• ²³
  C.F. Loeffler, W.J. Mansur, H.M. Barcelos & A. Bulcão. Solving Helmholtz problems with the boundary element method using direct radial basis function interpolation. Engineering Analysis with Boundary Elements, 61 (2015), 218-225.
• ²⁴
  R. Lo¨hner, K. Morgan & O.C. Zienkiewicz. The solution of non-linear hyperbolic equation systems by the finite element method. International Journal for Numerical Methods in Fluids, 4(11) (1984), 1043-1063.
• ²⁵
  K.J. Marfurt. Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave equations. Geophysics, 49(5) (1984), 533-549.
• ²⁶
  MATLAB. “version 8.5.0 (R2015a)”. The MathWorks Inc., Natick, Massachusetts (2015).
• ²⁷
  J.L. Meek. “Computer methods in structural analysis”. CRC Press (2017).
• ²⁸
  P. Moon & D.E. Spencer. “Field theory handbook, including coordinate systems, differential equations and their solutions”. Springer (1971).
• ²⁹
  N.S. Neves. “Modelo computacional avançado para análise de estruturas sob ação de gradientes térmicos”. Dissertação de mestrado, Universidade Federal do Espírito Santo, Vitória, ES (2019).
• ³⁰
  N.S. Neves. “Uma introdução ao uso da integral de Duhamel em sistemas dinâmicos estruturais”, volume 2. Estudos (inter) multidisciplinares nas engenharias, Ponta Grossa PR (2019), chapter 18, p. 190-201.
• ³¹
  N.S. Neves, M.S. Azevedo, C.B. Barcelos Filho, V.P. Silva & I. Pierin. Estudo térmico de pilares mistos de aço e concreto de seção circular em situação de incêndio. REA Revista da Estrutura de Aço, 9(2) (2020), 122-140.
• ³²
  N.S. Neves , M.S. Azevedo , R.S. Camargo & V.P. Pinheiro. Análise térmica bidimensional de perfil de aço sujeita a elevadas temperaturas. In “Anais do X Encontro Científico de Física Aplicada” . Blucher Proceedings (2019), p. 71-73.
• ³³
  N.S. Neves , M.S. Azevedo , R.S. Camargo & V.P. Pinheiro. Application of the finite element method for two-dimensional thermal analysis of steel-concrete composite structures in fire. In “18th Brazilian Congress of Thermal Sciences and Engineering ENCIT”. ABCM (2020), p. 1-6.
• ³⁴
  N.S. Neves & V.P. Pinheiro. Determinação do campo de tensao efetiva em um problema vetorial elástico-linear via método de elementos finitos. In “Anais do X Encontro Cientıfico de Fısica Aplicada”. Blucher Proceedings (2019), p. 85-87.
• ³⁵
  N.S. Neves , V.P. Pinheiro, D.C. Moura, C.F. Loeffler & N.L. Bodart. Uma estratégia didática para engenharia estrutural baseada na análise numérica de barras sob torção. In “As Diversidades de Debates na Pesquisa em Matemática”, volume 1 of 1. ATENA Editora, Ponta Grossa (PR) (2019), chapter 18, p. 206-218.
• ³⁶
  N.M. Newmark. A method of computation for structural dynamics. Journal of the engineering mechanics division, 85(3) (1959), 67-94.
• ³⁷
  S. Orivuori. Efficient method for solution of nonlinear heat conduction problems. IJNME, 14(10) (1979), 1461-1476.
• ³⁸
  D. Owen & F. Damjanić. The stability of numerical time integration techniques for transient thermal problems with special reference to reduced integration effects. In “Num. Meth. in Thermal Problems Conf. Proc.”, volume 2. Pineridge Press (1981), p. 487-505.
• ³⁹
  I. Pierin , V. Silva & H. La Rovere. Thermal analysis of two-dimensional structures in fire. Revista IBRACON de Estruturas e Materiais, 8(1) (2015), 25-36.
• ⁴⁰
  D. Pires, R.C. Barros, P.A.S. Rocha & R.A.d.M. Silveira. Thermal analysis of steel-concrete composite cross sections via CS-ASA/FA. REM-International Engineering Journal, 71(2) (2018), 149-157.
• ⁴¹
  G.A. Ramirez & J.T. Oden. Finite element technique applied to heat conduction in solids with temperature dependent thermal conductivity. International Journal for Numerical Methods in Engineering, 7(3) (1973), 345-355.
• ⁴²
  J.N. Reddy. An introduction to the finite element method. New York, 27 (1993).
• ⁴³
  J.B.C. Silva, E.C. Romão & L.F.M. de Moura. A comparison of time discretization methods in the solution of a parabolic equation. In “7th Brazilian Conference on Dynamics, Control and Applications” (2008).
• ⁴⁴
  A.N. Tikhonov & A.A. Samarskii. “Equations of mathematical physics”. Courier Corporation (2013).
• ⁴⁵
  P.M.M. Vila Real. “Modelação por elementos finitos do comportamento térmico e termo-elástico de sólidos sujeitos a elevados gradientes térmicos”. Ph.D. thesis, Tese de doutoramento, Universidade do Porto (1988).
• ⁴⁶
  M.X. Xiong, Y. Wang & J.R. Liew. Evaluation on thermal behavior of concrete-filled steel tubular columns based on modified finite difference method. Advances in Structural Engineering, 19(5) (2016), 746-761.

Datas de Publicação

• Publicação nesta coleção
  04 Abr 2022
• Data do Fascículo
  Jan-Mar 2022

Histórico

• Recebido
  04 Set 2020
• Aceito
  16 Jul 2021