Elastic dynamics analysis in steel tubular towers for horizontal axis wind turbines

Latest revision as of 13:55, 30 November 2020

Resumo

Este artigo tem por objetivo a análise dinâmica elástica de uma torre tubular de aço para aerogerador de eixo horizontal utilizando elementos finitos de barra (código próprio) e elementos finitos de casca e sólidos (modelagem feita no ANSYS). Para tal, apresenta-se uma fundamentação teórica que consiste no desenvolvimento da equação de movimento a partir da equação de Lagrange. Em seguida, mostra-se o método de superposição modal para resolução da equação de movimento matricial e obtenção da resposta dinâmica da torre. Na segunda parte deste artigo, apresenta-se o desenvolvimento para montagem das matrizes de massa e de amortecimento da torre modelada com código próprio em elementos finitos de barra. Na parte de resultados, inicialmente, expõem-se as matrizes de rigidez, de massa e de amortecimento para certo nível de discretização da torre modelada com elementos finitos de barra. Em seguida, foram mostrados os resultados da análise modal feita tanto para o modelo com elementos finitos de casca, com e sem base flexível, quanto para o modelo de elementos finitos de barra; em que se apresentam os modos e as frequências de vibração da torre em todos os modelos. Por fim, foram analisadas as respostas da torre: na direção do fluxo de vento, quando submetida por um vetor de forças ressonantes, com o 1º modo de vibração desta, que representam a parcela flutuante do vento; e, na direção transversal ao fluxo de vento em decorrência do desprendimento cadenciado de vórtices de von Kárman. Os resultados da análise modal feita, tanto para o modelo com elementos finitos de casca, com e sem base flexível, quanto para o modelo de elementos finitos de barra mostram-se similares, principalmente com relação ao 1º modo de vibração da torre. Portanto, o modelo constituído por Elementos Finitos (EF) de barra pode ser utilizado como representativo do comportamento dinâmico da torre quando esta é submetida às excitações ressonantes à sua frequência fundamental.

Palavras-chave: Torre do aerogerador, dinâmica estrutural, análise modal, análise transiente

Abstract

This article aims to obtain the dynamic elastic analysis of a steel tube tower for horizontal axis wind turbine using finite elements of bar (own code) and finite elements of bark and solids (modeling in ANSYS). For this, a theoretical foundation is presented that consists in the development of the movement equation from the Lagrange equation. Next, the modal overlay method is presented for solving the matrix motion equation and for obtaining the dynamic response of the tower. In the second part of this article, the development is presented for the assembly of the mass and damping matrices of the tower modeled with its own code in finite elements of the bar. In the part of results, initially, the rigidity, mass and damping matrices are exposed for a certain level of discretization of the modeled tower with finite elements of bar. Then, the results of the modal analysis were done for both the finite element shell model, with and without flexible base, and the finite element model of the bar, in which the modes and the frequencies of vibration of the tower in all the models are presented. Finally, the responses of the tower were analyzed: in the direction of the wind flow, when submitted by a vector of resonant forces to its 1st mode of vibration, representing the floating part of the wind; and, in the transverse direction to the wind flow due to the von Kárman cadenced vortex detachment. The results of the modal analysis made for both the finite element shell model, with and without flexible base, and the finite element model of the bar are similar, especially in the first mode of vibration of the tower. Therefore, the model represented by finite elements (FE) of bar can be used as representative of the dynamic behavior of the tower when it is subjected to the resonant excitations at its fundamental frequency.

Keywords: Wind turbine tower, structural dynamics, modal analysis, transient analysis

1. Introdução

No início do século XXI, teve-se um crescimento acelerado na implantação de aerogeradores, onshore e offshore, de porte crescente com torres cada vez mais altas [1]. O desenvolvimento, o comércio e a instalação de aerogeradores no mundo se desenvolveram rapidamente, de forma que a geração de energia a partir de termoelétricas, usina nucleares e hidrelétricas tenha sido complementada e/ou substituída pela produção daqueles equipamentos.

A geração de energia elétrica por meio de turbinas eólicas constitui uma alternativa para diversos níveis de demanda no Brasil. As pequenas centrais podem suprir pequenas localidades distantes da rede de distribuição; já às centrais de grande porte têm potencial para atender uma significativa parcela do Sistema Interligado Nacional (SIN) com importantes ganhos. Especificamente no Nordeste brasileiro (especialmente dos estados da Paraíba, Rio Grande do Norte, Ceará, Piauí e Pernambuco), o desenvolvimento da produção de energia eólica se deu de maneira promissora nos últimos anos, pois diversas usinas eólicas estão em operação e em fase de implantação, fazendo com que a geração de energia elétrica de origem eólica tenha crescido exponencialmente na última década [2].

Aliado ao exposto acima, a evolução do tamanho dos aerogeradores, cada vez mais pesados e potentes, torna necessária a instalação destes equipamentos sob a ação de ventos mais intensos e contínuos, fazendo com que as dimensões das torres destes aerogeradores estejam sendo incrementadas [3]. Particularmente, a altura da torre é um parâmetro essencial para captação de ventos estáveis de grande altura; entretanto, o custo da torre, que pode superar 20% do custo total do gerador eólico [4], faz com que o aumento de altura represente uma desvantagem. Além disto, o transporte, a montagem e a posta em operação da torre tornam-se mais custosos.

Com o aumento da altura das torres para aerogeradores os efeitos de vibração a que os componentes dos aerogeradores ficam submetidos são amplamente aumentados, de forma que, nas últimas décadas, pesquisadores [5-11] e fabricantes de componentes para aerogeradores têm buscado soluções para avaliar vibrações oriundas do funcionamento dos componentes mecânicos e dos efeitos eólicos e sísmicos aos quais estas estruturas são submetidas. Tais ações podem provocar vibrações excessivas, causando danos e comprometendo a integridade estrutural dos componentes do aerogerador. Mesmo em situações em que estas estruturas possam suportar a ação das cargas dinâmicas, sem obrigatoriamente sofrer danos estruturais, não se pode desprezar os efeitos de fadiga dos materiais constituintes. É importante mencionar que existem dados técnicos relevantes desenvolvidos pelas fabricantes de turbinas eólicas que auxiliariam às pesquisas acadêmicas. No entanto, tais informações são preservadas, evitando a divulgação de conteúdos que podem ser acessados por concorrentes, ainda que exista a patente. Os estudos oriundos da academia são, portanto, a única fonte pública possível para fomentar contribuições a serem amplamente compartilhadas.

Assim, a crescente utilização dos aerogeradores modernos e sua importância para a geração de energia elétrica reforçam a necessidade de estudos que busquem melhorar o desempenho de seus componentes. Especificamente, a torre que dá suporte ao aerogerador é um elemento essencial para seu funcionamento. Fatores como a esbeltez e a rigidez da torre influenciam diretamente no comportamento estrutural, na estabilidade e na dinâmica do aerogerador. Por isto, este trabalho propõe uma contribuição adicional analisando-se o comportamento dinâmico de uma torre de aço tubular com 120 m de altura destinada a suportar um aerogerador com potência nominal de 3,2 MW.

2. Fundamentação teórica

2.1. Equação de movimento em elementos finitos

A equação de movimento da torre pode ser estabelecida a partir da equação de Lagrange, que consiste em um equacionamento formal baseado em princípios energéticos (via cálculo variacional). Para utilização desta equação será considerado que a estrutura foi discretizada, de forma que os deslocamentos da estrutura passam a ser descritos em função dos deslocamentos nodais dos diversos elementos finitos que a compõem. A equação de Lagrange na forma matricial é dada por:

(1)

na qual: ${\textstyle t}$ é a variável temporal; ${\textstyle {E}_{c}}$ é o escalar da energia cinética da estrutura; ${\textstyle U}$ é o escalar da energia potencial de deformação da estrutura; ${\textstyle \Omega }$ é o escalar do potencial das forças externas conservativas aplicadas fora dos nós da estrutura; ${\textstyle \left\{{\boldsymbol {D}}\right\}=}$ $\left\{{\boldsymbol {D}}(t)\right\}$ é o vetor de funções temporais de deslocamentos dos graus de liberdade da estrutura no referencial global de coordenadas; ${\textstyle \left\{{\overset {\cdot }{\boldsymbol {D}}}\right\}=}$ $\left\{{\overset {\cdot }{\boldsymbol {D}}}(t)\right\}$ é o vetor de funções temporais de velocidades dos graus de liberdade da estrutura no referencial global de coordenadas; ${\textstyle \left\{{\boldsymbol {F}}\right\}=}$ $\left\{{\boldsymbol {F}}(t)\right\}$ é o vetor de funções temporais de forças aplicadas nos nós da estrutura no referencial global de coordenadas; e, ${\textstyle \left\{{\boldsymbol {F}}_{NC}\right\}=}$ $\left\{{\boldsymbol {F}}_{NC}(t)\right\}$ é o vetor de funções temporais de forças não conservativas no referencial global de coordenadas.

Para uma estrutura com comportamento linear físico, a energia cinética ${\textstyle {E}_{c}}$ é dada por:

(2)

a energia potencial de deformação ${\textstyle U}$ é dada por:

(3)

e o potencial das forças externas conservativas é:

(4)

nas quais: ${\textstyle \left[{\boldsymbol {M}}_{E}\right]}$ é a matriz de massa consistente da estrutura no referencial global de coordenadas, a qual recebe essa denominação de consistente, uma vez que, as matrizes de massa dos elementos finitos são obtidas utilizando-se as mesmas funções de forma ${\textstyle \left[{\boldsymbol {N}}\right]}$ que gerarão a matriz de rigidez dos elementos finitos [12-14]; ${\textstyle \left[{\boldsymbol {K}}_{E}\right]}$ é a matriz de rigidez tangencial da estrutura no referencial global de coordenadas; e, ${\textstyle \left\{{\boldsymbol {R}}_{E}^{F}\right\}=}$ $\left\{{\boldsymbol {R}}_{E}^{F}(t)\right\}$ é o vetor de funções temporais das reações de extremo fixo da estrutura no referencial global de coordenadas.

Ademais, as forças não conservativas podem ser internas ou externas, ou seja, resultantes da deformação da estrutura (internas) ou de forças diretamente aplicadas (externas). Assim, estabelece-se:

(5)

em que: ${\textstyle \left\{{\boldsymbol {F}}_{NE}\right\}=\left\{{\boldsymbol {F}}_{NE}(t)\right\}}$ é o vetor de funções temporais das forças não conservativas externas no referencial global de coordenadas; ${\textstyle \left[{\boldsymbol {C}}_{E}\right]}$ é matriz de amortecimento da estrutura no referencial global de coordenadas, que resulta das forças não conservativas internas à estrutura.

Substituindo as Eqs.(2), (3), (4) e (5) em (1), obtém-se a equação de movimento da estrutura:

(6)

na qual, ${\textstyle \left\{{\boldsymbol {F}}_{\boldsymbol {E}}\right\}=}$ $\left\{{\boldsymbol {F}}_{\boldsymbol {E}}(t)\right\}$ é o vetor de funções temporais das forças nodais da estrutura no referencial global de coordenadas, dado por:

(7)

2.2. Método da superposição modal

Uma forma amplamente utilizada para obter a resposta dinâmica de uma estrutura é feita a partir da superposição das respostas de seus modos de vibração. Para isso, entretanto, é necessário que se consiga desacoplar as equações de movimento, de maneira a reduzir o problema à resolução de n (número de graus de liberdade dinâmicos do problema) equações de movimento com apenas 1 grau de liberdade, que terão suas soluções superpostas.

A aplicação do método da superposição é realizada por meio da seguinte marcha de cálculo:

1) Diagonalizar a matriz de massa ${\textstyle \left[{\boldsymbol {M}}_{E}\right]}$ :

Calcula-se a matriz de autoversores ${\textstyle \left[{\boldsymbol {\varphi }}\right]}$ da matriz de massa ${\textstyle \left[{\boldsymbol {M}}_{E}\right]}$ para diagonalizá-la mediante uma transformação linear ortogonal, uma vez que ${\textstyle \left[{\boldsymbol {M}}_{E}\right]}$ é simétrica:

(8)

na qual, ${\textstyle \left[{\boldsymbol {M}}^{\mathit {\boldsymbol {\ast }}}\right]}$ é a matriz de massa diagonalizada da estrutura, que contém os autovalores de ${\textstyle \left[{\boldsymbol {M}}_{E}\right]}$ . Com esta transformação:

(9)

a equação matricial de movimento não amortecido passa a ter a seguinte forma:

(10)

a qual está escrita no referencial que diagonaliza a matriz de massa. Sendo:

(11)

2) Calcular a matriz dinâmica inversa ${\textstyle \left[{\boldsymbol {DI}}\right]}$ :

(12)

a qual é calculada desta maneira de forma a ficar simétrica e reduzir o custo computacional do sistema analisado. Com mais uma transformação:

(13)

a equação matricial de movimento passa a ter a seguinte forma:

(14)

a qual está escrita no referencial que transforma a matriz de massa original em uma matriz identidade ${\textstyle \left[{\boldsymbol {I}}\right]}$ de ordem n. Sendo:

(15)

3) Diagonalizar a matriz dinâmica inversa ${\textstyle \left[{\boldsymbol {DI}}\right]}$ :

Calcula-se a matriz de autoversores ${\textstyle \left[{\overline {\boldsymbol {\psi }}}\right]}$ da matriz de dinâmica inversa ${\textstyle \left[{\boldsymbol {DI}}\right]}$ para diagonalizá-la mediante uma transformação linear ortogonal, uma vez que ${\textstyle \left[{\boldsymbol {DI}}\right]}$ é simétrica:

(16)

na qual, ${\textstyle \left[{\boldsymbol {\omega }}^{2}\right]}$ é uma matriz diagonal que contém as frequências angulares quadradas referentes aos modos de vibração da estrutura (que contém os autovalores de ${\textstyle \left[{\boldsymbol {DI}}\right]}$ ), obtida a partir da seguinte transformação ortogonal:

(17)

Assim, a matriz ${\textstyle \left[{\boldsymbol {\omega }}^{2}\right]}$ contém a solução do problema padrão de autovalores e autovetores dado pela Eq.(14). É importante que a matriz de frequências angulares ${\textstyle \left[{\boldsymbol {\omega }}^{2}\right]}$ seja arranjada de maneira a terem-se as frequências organizadas da menor para maior, assim como a matriz de transformação ${\textstyle \left[{\overline {\boldsymbol {\psi }}}\right]}$ seja correspondentemente organizada. Como são calculadas n frequências, é necessário que as frequências menores, que em geral são as mais importantes na contribuição modal, sejam visualizadas prioritariamente.

4) Calcular a matriz modal ponderada ${\textstyle \left[{\boldsymbol {\Phi }}\right]}$ :

(18)

5) Calcular o vetor de forças no referencial denominado generalizado, em que a matriz ${\textstyle \left[{\boldsymbol {DI}}\right]}$ fica diagonalizada:

(19)

6) Reescrever as equações de movimento (agora desacopladas) no referencial denominado generalizado:

(20)

ou na forma indicial:

(21)

na qual, j é o índice correspondente ao j-ésimo modo de vibração.

7) Solucionar o conjunto de equações diferenciais ordinárias (EDOs) considerando-se as condições iniciais do problema em questão:

(22)

nas quais, ${\textstyle \left\{{\overline {\boldsymbol {q}}}_{0}\right\}}$ e ${\textstyle \left\{{\boldsymbol {D}}_{0}\right\}}$ são os vetores de deslocamentos iniciais nos referenciais generalizado e original, respectivamente; e, ${\textstyle \left\{{\overset {\cdot }{\overline {\boldsymbol {q}}}}_{0}\right\}}$ e ${\textstyle \left\{{\overset {\cdot }{\boldsymbol {D}}}_{0}\right\}}$ são os vetores de velocidades iniciais nos referenciais generalizado e original, respectivamente.

8) Calcular a resposta da estrutura:

Após a resolução das EDOs, retorna-se ao referencial original aplicando a combinação modal, conforme:

(23)

Sendo a matriz modal ${\textstyle \left[{\overline {\boldsymbol {\Phi }}}\right]}$ dada por:

(24)

em que os modos de vibração são ortogonais entre si.

9) Calcular a matriz modal normalizada da estrutura:

É necessário normalizar a matriz modal para que toda informação referente aos modos de vibração provenham das funções de tempo ${\textstyle {q(t)}_{j}}$ que ponderam os modos de vibração. Logo, as formas modais normalizadas em notação indicial são dadas por:

(25)

em que, i representa o i-ésimo grau de liberdade dinâmico e ${\textstyle {{\eta }^{0}}_{j}}$ o módulo do j-ésimo modo de vibração:

(26)

assim, a superposição modal é feita em um novo referencial generalizado ${\textstyle {q}_{j}}$ :

(27)

da seguinte forma:

(28)

ou na seguinte forma, em que percebe-se claramente a superposição dos modos para obtenção da resposta da estrutura, de forma de cada modo de vibração tem uma contribuição (ponderação) para resposta total da estrutura:

(29)

na qual, ${\textstyle {\left\{{\boldsymbol {\Phi }}^{0}\right\}}_{j}}$ é o j-ésimo modo de vibração.

2.3. Matriz de amortecimento

Em geral, a magnitude das forças de amortecimento é inferior às forças de inércia e de rigidez. Entretanto, é difícil de especificar as forças de fricção desenvolvidas nos componentes estruturais, sendo mais significativo especificar as razões de amortecimento modal para construir a matriz de amortecimento da estrutura [14-16]. Assim, para que o método da superposição modal possa ser utilizado, a matriz de amortecimento deve atender a condição de ortogonalidade, ou seja, deve ser possível diagonalizar a matriz de amortecimento no mesmo referencial em que as matrizes de massa e de rigidez da estrutura também ficam diagonalizadas. E com isso, continua a ser possível desacoplar as equações de movimento dinâmico da estrutura.

Uma forma classicamente utilizada para considerar o amortecimento é obtida a partir do intitulado amortecimento proporcional ou amortecimento de Rayleigh, conforme:

(30)

na qual, os coeficientes ${\textstyle {a}_{0}}$ (em ${\textstyle {s}^{-1}}$ no SI) e ${\textstyle {a}_{1}}$ (em ${\textstyle s}$ no SI) podem ser determinados a partir de duas razões de amortecimento ( ${\textstyle {\xi }_{i}}$ e ${\textstyle {\xi }_{k}}$ ) relacionadas a duas frequências angulares de vibração distintas ( ${\textstyle {\omega }_{i}}$ e ${\textstyle {\omega }_{k}}$ ), montando-se o seguinte sistema (no referencial generalizado):

\left\{{\begin{matrix}{{\left\{{\boldsymbol {\Phi }}\right\}}_{i}}^{T}\left[{\boldsymbol {C}}_{E}\right]{\left\{{\boldsymbol {\Phi }}\right\}}_{i}={a}_{0}+{a}_{1}{{\omega }_{i}}^{2}=2{\,\xi }_{i}{\,\omega }_{i}\\{{\left\{{\boldsymbol {\Phi }}\right\}}_{k}}^{T}\left[{\boldsymbol {C}}_{E}\right]{\left\{{\boldsymbol {\Phi }}\right\}}_{k}={a}_{0}+{a}_{1}{{\omega }_{k}}^{2}=2{\,\xi }_{k}{\,\omega }_{k}\end{matrix}}\right.

(31)

donde, calcula-se os coeficientes por:

\left\{{\begin{matrix}{a}_{0}\\{a}_{1}\end{matrix}}\right\}=2{\left[{\begin{matrix}{\frac {1}{{\omega }_{i}}}&{\omega }_{i}\\{\frac {1}{{\omega }_{k}}}&{\omega }_{k}\end{matrix}}\right]}^{-1}\left\{{\begin{matrix}{\xi }_{i}\\{\xi }_{k}\end{matrix}}\right\}

(32)

então, com os coeficientes ${\textstyle {a}_{0}}$ e ${\textstyle {a}_{1}}$ calculados, determina-se a matriz de amortecimento ${\textstyle \left[{\boldsymbol {C}}_{E}\right]}$ e todas as demais razões de amortecimento ${\textstyle {\xi }_{j}}$ do j-ésimo modo de vibração da estrutura. Logo, incluindo a parcela do amortecimento à Eq.(21), obtém-se:

(33)

que pode ser resolvida para cada j-ésimo modo de vibração.

3. Modelo do projeto da torre

Foi considerada uma torre tubular de aço S355J2 [17] que dá suporte a um aerogerador no padrão SWT-3.2-113 [18], conforme características especificadas na Tabela 1.

Tabela 1. Dados do padrão do aerogerador utilizado [18]

TIPO DE PARÂMETRO
Classe segundo IEC (International Electrotechnical Commission)	IIA
Potência nominal (MW)	3,2
Diâmetro do rotor (m)	113,0
Comprimento da pá (m)	55,0
Área varrida pelo rotor (m²)	10000
Altura do cubo do rotor (m)	79,5 – 142,0 (usou-se 122,5 m)
Regulação de potência	Ângulo de passo regulado
Energia elétrica produzida anualmente a 8,5 m/s	14402 MWh
Peso da nacele (tf)	78
Peso do rotor (tf)	67

O modelo de torre tubular de aço analisado neste trabalho foi pautado no projeto elaborado em [19], apresentado na Figura 1.

Figura 1. Esquema do modelo de projeto da torre utilizado para análises (sem escala)

4. Modelagem da torre

Foi implementado no software Mathcad 14 os parâmetros do modelo de elementos finitos de barra da torre. Foram utilizados elementos finitos de barra com 8 graus de liberdade (4 graus de liberdade por nó: translações axiais e transversais, rotação flexional e rotação torcional) para modelagem da torre feita com código próprio.

No Sistema Global de Coordenadas (SGC), o equilíbrio dinâmico da torre pode ser estabelecido pelas condições de equilíbrio das ações em cada nó da estrutura, com relação a cada grau de liberdade, devido às reações elásticas dos elementos que concorrem neste nó, bem como as possíveis ações externas aplicadas diretamente aos nós. Assim, a expressão de equilíbrio de momentos de flexão é dada por:

${}{\lambda }^{(j)}{\ddot {\theta }}_{j-1}+\left({{\beta }_{f}}^{(j)}+{{\beta }_{i}}^{(j+1)}\right)\,{\ddot {\theta }}_{j}+$ ${\lambda }^{(j+1)}\,{\ddot {\theta }}_{j+1}-{{\tau }_{i}}^{\left(j\right)}{\ddot {u}}_{j-1}-$ $\left({{\gamma }_{i}}^{\left(j+1\right)}+{{\gamma }_{f}}^{\left(j\right)}\right)\,{\ddot {u}}_{j}-$

${-{{\tau }_{f}}^{\left(j+1\right)}{\ddot {u}}_{j+1}+}{{a}_{T}}^{(j)}{\theta }_{j-1}+$ $\left({{k}_{f,T}}^{(j)}+{{k}_{i,T}}^{(j+1)}\right)\,{\theta }_{j}+{{a}_{T}}^{(j+1)}\,{\theta }_{j+1}+$ ${{b}_{f,T}}^{(j)}{u}_{j-1}+$

$+\left({{b}_{i,T}}^{\left(j+1\right)}-{{b}_{f,T}}^{\left(j\right)}\right)\,{u}_{j}-$ ${{b}_{i,T}}^{\left(j+1\right)}{u}_{j+1}=-{{M}_{i}^{F}}^{\left(j+1\right)}-{{M}_{f}^{F}}^{(j)}$

(34)

a qual é válida para os nós e elementos ${\textstyle j=1,\,\ldots ,\,n-}$ $1$ . Para o nó n, a expressão de equilíbrio de momento de flexão fica:

${}{\lambda }^{(n)}{\ddot {\theta }}_{n-1}+\left({{\beta }_{f}}^{(n)}+{I}_{nac}^{fl}\right)\,{\ddot {\theta }}_{n}-$ ${{\tau }_{i}}^{\left(n\right)}{\ddot {u}}_{n-1}-{{\gamma }_{f}}^{\left(n\right)}\,{\ddot {u}}_{n}+$

$+\,{{a}_{T}}^{(n)}{\theta }_{n-1}+{{k}_{f,T}}^{(n)}\,{\theta }_{n}+{{b}_{f,T}}^{(n)}{u}_{n-1}-$ ${{b}_{f,T}}^{\left(n\right)}\,{u}_{n}=-{{M}_{f}^{F}}^{\left(n\right)}+{M}_{H}\left(t\right)$

(35)

na qual: ${\textstyle {M}_{H}\left(t\right)}$ é a função temporal do momento de flexão aplicado ao topo da torre; ${\textstyle {I}_{nac}^{fl}}$ é o momento de inércia da nacele em relação ao eixo que passa pelo diâmetro da seção transversal de topo da torre (eixo em torno do qual tem-se a rotação ${\textstyle {\theta }_{n}}$ ). Considerando-se que a nacele tem formato paralelepipédico e massa específica constante equivalente ${\textstyle {\rho }_{nac}}$ e que suas dimensões horizontal e vertical ao plano de análise são, respectivamente, ${\textstyle {a}_{nac}}$ (12 m) e ${\textstyle {h}_{nac}}$ (5 m), o momento de inércia é dado por:

(36)

sendo, ${\textstyle {M}_{nac}}$ a massa da nacele, do rotor e das pás (todos equipamentos locados no topo da torre), tomada igual a 200 t neste trabalho.

A expressão de equilíbrio de forças horizontais fica:

${}{{-\tau }_{f}}^{(j)}{\ddot {\theta }}_{j-1}-\left({{\gamma }_{i}}^{\left(j+1\right)}+\right.$ $\left.{{\gamma }_{f}}^{\left(j\right)}\right)\,{\ddot {\theta }}_{j}-{{\tau }_{i}}^{\left(j+1\right)}\,{\ddot {\theta }}_{j+1}+$ ${\sigma }^{\left(j\right)}{\ddot {u}}_{j-1}+\left({{\alpha }_{i}}^{\left(j+1\right)}+\right.$ $\left.{{\alpha }_{f}}^{\left(j\right)}\right)\,{\ddot {u}}_{j}+$

$+{\sigma }^{\left(j+1\right)}{\ddot {u}}_{j+1}{-}{{b}_{i,T}}^{\left(j\right)}{\theta }_{j-1}+$ $\left({{b}_{i,T}}^{\left(j+1\right)}-{{b}_{f,T}}^{\left(j\right)}\right)\,{\theta }_{j}+$ ${{b}_{f,T}}^{\left(j+1\right)}\,{\theta }_{j+1}-{{t}_{T}}^{\left(j\right)}{u}_{j-1}+$

$+\left({{t}_{T}}^{\left(j\right)}+{{t}_{T}}^{\left(j+1\right)}\right)\,{u}_{j}-$ ${{t}_{T}}^{\left(j+1\right)}{u}_{j+1}={{V}_{i}^{F}}^{\left(j+1\right)}+{{V}_{f}^{F}}^{(j)}$

(37)

também válida para os nós e elementos ${\textstyle j=1,\,\ldots ,\,n-}$ $1$ . Para o nó n, a expressão de equilíbrio de forças horizontais fica:

${}{{-\tau }_{f}}^{(n)}{\ddot {\theta }}_{n-1}-{{\gamma }_{f}}^{\left(n\right)}\,{\ddot {\theta }}_{n}+$ ${\sigma }^{\left(n\right)}{\ddot {u}}_{n-1}+\left({{\alpha }_{f}}^{\left(j\right)}+\right.$ $\left.{M}_{nac}\right)\,{\ddot {u}}_{n}-$

${-}{{b}_{i,T}}^{\left(n\right)}{\theta }_{n-1}-{{b}_{f,T}}^{\left(n\right)}\,{\theta }_{n}-$ ${{t}_{T}}^{\left(n\right)}{u}_{n-1}+{{t}_{T}}^{\left(n\right)}\,{u}_{n}={{V}_{f}^{F}}^{\left(n\right)}+$ ${F}_{H}(t)$

(38)

na qual, ${\textstyle {F}_{H}(t)\,}$ é a função temporal da força horizontal aplicada ao topo da torre.

A expressão de equilíbrio de forças verticais é dada por:

${}\,{\eta }^{(j)}{\ddot {v}}_{j-1}+\left({{\nu }_{f}}^{\left(j\right)}+{{\nu }_{i}}^{\left(j+1\right)}\right){\ddot {v}}_{j}+$ ${\eta }^{(j+1)}{\ddot {v}}_{j+1}{{-r}_{T}}^{\left(j\right)}{v}_{j-1}+\left({{r}_{T}}^{\left(j+1\right)}+\right.$ $\left.{{r}_{T}}^{\left(j\right)}\right)\,{v}_{j}-$

$-{{r}_{T}}^{\left(j+1\right)}{v}_{j+1}=-{{N}_{i}^{F}}^{\left(j+1\right)}-{{N}_{f}^{F}}^{(j)}$

(39)

a qual é válida para os nós e elementos ${\textstyle j=1,\,\ldots ,\,n-}$ $1$ . Para o nó n, a expressão de equilíbrio de forças verticais fica:

${}\,{\eta }^{\left(n\right)}{\ddot {v}}_{n-1}+\left({{\nu }_{f}}^{\left(n\right)}+\right.$ $\left.{M}_{nac}\right){\ddot {v}}_{n}{{-r}_{T}}^{\left(n\right)}{v}_{n-1}+{{r}_{T}}^{\left(n\right)}\,{v}_{n}=$ $-{{N}_{f}^{F}}^{\left(n\right)}+P(t)$

(40)

na qual, ${\textstyle P(t)\,}$ é a função temporal da força vertical aplicada ao topo da torre.

A expressão de equilíbrio de momentos de torção é dada por:

${}{\delta }^{\left(j\right)}{\ddot {\varphi }}_{j-1}+\left({{\mu }_{f}}^{\left(j\right)}+\right.$ $\left.{{\mu }_{i}}^{\left(j+1\right)}\right)\,{\ddot {\varphi }}_{j}+{\delta }^{\left(j+1\right)}{\ddot {\varphi }}_{j+1}{-{\hat {s}}}^{\left(j\right)}{\varphi }_{j-1}+$ $\left({\hat {s}}^{\left(j+1\right)}+{\hat {s}}^{\left(j\right)}\right)\,{\varphi }_{j}-$

$-{\hat {s}}^{\left(j+1\right)}{\varphi }_{j+1}=-{{T}_{i}^{F}}^{\left(j+1\right)}-$ ${{T}_{f}^{F}}^{(j)}$

(41)

também válida para os nós e elementos ${\textstyle j=1,\,\ldots ,\,n-}$ $1$ . Para o nó n, a expressão de equilíbrio de momentos de torção fica:

${}{\delta }^{\left(n\right)}{\ddot {\varphi }}_{n-1}+\left({{\mu }_{f}}^{\left(n\right)}+\right.$ $\left.{I}_{nac}^{tr}\right)\,{\ddot {\varphi }}_{n}{-{\hat {s}}}^{\left(n\right)}{\varphi }_{n-1}+$ ${\hat {s}}^{\left(n\right)}\,{\varphi }_{n}=-{{T}_{f}^{F}}^{\left(n\right)}+$ $T(t)\,$

(42)

na qual, ${\textstyle T(t)\,}$ é a função temporal do momento de torção aplicado ao topo da torre; ${\textstyle {I}_{nac}^{tr}}$ é o momento de inércia da nacele em relação ao eixo vertical da torre (eixo em torno do qual têm-se as rotações ${\textstyle {\varphi }_{j}}$ ) dado por:

(43)

em que, ${\textstyle {b}_{nac}}$ é a largura da nacele (dimensão perpendicular ao plano de análise) igual a 4 m.

As Eqs.(34), (35), (37), (38), (39), (40), (41) e (42) podem ser concatenadas para a forma matricial no SGC da torre, conforme a Eq.(6).

Em certas análises dinâmicas, não é conveniente considerar todos os graus de liberdade possíveis da estrutura. Neste estudo, por exemplo, para análise dinâmica da torre, se pretende estudar apenas as vibrações laterais. Portanto, para este caso específico, os graus de liberdade de interesse são aqueles que definem o movimento lateral da torre (designados graus de liberdade dinâmicos), podendo ser omitidos da análise os demais graus de liberdade, sem, contudo, deixar de se considerar a influência destes no comportamento global da estrutura. Assim, particionando-se a equação matricial de equilíbrio dinâmico, de maneira a separar os deslocamentos transversais à torre ( ${\textstyle {u}_{j}}$ ) dos demais, tem-se:

\left[{\begin{matrix}\left[{\boldsymbol {M}}_{\Delta \Delta }\right]&\left[{\boldsymbol {M}}_{\Delta u}\right]\\\left[{\boldsymbol {M}}_{u\Delta }\right]&\left[{\boldsymbol {M}}_{uu}\right]\end{matrix}}\right]\,\left\{{\begin{matrix}\left\{{\ddot {\boldsymbol {\Delta }}}\right\}\\\left\{{\ddot {\boldsymbol {u}}}\right\}\end{matrix}}\right\}+

(44)

em que: ${\textstyle \left\{{\boldsymbol {\Delta }}\right\}}$ é o vetor de deslocamentos verticais, de rotações torcionais e de rotações flexionais da torre (ou seja, os graus de liberdade condensados) e ${\textstyle \left\{{\boldsymbol {u}}\right\}}$ é o vetor de deslocamentos horizontais.

Considerando-se que as forças e os momentos inerciais e de amortecimento associados aos deslocamentos axiais e às rotações torcionais e flexionais sejam desprezíveis, a Eq.(44) pode ser escrita da seguinte forma:

(45)

de onde:

$\left[{\boldsymbol {M}}_{uu}\right]\left\{{\ddot {\boldsymbol {u}}}\right\}+$ $\left[{\boldsymbol {K}}_{u\Delta }\right]\left\{{\boldsymbol {\Delta }}\right\}+$ $\left[{\boldsymbol {K}}_{uu}\right]\left\{{\boldsymbol {u}}\right\}=$ $\left\{{\boldsymbol {F}}_{u}\right\}$

(46)

isolando-se o vetor de deslocamentos ${\textstyle \left\{{\boldsymbol {\Delta }}\right\}}$ da primeira:

(47)

e substituindo-se esta na segunda das Eqs.(46), obtém-se:

(48)

a partir da qual define-se a matriz de rigidez global condensada ${\textstyle \left[{\overline {\boldsymbol {K}}}_{u}\right]}$ [15, 16]:

(49)

e o vetor de forças transversais à torre condensado ${\textstyle \left\{{\overline {\boldsymbol {F}}}_{u}\right\}}$ :

(50)

Logo, a equação de equilíbrio dinâmico condensada fica:

(51)

Por fim, incorporam-se as forças de amortecimento, utilizando-se as matrizes de massa ${\textstyle \left[{\boldsymbol {M}}_{uu}\right]}$ e de rigidez ${\textstyle \left[{\overline {\boldsymbol {K}}}_{u}\right]}$ e o conceito de amortecimento proporcional, tem-se:

(52)

A modelagem no software ANSYS para análise estrutural e para o projeto da torre tubular de aço foi elaborada mediante o método dos elementos finitos (MEF), considerando materiais de comportamento elástico linear, do ponto de vista físico, e não linear, do ponto de vista geométrico. Para este fim, foi criado um modelo de elementos finitos no software ANSYS r.14.5 [20], no qual se considerou a torre engastada na base (Figura 2a) com 7272 elementos de casca, designado por SHELL 181 (Figura 3a), com 4 nós e 6 graus de liberdade por nó. A nacele foi modelada com elementos finitos sólidos tetraédricos, designados por SOLID187 (Figura 3b), com 10 nós e 3 graus de liberdade de translação por nó; de forma considerar que a nacele é formada por uma massa uniforme. O motivo que levou à utilização de um modelo em elementos finitos detalhado e outro em elementos de barra, simplificado, portanto, foi a necessidade de avaliar a confiabilidade e a precisão dos resultados numéricos obtidos.

O modelo com elementos finitos de casca foi complementado simulando-se a torre em conjunto com sua fundação (Figura 2b). Para tal, a sapata foi modelada com 11766 elementos sólidos tetraédricos, designados por SOLID 186 (Figura 3c), com 20 nós e 3 graus de liberdade de translação por nó. Além disto, com o objetivo de avaliar a interação solo-estrutura, a reação elástica do solo foi modelada com 2145 elementos de mola com rigidez axial, colocados na base da sapata e designados por COMBIN 14 (Figura 3d). A rigidez destes elementos foi avaliada a partir do valor médio do coeficiente de reação vertical, de uma areia com densidade relativa média, proposto por [21]. Assim, o valor do coeficiente reação vertical do solo, que é igual a 45023 kN/m³, foi multiplicado pela área de influência de cada nó da base da sapata que está em contato com o terreno.


(a) Engastado na base: EF de casca	(b) Com base flexível: EF de casca, sólido e de rigidez axial
Figura 2. Configuração dos modelos analisados


(a) Elemento SHELL 181	(b) Elemento SOLID 187

(c) Elemento SOLID 186	(d) Elemento COMBIN 14
Figura 3. Representação geométrica dos elementos finitos utilizados [20]

5. Resultados e discussão

5.1 Análise modal

Para obter as matrizes de massa, amortecimento e rigidez do modelo de torre é necessário estabelecer o número de elementos finitos a ser discretizada a torre. Devido à extensão das matrizes, só serão apresentadas abaixo as matrizes condensadas aos graus de liberdade dinâmicos (deslocamentos transversais à torre) com n = 8 (número de EF da torre), conforme Eqs.(53), (54) e (55):

(53)

em unidades do SI.

(54)

em unidade do SI.

(55)

em unidades do SI. Para montagem da matriz de amortecimento foi considerado uma razão de amortecimento ${\textstyle \xi }$ igual a 0,008 (conforme indicações de [22,23]) para o modo fundamental e para o quinto modo de vibração, que é considerado como o modo mais alto que contribui significativamente para a resposta da torre. Assim, as constantes ${\textstyle {a}_{0}}$ e ${\textstyle {a}_{1}}$ do amortecimento de Rayleigh resultam iguais a 2,940 x 10^-2 s^-1 e 1,441 x 10^-4 s, respectivamente.

Partindo das matrizes de massa ${\textstyle \left[{\boldsymbol {M}}_{E}\right]}$ e de rigidez ${\textstyle \left[{\boldsymbol {K}}_{E}\right]}$ com n = 16 e resolvendo o problema de autovalores e de autovetores, inerente à questão da vibração livre, obtêm-se as frequências naturais e os respectivos modos de vibração não amortecidos da torre (Tabela 2 e Figura 4).

Tabela 2. Dados de vibração livre da torre.

Ordem	1ª	2ª	3ª	4ª	5ª	6ª
Natureza	Flexão	Flexão	Torção	Flexão	Flexão	Torção
Frequência (Hz)	0,29745	1,72847	3,98035	4,52973	8,05650	8,75028
Período (s)	3,36187	0,57855	0,25123	0,22076	0,12412	0,11428
Frequência angular (rad/s)	1,86896	10,86030	25,00930	28,46116	50,62049	54,97963

Figura 4. Deslocamentos transversais à torre dos cinco primeiros modos de vibração

Neste caso, a frequência fundamental da torre (0,29745 Hz) é a mais preocupante quanto à possibilidade de haver ressonância tanto para o caso de vibração na direção do vento (along-wind) quanto para a vibração perpendicular a direção do vento (across wind) provocada pelo fenômeno de desprendimento cadenciado de vórtices (vortex shedding). Adicionalmente, comenta-se que a frequência fundamental é extremamente baixa resultando em um período de vibração de 3,362 s, tempo necessário a uma oscilação completa, o que corrobora a com a flexibilidade desta estrutura. Além disso, é mais provável que haja a excitação do modo fundamental, pois, por exemplo, assim como a velocidade do vento cresce com a altura em relação ao nível do solo, os deslocamentos horizontais à torre também aumentam com relação à altura no 1º modo de vibração e também se tem uma massa concentrada no topo (nacele), da mesma ordem de grandeza da massa da torre.

O primeiro modelo para análise modal via ANSYS foi discretizado com EF de casca, no qual a torre foi considerada engastada na base. O método de extração dos modos de vibração foi o Block Lanczos [24], a partir do qual foram extraídos os 18 primeiros modos e frequências de vibração da torre. Foram extraídos os 18 primeiros modos de vibração para se detectar os 6 modos de vibração obtidos no modelo discretizado com EF de barra, previamente apresentado.

A primeira e a segunda frequências são referentes aos primeiros modos de flexão nos planos XY e YZ, respectivamente (o eixo Y está posto na vertical no ANSYS). Há uma pequena diferença entre estas frequências devido à diferente disposição geométrica da nacele em relação aos planos XY e YZ. A terceira e a quarta frequências referem-se a outros dois modos de flexão nos planos XY e YZ, respectivamente. Os modos de vibração 5 e 6 são modos de oscilação ovais (modos de ovalização) com mesma frequência, que não são captados na análise com EF de barra. Estes tipos de modos surgem devido à reduzida espessura da parede da torre, entretanto, podem ser eliminados ou posicionados em frequências superiores (de menor importância para a resposta da estrutura) com a utilização de enrijecedores transversais à torre (anéis de rigidez). A sétima frequência é referente ao primeiro modo de torção em torno do eixo Y. Os modos 8 e 9 também são modos de flexão nos planos XY e YZ, respectivamente. Tornam a aparecer modos de ovalização nas frequências 10, 11, 12, 13, 14 e 15; estes são modos de vibrações locais acoplados, pois enquanto em um modo de vibração de flexão ou de torção globais tem-se deslocamento segundo um determinado eixo, nos modos de oscilação ovais há deslocamentos em mais de um eixo coordenado. Destes modos de ovalização, os seguintes pares têm a mesma frequência: 10º e 11º; 12º e 13º; e, 14º e 15º. Finalmente, as frequências 16 e 18 correspondem a outros modos de flexão nos planos XY e YZ, respectivamente; e a frequência 17 a mais um modo de torção em torno do eixo Y.

Na Figura 5 têm-se as representações gráficas dos modos de vibração descritos acima. Para os modos de vibração flexionais e torcionais, em que se tem deslocamento/rotação preponderante segundo um determinado eixo, os deslocamentos modais foram plotados segundo o respectivo eixo coordenado, a saber: translação em X para o 1º, 3º, 8º e 16º modos; translação em Z para o 2º, 4º, 9º e 18º modos; e, rotação em torno de Y para o 7º e 17º modos. Adicionalmente, para os modos de ovalização (5º, 6º, 10º, 11º, 12º, 13º, 14º e 15º), plotaram-se o módulo vetorial dos deslocamentos translacionais, uma vez que, nestes modos, têm-se deslocamentos pronunciados em todas as direções coordenadas.


Plano XY.	Plano YZ.	Legenda (m).	Plano XY.	Plano YZ.	Legenda (m).
(a) Modo 1 (1º modo de flexão XY – 0,296100 Hz).			(b) Modo 2 (1º modo de flexão YZ – 0,296357 Hz).

Plano XY.	Plano YZ.	Legenda (m).	Plano XY.	Plano YZ.	Legenda (m).
(c) Modo 3 (2º modo de flexão XY – 1,67352 Hz).			(d) Modo 4 (2º modo de flexão YZ – 1,69196 Hz).

Perspectiva.	Vista superior.	Legenda (m).	Perspectiva.	Vista superior.	Legenda (m).
(e) Modo 5 (1º modo de ovalização – 3,91390 Hz).			(f) Modo 6 (2º modo de ovalização – 3,91390 Hz).

Perspectiva.		Vista superior.		Legenda (rad).
(g) Modo 7 (1º modo de Torção Y – 4,03033 Hz).

Plano XY.	Plano YZ.	Legenda (m).	Plano XY.	Plano YZ.	Legenda (m).
(h) Modo 8 (3º modo de flexão XY – 4,46417 Hz).			(i) Modo 9 (3º modo de flexão YZ – 4,61078 Hz).

Perspectiva.		Vista superior.		Legenda (m).
(j) Modo 10 (3º modo de ovalização – 4,72284 Hz).

Perspectiva.		Vista superior.		Legenda (m).
(k) Modo 11 (4º modo de ovalização – 4,72284 Hz).

Perspectiva.		Vista superior.		Legenda (m).
(l) Modo 12 (5º modo de ovalização – 5,61856 Hz).

Perspectiva.		Vista superior.		Legenda (m).
(m) Modo 13 (6º modo de ovalização – 5,61866 Hz).

Perspectiva.		Vista superior.		Legenda (m).
(n) Modo 14 (7º modo de ovalização – 7,27205 Hz).

Perspectiva.		Vista superior.		Legenda (m).
(o) Modo 15 (8º modo de ovalização 7,27229 Hz).

Plano XY.	Plano YZ.	Legenda (m).	Plano XY.	Plano YZ.	Legenda (m).
(p) Modo 16 (4º modo de flexão XY – 8,35388 Hz).			(q) Modo 18 (4º modo de flexão YZ – 8,93076 Hz).

Perspectiva.		Vista superior.		Legenda (rad).
(r) Modo 17 (2º modo de torção Y – 8,87291 Hz).
Figura 5. Modos de vibração da torre engastada na base e modelada com EF de casca

No segundo modelo para análise modal via ANSYS, a torre foi discretizada com EF de casca e a fundação foi modelada com EF sólidos tetraédricos. Para este modelo do sistema da torre, foram obtidos os 50 primeiros modos de vibração. Foram extraídos os 50 primeiros modos de vibração para se detectar o 2ª modo de vibração torcional, uma vez que, este passou para uma frequência mais alta, quando comparado com vários outros modos de ovalização e flexionais que tinham frequências superiores a 2ª frequência torcional no modelo engastado na base. O nível de flexibilidade da fundação foi considerado de maneira semelhante ao feito na análise de estabilidade da torre exposta no capítulo 3; entretanto, o ponto de equilíbrio do conjunto torre-fundação apoiado sobre o solo deformável foi estabelecido utilizando-se a velocidade média de vento ( ${\textstyle {\overline {V}}}$ , definida em [22]), de maneira que, nesta situação, 10% área da base da sapata está sem contato com o solo (220 dos 2145 nós da base da sapata).

Assim como no modelo engastado na base, a primeira e a terceira e a segunda e a quarta frequências são referentes ao 1º e ao 2º modos de flexão no plano XY e YZ, respectivamente. Além disso, também há uma pequena diferença entre as 1ª e 2ª e entre a 3ª e 4ª frequências devido à diferente disposição geométrica da nacele em relação aos planos XY e YZ. Os pares de modos de vibração com a mesma frequência 5 e 6, 7 e 8 e 9 e 10 são modos de oscilação ovais (modos de ovalização). A 11ª frequência é referente ao 1º modo de torção em torno do eixo Y. Percebe-se, em comparação ao modelo com base engastada, que quatro modos de ovalização (3º, 4º, 5º e 6º modos de ovalização) passam a aparecer antes do 1º modo torcional da torre e que o 5º par de modos de ovalização (9º e 10º modos de ovalização) encontra-se entre as 18 primeiras frequências da estrutura com base flexível. Os modos 12 e 13 são os 3º’s modos de flexão nos planos XY e YZ, respectivamente. Tornam a aparecer pares de modos de ovalização com mesmas frequências: 14 e 15; 17 e 18. Finalmente, as frequências 16 e 19 correspondem a outros modos de flexão nos planos XY e YZ, respectivamente; e apenas a frequência 36 corresponde ao 2º modo de torção em torno do eixo Y. Entre a frequência 20 e 35 tem-se os seguintes modos, em ordem: 6º par de modos de ovalização (20º e 21º modos); 7º par de modos de ovalização (22º e 23º modos); 8º par de modos de ovalização (24º e 25º modos); 9º par de modos de ovalização (26º e 27º modos); 5º modo de flexão XY (28º modo); 10º par de modos de ovalização (29º e 30º modos); 11º par de modos de ovalização (31º e 32º modos); 12º par de modos de ovalização (33º e 34º modos); 5º modo de flexão YZ (35º modo).

Na Figura 6 são apresentadas as representações gráficas de alguns dos modos de vibração descritos acima. Para os modos de vibração flexionais e torcionais, em que se tem deslocamento preponderante segundo um determinado eixo, plotaram-se os deslocamentos modais segundo o respectivo eixo coordenado, a saber: translação em X para o 1º, 3º, 12º e 16º modos; translação em Z para o 2º, 4º, 13º e 19º modos; e, rotação em torno de Y para o 11º e 36º modos. Adicionalmente, para os modos de ovalização (5º, 6º, 7º, 8º, 9º, 10º, 14º, 15º, 17º e 18º), plotaram-se o módulo vetorial dos deslocamentos translacionais, uma vez que, nestes modos, têm-se deslocamentos pronunciados em todas as direções coordenadas.


Plano XY.	Plano YZ.	Legenda (m).	Plano XY.	Plano YZ.	Legenda (m).
(a) Modo 1 (1º modo de flexão XY – 0,282205 Hz).			(b) Modo 2 (1º modo de flexão YZ – 0,283639 Hz).

Plano XY.	Plano YZ.	Legenda (m).	Plano XY.	Plano YZ.	Legenda (m).
(c) Modo 3 (2º modo de flexão XY – 1,58747 Hz).			(d) Modo 4 (2º modo de flexão YZ – 1,61441 Hz).

Perspectiva.	Vista superior.	Legenda (m).	Perspectiva.	Vista superior.	Legenda (m).
(e) Modo 5 (1º modo de ovalização – 2,34283 Hz).			(f) Modo 6 (2º modo de ovalização – 2,34301 Hz).

Perspectiva.		Vista superior.		Legenda (rad).
(g) Modo 7 (3º modo de ovalização – 2,75494 Hz).

Perspectiva.		Vista superior.		Legenda (rad).
(h) Modo 8 (4º modo de ovalização – 2,75516 Hz).

Perspectiva.		Vista superior.		Legenda (rad).
(i) Modo 9 (5º modo de ovalização – 3,96102 Hz).

Perspectiva.		Vista superior.		Legenda (rad).
(j) Modo 10 (6º modo de ovalização – 3,96172 Hz).

Perspectiva.		Vista superior.		Legenda (rad).
(k) Modo 11 (1º modo de torção Y – 4,00274 Hz).

Plano XY.	Plano YZ.	Legenda (m).	Plano XY.	Plano YZ.	Legenda (m).
(l) Modo 12 (3º modo de flexão XY – 4,25413 Hz).			(m) Modo 13 (3º modo de flexão YZ – 4,41612 Hz).

Perspectiva.		Vista superior.		Legenda (m).
(n) Modo 14 (7º modo de ovalização – 5,96165 Hz).

Perspectiva.		Vista superior.		Legenda (m).
(o) Modo 15 (8º modo de ovalização – 5,96201 Hz).

Plano XY.	Plano YZ.	Legenda (m).	Plano XY.	Plano YZ.	Legenda (m).
(p) Modo 16 (4º modo de flexão XY – 7,93075 Hz).			(q) Modo 19 (4º modo de flexão YZ – 8,54675 Hz).

Perspectiva.		Vista superior.		Legenda (rad).
(r) Modo 17 (9º modo de ovalização – 8,30957 Hz).

Perspectiva.		Vista superior.		Legenda (rad).
(s) Modo 18 (10º modo de ovalização – 8,31100 Hz).

Perspectiva.		Vista superior.		Legenda (rad).
(t) Modo 37 (2º modo de torção Y – 13,7310 Hz).
Figura 6. Modos de vibração da torre engastada na base e modelada com EF de casca

Percebe-se a influência da flexibilidade da base da torre, no sentido de diminuir as frequências, em praticamente todos os modos de vibração; tendo como exceção a 2ª frequência torcional em torno do eixo Y. Além disso, a adição da sapata contribui com um aumento considerável da massa do modelo, o que naturalmente diminui o valor das frequências; embora a sapata esteja situada nas cotas mais baixas, fato que impede uma diminuição acentuada de tais frequências.

Outra questão observada no modelo com base flexível é notada no 8º modo, por exemplo, no qual há rotações e deslocamentos diferenciados em relação ao respectivo modo (3º modo de flexão YZ) do modelo engastado na base, tanto ao longo do tubo da torre quanto na nacele, que não mais rotaciona apenas em torno do eixo X. Uma característica a mais, relacionada ao modelo com base flexível, é a inexistência de simetria dos modos de ovalização em relação aos planos coordenados XY e YZ, quando comparados com os respectivos modos de ovalização do modelo engastado na base.

Na Figura 7, na qual estão plotados os deslocamentos verticais da sapata da torre (direção Y), observa-se a rotação da base da torre devido à flexibilidade da fundação.


(a) Sapata do 1º modo de vibração

(b) Sapata do 3º modo de vibração
Figura 7. Deslocamentos verticais da sapata (m)

Por fim, expõe-se na Tabela 3 um resumo das frequências de vibração obtidas nos três modelos estudados (EF de barra; EF de casca com base engastada; e, EF de casca com base flexível).

Tabela 3. Frequências de vibração dos modelos estudados.

Natureza dos Modos	Eixos	Frequências (Hz)
Natureza dos Modos	Eixos	EF de barra	EF de casca com base engastada	EF de casca com base flexível
1º Flexão	XY	0,29745	0,296100	0,282205
1º Flexão	YZ	0,29745	0,296357	0,283639
2º Flexão	XY	1,72847	1,67352	1,58747
2º Flexão	YZ	1,72847	1,69196	1,61441
1º Torção	Y	3,98035	4,03033	4,00274
3º Flexão	XY	4,52973	4,46417	4,25413
3º Flexão	YZ	4,52973	4,61078	4,41612
4º Flexão	XY	8,05650	8,35388	7,93075
4º Flexão	YZ	8,05650	8,93076	8,54675
2º Torção	Y	8,75028	8,87291	13,7310

5.2. Análise transiente

As amplitudes dos termos do vetor de forças aplicadas horizontalmente à torre foram calculadas por meio dos parâmetros estabelecidos em [22], a saber: parcela flutuante de vento estabelecida segundo a diferença entre a velocidade de vento estática (velocidade básica do vento medida em um intervalo de tempo de 3 s) e a velocidade média do vento (medida em um intervalo de tempo de 10 min), conforme Eq.(56); e uma força horizontal proveniente do sistema nacele-rotor aplicada ao topo da torre [25]. Esta parcela flutuante e a força aplicada ao topo foram adotadas segundo um carregamento harmônico considerado ressonante com o modo fundamental de vibração da torre. Desta forma, considerou-se que a torre vibra em torno da posição de equilíbrio que é determinada pelo carregamento obtido a partir da velocidade média de vento. Logo:

(56)

na qual: ${\textstyle {q}_{w}\left(x\right)}$ é a função em relação ao nível geral do terreno da pressão dinâmica da parcela flutuante do vento em N/m²; ${\textstyle {\rho }_{ar}}$ é a massa específica do ar em condições ambientais normais [26], tomada igual a 1,225 kg/m³; ${\textstyle {V}_{3s}(x)}$ é a função em relação ao nível geral do terreno da velocidade de vento medida em um intervalo de tempo de 3 s; ${\textstyle {V}_{10min}(x)}$ é a função em relação ao nível geral do terreno da velocidade de vento medida em um intervalo de tempo de 10 min.

Assim, na Eq.(57) abaixo, se apresenta o vetor condensado das forças referentes à parcela flutuante do vento consideradas harmônicas ressonantes, com n = 8:

\left\{{\overline {\boldsymbol {F}}}_{u}\right\}=F_{u}=\left({\begin{matrix}4.374\times 10^{3}\\5.95\times 10^{3}\\6.417\times 10^{3}\\6.696\times 10^{3}\\6.797\times 10^{3}\\6.544\times 10^{3}\\7.251\times 10^{3}\\1.137\times 10^{5}\end{matrix}}\right)sen\,{\boldsymbol {(}}{\omega }_{1}t{\boldsymbol {)}}

(57)

em unidades do SI. Percebe-se que a amplitude da força aplicada ao topo da torre tem ordem de grandeza 100 vezes maior que as demais amplitudes, o que corrobora com a excitação do modo fundamental da torre. Comenta-se ainda, que não foi considerado o efeito aleatório da ação de vento na determinação do vetor de forças acima (Eq.(57)), mas apenas a intensidade da ação de vento que corresponde à sua parcela flutuante.

Considerando-se nulos os vetores de condições inicias (deslocamentos e velocidades), o gráfico de deslocamentos do topo da torre, quando submetida ao vetor de forças ressonantes da Eq.(57), é mostrado na Figura 8.

${\textstyle {u}_{n}}$ (m)
	$t\,(s)$
Figura 8. Gráfico de deslocamento no topo da torre ${\textstyle {u}_{n}}$ versus tempo

Percebe-se na Figura 8, que o deslocamento cresce quase que infinitamente com o tempo, tendo valor máximo aproximadamente igual a 6 m, pois se trata de uma vibração forçada amortecida excitada por uma carga ressonante. Tal fato leva a torre a níveis de deslocamentos inaceitáveis do ponto de vista de projeto e, além disso, para um deslocamento no topo da ordem de 2 m, ter-se-ia o início do comportamento de não linearidade física na torre (entretanto, o comportamento implementado considerou o material em regime elástico linear, ou seja, não foi feita análise física não linear). Entretanto, é importante mencionar que a carga aplicada além de ser ressonante com o 1º modo da torre, está atuando em todo o tempo de análise, o que leva, inevitavelmente, a estrutura a deslocamentos elevados.

Desta forma, resolveu-se realizar uma segunda análise que consiste no fato de avaliar os deslocamentos quando a excitação atua em parte do tempo de análise. Além disso, é considerado que durante a atuação do carregamento a torre vibra em torno da posição de equilíbrio estático definida pela atuação da força eólica média (estabelecida pela velocidade média de vento ${\textstyle {\overline {V}}}$ ), ou seja, a parcela flutuante do vento, adotada como harmônica, incide em torno da sua parcela média. Na Figura 9 é apresentada a função dos deslocamentos do topo da torre (em torno do deslocamento estático) quando esta é submetida ao vetor de forças da Eq.(57) durante 5 s (semelhante ao efeito de uma rajada). Após o tempo de atuação da força ( ${\textstyle {t}_{F}}$ ), a torre passa a movimentar-se sob vibração livre amortecida em torno da posição de deslocamento nulo; entretanto, mesmo após 100 s de análise, a torre continua a vibrar com amplitude relativamente elevada (15 cm), acarretando problemas de fadiga do material.

${\textstyle {{u}^{FL}}_{n}}$ (m)
	$t\,(s)$
Figura 9. Gráfico de deslocamento no topo da torre ${\textstyle {{u}^{FL}}_{n}}$ versus tempo

Será investigado nesta seção o efeito dinâmico de desprendimento cadenciado de vórtices (vórtices de von Kárman), que causa movimentos transversais à direção do vento (across-wind) em certa frequência de desprendimento de um par de vórtices, potencialmente preocupante quando esta se iguala a uma das frequências naturais da estrutura, dentro da faixa de velocidade de vento esperada. Mais uma vez, assim como no caso da vibração na direção do vento (along-wind), tem-se especial atenção à possibilidade de excitação do 1º modo de vibração da torre devido ao perfil de velocidade de vento, que obedece a uma função crescente em relação ao nível geral do terreno.

Dependendo do número de Reynolds do fluxo, Eq.(58), referido ao diâmetro externo da torre (dimensão característica da torre), observam-se três regiões características [23]: subcrítica (300 < Re < 1,5 10⁵), transicional (1,5 10⁵ < Re < 3,5 10⁶) e supercrítica (3,5 10⁵ < Re).

(58)

na qual: ${\textstyle Re}$ é o número adimensional de Reynolds; ${\textstyle {V}_{w}}$ é a velocidade do fluxo de vento; ${\textstyle D}$ é o diâmetro externo do tubo da torre; e, ${\textstyle \nu }$ é a viscosidade cinemática do fluido, que, neste caso, é igual a 70000^-1 m²/s para o ar [22]. Para o perfil de velocidades média do vento ${\textstyle {\overline {V}}}$ em relação ao nível geral do terreno, o vetor de número de Reynolds para torre com n = 8, fica:

(59)

A velocidade de vento, para a qual a frequência de desprendimento de um par de vórtices coincide com uma frequência natural (neste caso, será a 1ª frequência de vibração da torre), é chamada de velocidade crítica ${\textstyle {V}_{cr}}$ e é dada por:

(60)

em que, ${\textstyle St}$ é o número adimensional de Strouhal definido por:

(61)

na qual, ${\textstyle {f}_{s}}$ é a frequência de desprendimento de um par de vórtices. Segundo [22], St = 0,28 para Re > 10⁶. Logo, o vetor de velocidades críticas para a torre com n = 8 é dado por:

(62)

em unidades do SI. Este vetor apresenta valores distintos e decrescentes de velocidade ao longo da altura da torre, de forma que não é possível excitar o desprendimento de vórtices ressonante em toda a altura da torre. Segundo [26], o desprendimento de vórtices ocorre em trechos ao longo do comprimento da torre formando as chamadas células de vórtices, sendo constante a frequência de desprendimento em cada célula. Além disso, os efeitos dinâmicos do desprendimento cadenciado de vórtices são possíveis se a velocidade crítica for igual ou inferior à máxima velocidade média de vento, que, neste caso, para o topo da torre é igual a 39,7 m/s. Desta forma, a velocidade crítica é atingida primeiro no topo da torre à medida que o perfil de velocidades cresce desde zero na base até atingir o valor crítico no topo da torre.

Segundo [23], devido ao fato de o processo de desprendimento cadenciado de vórtices ser aproximadamente senoidal, é razoável modelar o fenômeno de desprendimento de vórtices por um carregamento transversal senoidal distribuído por unidade de comprimento de torre, ${\textstyle {q}_{L}}$ (em N/m), conforme:

(63)

na qual: ${\textstyle {C}_{L}}$ é coeficiente adimensional de sustentação tomado igual a 0,5, a favor da segurança, para a faixa de número de Reynolds em análise [23,27]; ${\textstyle {\omega }_{s}}$ é a frequência angular de desprendimento de um par de vórtices de von Kárman, em rad/s, que varia em relação ao nível geral do terreno, pois o carregamento utilizado refere-se a situação em que a velocidade de vento se iguala a velocidade crítica apenas no topo da torre (condição para possível excitação do modo de vibração fundamental). Na Eq.(64) tem-se o vetor de frequências angulares de desprendimento de um par de vórtices quando a torre é discretizada em EF de barra com n = 8:

(64)

em unidades do SI.

O vetor de amplitudes das forças transversais de sustentação ${\textstyle \left\{{\boldsymbol {F}}_{L}\right\}}$ , que têm o efeito do desprendimento de vórtices de von Kárman, com n = 8, é dado por:

(65)

em unidade do SI.

Com a aplicação do referido vetor de forças tem-se a resposta da torre. Assim, na Figura 10 é plotado o gráfico de deslocamentos do topo da torre quando sujeita ao vetor de cargas transversais que simulam a ação do desprendimento de vórtices de Kárman.

${\textstyle {{u}^{v}}_{n}}$ (m)
	$t\,(s)$
Figura 10. Gráfico de deslocamento no topo da torre ${\textstyle {{u}^{v}}_{n}}$ versus tempo

Percebe-se na Figura 10 que há um crescimento dos deslocamentos com tempo, mas isso ocorre até certo ponto (abaixo de 5 mm). Apesar das forças aplicada ao topo ser ressonante, as demais forças aplicadas ao longo do comprimento da torre não são, ou seja, há falta de sincronismo das forças excitadoras, o que impede que os deslocamentos cresçam infinitamente como é o caso quando tem-se todas as forças aplicadas ressonantes. Adicionalmente, a intensidade das forças transversais é baixa quando comparada com as forças postas na direção do vento do item anterior, de maneira a ter-se baixa intensidade de deslocamentos no caso estudado nesta seção.

6. Conclusões

Os resultados da análise modal feita, tanto para o modelo com elementos finitos de casca, com e sem base flexível, quanto para o modelo de elementos finitos de barra mostram-se similares, principalmente com relação ao 1º modo de vibração da torre. Portanto, o modelo representado por EF de barra pode ser utilizado como representativo do comportamento dinâmico da torre quando esta é submetida às excitações ressonantes à frequência fundamental da torre.

Com isso, a ação da parcela flutuante do vento e a excitação que simula o desprendimento cadenciado de vórtices de von Kárman foram modeladas por forças harmônicas ressonantes, considerando a tendência de que haja a excitação do modo fundamental da torre. A intensidade crescente das forças do vento com a altura em relação ao nível do solo, o aumento dos deslocamentos horizontais à torre com relação à altura no 1º modo de vibração e a massa concentrada no topo (nacele), da mesma ordem de grandeza da massa da torre, são fatores que contribuem para excitação do 1º modo de vibração da torre.

Adicionalmente, a influência da flexibilidade da fundação nos modos de vibração da torre foi investigada. Foram extraídos os 50 primeiros modos de vibração do modelo com base flexível para se detectar o 2ª modo de vibração torcional. Este passou para uma frequência mais alta, quando comparados aos outros modos locais e flexionais que tinham frequências superiores a 2ª frequência torcional do modelo engastado na base. Assim, com a consideração da flexibilidade da base, quase a totalidade dos valores das frequências de vibração diminuiu e alguns modos de vibração locais superiores passaram a ficar alocados em frequências mais baixas.

Conclui-se, então, que os resultados desta pesquisa envolvem contribuições de interesse prático imediato, uma vez que, se pretende desenvolver subsídios para análises estruturais e de controle de vibrações, além de projetos de torres e fundações para aerogeradores de multi-megawatt de potência nominal, a serem implantados no Brasil (especificamente, no estado de Pernambuco).

Referências

[1] Engström S., Lyrner T., Hassanzadeh M., Stalin T., Johansson J. Tall towers for large wind turbines. Estocolmo, Elforsk, 2010.

[2] Barifouse R., Schreiber M. Como o Nordeste virou principal polo da energia eólica no Brasil. BBC BRASIL, São Paulo e Brasília, 13 nov. 2015. Disponível em: [<http://www.bbc.com/portuguese/noticias/2015/11/151110_energia_eolica_nordeste_rb>. <http://www.bbc.com/portuguese/noticias/2015/11/151110_energia_eolica_nordeste_rb>.].

[3] Hansen, M.O.L., Sorensen J.N., Voutsinas S., Sorensen N., Madsen H.A. State of the art in wind turbine aerodynamics and aeroelasticity. Prog. in Aerosp. Sci., 42(4):285–330, 2006.

[4] Hau, E. Wind turbines: Fundamentals, technologies, application, economics. Munich (Germany), 2nd ed., Springer, 2006.

[5] Bazeos N., Hatzigeorgiou G.D., Hondros I.D., Karamaneas H., Karabalis D.L., Beskos D.E. Static, seismic and stability analyses of a prototype wind turbine steel tower. Eng. Struct., 24(8):1015-1025, 2002.

[6] Sirqueira A.S. Comportamento estrutural de torres de aço para suporte de turbinas eólicas. 112f. Dissertação (Mestrado em Engenharia Civil) – Faculdade de Tecnologia e Ciências, Universidade do Estado do Rio de Janeiro, Rio de Janeiro, 2008.

[7] Wang J., Qin D., Lim T.C. Dynamic analysis of horizontal axis wind turbine by thin-walled beam theory. J. of Sound and Vib., 329(17):3565-3586, 2010.

[8] Dellezzopolles Jr. C.F. Análise dinâmica de torres de energia eólica. 83 f. Dissertação (Mestrado em Estruturas e Construção Civil) – Faculdade de Tecnologia, Universidade de Brasília, Brasília, 2011.

[9] Adhikari S., Bhattacharya S. Dynamic analysis of wind turbine towers on flexible foundations. Shock. and Vib., 19:37-56, 2012.

[10] Oliveira L.F.M.P. Análise do comportamento dinâmico de torres de geradores eólicos. 96 f. Dissertação (Mestrado em Engenharia Civil) – Universidade do Porto, Faculdade de Engenharia da Universidade do Porto, Porto, 2012.

[11] Avila S.M., Shuzu, M.A.M., Pereira W.M., Santos L.S., Morais M.V.G., Prado Z.J.G. Numerical modeling of the dynamic behavior of a wind turbine tower. Proceeding of International Conference on Vibration Problems (ICOVP), 11, Lisboa, 2013.

[12] Cook R.D., Malkus D.S., Plesha M.E., Witt, R.J. Concepts and applications of finite element analysis. 4th ed. Madison, John Wiley & Sons, 2002.

[13] Huebner K.H., Dewhirst D.L., Smith D.E., Byrom T.G. The finite element method for engineers. 4th ed. John Wiley & Sons, 2001.

[14] Humar J.L. Dynamics of structures. Ottawa (Canada), A.A. Balkema Publishers, 2nd ed., Lisse, 2002.

[15] Clough R.W., Penzien J. Dynamics of structures. Berkeley, Computers and Structures, 3rd ed., 2003.

[16] Chopra A.K. Dynamics of structures: Theory and applications to earthquake engineering. Practice Hall, 4th ed., 2012.

[17] European Committee For Standardization. EN 10025-2: Hot rolled products of structural steels – Part 2: Technical delivery conditions for non-alloy structural steels. Brussels, 2004.

[18] Siemens. Siemens D3 platform – 3.0-MW and 3.2 – MW direct drive wind turbines: Reduced complexity, increased profitability. Erlangen, Germany, 2014. Disponível em: <[https https]>. ://www.energy.siemens.com/br/pool/hq/power-generation/renewables/wind-power/platform20brochures/D3%20Onshore%20brochure_ENGLISH_Apr2014_WEB.pdf>.

[19] Lima D.M. Análise da estabilidade elástica, análise dinâmica e controle de vibração em torres tubulares de aço para aerogeradores de eixo horizontal. Tese (Doutorado em Engenharia Civil), Universidade Federal de Pernambuco, Recife, 2018.

[20] Swanson Analysis Systems Inc. ANSYS Mechanical User's Guide. Release 14.5. Canonsburg (Pennsylvania, USA): South pointe, 275 Technology Drive, PA 15317, 2012.

[21] Terzaghi K. Evaluation of coefficient of subgrade reaction. Geotechnique, 5(4):297-326, 1955.

[22] Associação Brasileira De Normas Técnicas – Norma Brasileira. NBR 6123: Forças devidas ao vento em edificações. Rio de Janeiro, 66 pp., 1988.

[23] Blevins, R.D. Flow-induced vibration. Malabar (Florida), Krieger Publishing Company, 2nd ed., 2001.

[24] Lanczos C. An efficient method for the solution of the eigenvalue problem of linear differential and integral operators. J. of Res. of Natl. Bur. of Stand., 45(4):255-282, 1950.

[25] Asibor A.I., Garcia J.R., Ramos M.C., Silva E.C.M., Araújo A.M. Wind turbine performance and loading calculations using aero elastic modelling. Proceeding of International Congress of Mechanical Engineering (ABCM), 23, Rio de Janeiro, 2015.

[26] Associação Brasileira De Normas Técnicas. NBR IEC 61400-1: Aerogeradores Parte 1: Requisitos de Projeto. Rio de Janeiro, 2008. 82p.

[27] Blessmann J. Introdução ao estudo das ações dinâmicas do vento. Porto Alegre, Editora da UFRGS, 2nd ed., 2005.