Análise da direção do vento no anemômetro de uma turbina eólica por meio da simulação em CFD christianoacb@al.insper.com.br danielcgc@insper.edu.br estersq@al.insper.edu.br felipeant@al.insper.edu.br Trabalho de Conclusão de Curso Relatório Versão Final do Projeto Final de Engenharia São Paulo - SP Junho 2022 Christiano Ayres de Carvalho Borges Daniel do Carmo Granja de Castro Ester dos Santos Quintino Felipe Aron Nudelman Tabacinik Análise da direção do vento no anemômetro de uma turbina eólica por meio da simulação em CFD Relatório (F) do Projeto Final de Engenharia Relatório apresentado ao curso de Engenharia, como requisito para o Trabalho de Conclusão de Curso. Professor Orientador: Prof. Dr. Paulo Roberto Bufacchi Mendes Mentores na Empresa: Augusto Cargnin Morcelli e Diego Schmitt Montero. Coordenador TCC/PFE: Prof. Dr. Luciano Pereira Soares São Paulo - SP Junho 2022 Christiano Ayres de Carvalho Borges Daniel do Carmo Granja de Castro Ester dos Santos Quintino Felipe Aron Nudelman Tabacinik Projeto Final de Engenharia apresentado ao programa Graduação em Engenharia (Computação/Mecânica/Mecatrônica) como requisito parcial para a obtenção do título de Bacharel em Engenharia. Orientador: Prof. Dr. Paulo Roberto Bufacchi Mendes Banca Examinadora ______________________________________ Paulo Roberto Bufacchi Mendes Insper ______________________________________ Frederico Augusto Alem Barbieri Insper ______________________________________ Vinicius Licks Insper Sumário RESUMO ............................................................................................................................................................... 4 1. INTRODUÇÃO ............................................................................................................................................ 7 1.1 ESCOPO DO PROJETO .................................................................................................................................. 8 1.2 RECURSOS.................................................................................................................................................. 8 1.3 CRONOGRAMA ........................................................................................................................................... 9 1.4 MAPEAMENTO DOS STAKEHOLDERS ......................................................................................................... 10 1.5 RISCOS ENVOLVIDOS ................................................................................................................................ 10 1.6 REVISÃO DO ESTADO DA ARTE ................................................................................................................. 12 2. METODOLOGIA ...................................................................................................................................... 18 2.1 CÁLCULO DA CAMADA LIMITE PARA GERAÇÃO DA MALHA PRÓXIMA ÀS SUPERFÍCIES SÓLIDAS ............... 18 2.2 DOMÍNIO .................................................................................................................................................. 20 2.3 CRIAÇÃO DA MALHA ................................................................................................................................ 22 2.4 CONDIÇÕES DE CONTORNO ...................................................................................................................... 23 2.5 PARÂMETROS PARA A SIMULAÇÃO .......................................................................................................... 23 2.6 SIMPLIFICAÇÕES DO MODELO CAD FORNECIDO PELA WEG .................................................................... 26 2.11 INDEPENDÊNCIA DE MALHA E BOAS PRÁTICAS APLICADAS. ..................................................................... 33 2.12 MODELOS DE TURBULÊNCIA .................................................................................................................... 46 2.13 ESCOAMENTO DO AR NA PÁ ..................................................................................................................... 46 2.14 SIMULAÇÕES ........................................................................................................................................... 48 2.15 DESIGN OF EXPERIMENTS (DOE) ............................................................................................................. 49 2.16 MATRIZ DE CORRELAÇÃO ........................................................................................................................ 53 3. RESULTADOS .......................................................................................................................................... 54 3.1 SIMPLIFICAÇÕES DO MODELO CAD ......................................................................................................... 54 3.2 SIMULAÇÕES ............................................................................................................................................ 63 3.3 VALIDAÇÃO DA MALHA APÓS A INDEPENDÊNCIA DE MALHA ................................................................... 65 3.2 CFL ......................................................................................................................................................... 73 3.3 MODELOS DE TURBULÊNCIA .................................................................................................................... 74 3.4 RESULTADO DAS SIMULAÇÕES DO DOE ................................................................................................... 77 3.5 EFEITO DA PASSAGEM DA PÁ EM FRENTE AO ANEMÔMETRO. ................................................................... 88 4. CONCLUSÃO E TRABALHOS FUTUROS .......................................................................................... 89 5. REFERÊNCIA BIBLIOGRÁFICA ......................................................................................................... 90 RESUMO Este projeto tem como objetivo analisar de forma numérica a direção do vento nos dois anemômetros da turbina eólica 4.2 da WEG para otimização do posicionamento da nacele (ângulo de yaw). O propósito é aumentar a captura de energia do vento pelas turbinas eólicas até a velocidade nominal pela correção da leitura feita pelos anemômetros. Para que isso ocorra, será necessário o uso do software STAR-CCM+ para realizar as simulações em Computational Fluid Dynamics (CFD). Palavras-chave: Turbina eólica; CFD; Otimização. 7 1. Introdução O projeto foi realizado juntamente com a WEG S.A, que inicialmente produzia motores elétricos, mas a partir da década de 80 ampliou suas atividades com a produção de componentes eletroeletrônicos, produtos para automação industrial, transformadores de força e distribuição, tintas líquidas e em pó, vernizes eletroisolantes e geradores [1]. Em 2013, soluções de turbinas eólicas foram introduzidas com sucesso e, a partir daí, foi realizado o projeto estratégico de um aerogerador, de 2,1 MW. Com a análise do seu desempenho, foi possível desenvolver o novo aerogerador com 4,2 MW de potência, que se mostrou mais adequado às condições de vento do Brasil e será o foco deste estudo. Para que um aerogerador, com suas partes evidenciadas na Figura 1, consiga gerar mais energia, ou seja, produzir a maior potência, é preciso que ele consiga extrair o máximo da energia do vento possível, isto ocorre quando as pás do rotor estão posicionadas perpendicularmente à direção do vento. Para fazer com que as pás acompanhem corretamente a direção do vento é necessária uma leitura, processamento das informações e controle de posição adequados. O objetivo do projeto é utilizar os anemômetros que estão localizados acima da nacele e atrás das pás do rotor para fazer a leitura da direção do vento e utilizá-la para a calibração do posicionamento das pás melhorando a eficiência da turbina eólica. Figura 1 - Vista lateral da turbina e descrição das partes. O vento que chega no anemômetro possui um direcionamento diferente do vento de entrada no rotor devido à turbulência gerada pela passagem do fluido pelas pás do rotor em movimento. Portanto, é válido questionar o porquê de o anemômetro não ser posicionado à frente do rotor já que assim ele leria a direção original do vento sem que houvesse nenhuma 8 variação. A resposta para esta pergunta é que para posicionar o anemômetro à frente da turbina seria necessário a instalação de um mastro com um tamanho próximo ao da torre e, fazer isso para cada uma delas é financeiramente inviável. Colocar o sensor fixado ao rotor faria com que ele também girasse e isso afetaria a leitura do sensor tal como posicionar atrás das pás. Portanto, é logico e financeiramente mais viável instalar o anemômetro atrás do rotor. 1.1 Escopo do projeto O projeto, engloba simulações numéricas que consideram diferentes direções horizontais e verticais do vento na entrada do rotor, velocidades do vento, rotações do rotor, densidades do ar e intensidades turbulentas que serão lidas pelo anemômetro. Portanto, análises da reação em cadeia e efeito esteira como geralmente ocorre em wind farms estão fora do escopo. As simulações feitas pelo grupo gerarão um modelo de calibração que poderá ser aplicado ao controlador das turbinas para melhorar a sua eficiência, sem necessidade de modificar a estrutura do equipamento. Entretanto, para que uma simulação seja utilizada dessa forma é necessário antes realizar a validação do modelo numérico com dados reais de algum experimento. A WEG poderia fornecer os dados experimentais de medição da direção do vento no anemômetro e antes da turbina, mas essa compilação levaria um tempo superior ao que o grupo poderia esperar. Ao buscar por outras instituições que pudessem fornecer os dados, o grupo entrou em contato com um engenheiro de sistemas de energia eólica distribuída da National Renewable Energy Lab (NREL) que informou não dispor dos dados experimentais necessários para avaliar a direção do vento nos anemômetros acima da nacele. Entretanto, informou que a Denmark Tecnology University (DTU) poderia ter os dados solicitados, mas nenhum dos contatos desta instituição retornou com as informações solicitadas. 1.2 Recursos O recurso mais importante do projeto é o tempo que os membros disponibilizaram ao projeto. Foi importante a colaboração de todos e a ajuda do orientador que possui muita experiência no ramo de CFD e conhecimento do STAR-CCM+, como também dos mentores da WEG. 9 Além disso, o grupo aproveitou ao máximo os recursos computacionais que o Insper possui, dado que as máquinas presentes nos laboratórios têm boa capacidade de processamento, o que auxiliou na agilidade das simulações. O projeto foi realizado com simulações computacionais no STAR-CCM+ e, como o Insper possui a licença do software, não houve nenhum custo para o projeto. 1.3 Cronograma Foi feito um cronograma, como pode ser observado na Figura 2, visando todos os trabalhos que foram realizados durante o projeto. Ele foi aprovado pela WEG e pelo orientador. Para o desenvolvimento do projeto, observa-se que as etapas que mais precisaram de dedicação de tempo foi a revisão bibliográfica e a simulação da turbina 4.2 da WEG. Figura 2 - Cronograma do PFE. A falta de disponibilidade de dados experimentais da NREL e da DTU, fez com que o grupo não realizasse essa etapa, como pode ser observado no cronograma e, utilizasse como base, o vídeo do canal "Applied Computational Fluid Dynamics" [3] que faz a simulação de uma turbina eólica de duas pás. As reuniões com a WEG foram acordadas ao término de cada reunião e via e-mail de acordo com a disponibilidade dos membros do grupo, coordenador e mentores da WEG. A etapa que se refere a simulação de turbina NREL/WEG não foi realizada por motivos de não haver dados para serem analisados até então. Houve a inserção do cronograma da análise dos dados experimentais enviados pela WEG em meados de Abril, como forma de validação. 10 As etapas que se referem ao relatório final e apresentações com o Insper e com a WEG, foram prolongadas de acordo com a data marcada para tais avaliações. 1.4 Mapeamento dos stakeholders Os principais stakeholders do projeto são os colaboradores da WEG, seguido dos alunos e orientador do Insper. É possível destacar a importância do projeto para a WEG como possível melhoria na eficiência da turbina, assim como para o Insper, na qual pode utilizar o projeto como fonte de estudo para próximos projetos. A matriz de responsabilidade de cada um dos principais contribuintes, foi evidenciada na Tabela 1. Tabela 1 - Matriz de responsabilidade. Função no Nome projeto Planejamento Gestão Escopo Execução Validação Christiano Ayres de Carvalho Borges Time de projeto R R R P A Daniel do Carmo Granja de Castro Time de projeto R R R P A Ester dos Santos Quintino Time de projeto R R R P A Felipe Aron Nudelman Tabacinik Time de projeto R R R P A Paulo Bufacchi Orientador C C C C C Mentor da Augusto empresa C C I I A Mentor da Diego empresa I I I I I P: participa da atividade I: deve ser informado C: deve ser consultado R: responsável pela atividade A: aprova 1.5 Riscos envolvidos Durante a discussão do projeto, foi feito a análise de risco para o seu sucesso. Para isso, foi subdividiu em categoria do risco, com medida qualitativa do impacto e da sua probabilidade de ocorrência, além da medida de redução através do plano de resposta. A partir disso, foram adotados 4 planos, sendo eles: • Eliminação – A partir da eliminação das causas que colocaria o projeto em risco. • Mitigação – A partir de ações para minimizar a probabilidade de ocorrer os problemas e instabilidades. 11 • Gerenciamento – A partir de medidas que viesse a impactar menos a qualidade e o desenvolvimento do projeto como um todo. • Empenho – A partir de força tarefa a fim de mitigar a ausência do membro naquele instante. É possível acompanhar os detalhes descritos acima na Tabela 2. Tabela 2 - Análise de risco. De acordo com a Tabela 2, no item 1, a falta de comunicação com a empresa impactaria altamente a realização do projeto, mas como desde o início houve uma boa comunicação com a WEG, através de reuniões e trocas de e-mail, foi colocado uma probabilidade muito baixa para esse item. E de fato, a comunicação ocorreu muito bem. No item 2, em relação as instabilidades com o STAR-CCM+, a situação do início do projeto era que o servidor de licença do Insper estava fora do ar e a VPN não estava disponível para acessar os computadores do Insper de fora da sua rede, sendo considerado como gravidade muito alta. O orientador conseguiu com a Siemens licenças provisórias para que fosse possível trabalhar com o STAR-CCM+ fora da faculdade. Cada arquivo de licença tem validade de cinco dias e foi utilizado 4 arquivos de licenças provisórias. O computador do Insper com mais processadores em meados de Abril ficou operacional para ser acessado de dentro da rede do Insper. Entretanto a VPN não foi possível de ser utilizada. Essa limitação, impactou o grupo de 12 obter e visualizar resultados em qualquer lugar. O grupo ficou limitado a ver os resultados das simulações que foram rodadas no Insper só quando estavam em campus, com a rede da faculdade conectada. No terceiro item, para que o escopo do projeto não sofresse grandes transformações ao longo do tempo foi preciso aplicar um gerenciamento a fim de mantê-lo como foi pré- estabelecido com a WEG. Para mitigar ainda mais o aumento do escopo, foi aplicado um planejamento de experimentos, para otimizar as simulações a serem desenvolvidas. No quarto item, em relação ao cenário global ainda há pandemia, e mesmo com uma diminuição do número de casos e aumento do número de pessoas vacinadas, ainda é provável que qualquer membro do grupo contraia a COVID-19, mas como o projeto apresenta a possibilidade de trabalho remoto este é um problema que não oferece um risco elevado. A solução será o isolamento social, medicação e descanso adequado e trabalhar o máximo possível remotamente até estar saudável novamente. Vale ressaltar que até o momento nenhum membro do grupo contraiu a doença. No quinto item, para que o grupo obtivesse uma boa curva de aprendizado, a solução foi realizar um trabalho em conjunto para compartilhar conhecimentos sobre o software STAR- CCM+ a fim de otimizar o seu manuseio e ter uma sensibilidade maior sobre causas e consequências dos problemas ou soluções que surgiram no desenvolvimento do projeto. Vale ressaltar que o trabalho não finaliza no plano de resposta. Uma etapa crucial foi o acompanhamento de cada risco, a documentação dos resultados para ser possível realizar análises do impacto causado, de maneira a focar na redução deste cada vez mais. 1.6 Revisão do estado da arte Pelo fato de o projeto ter como base a simulação computacional através do CFD, foi preciso buscar referências para aplicação dos parâmetros no software STAR-CCM+. Dentre eles pode-se citar, o tamanho do domínio, as características da malha, as condições de contorno, os modelos de turbulência, passo de tempo com estudo sobre o CFL, subrelaxamento, validação da malha e redução da quantidade de simulações a partir do DoE. Com base nas referências bibliográficas, foi realizado uma seleção dos artigos que forneciam dados para cada parâmetro. 13 • Domínio computacional cilíndrico: Em um estudo para investigação numérica do efeito de torres e naceles na esteira próxima de um modelo de turbina eólica de eixo horizontal foi utilizado o domínio (2D+8D), 6.7D, 6.7D, respectivamente em X, Y, Z [4], sendo D o diâmetro das pás da turbina. Outro estudo, para investigações numéricas e experimentais de previsões de Horizontal Axis Wind Turbine (HAWT) perto de esteira usando partículas, foi utilizado um domínio cilíndrico de 4.5D + 5.5D axialmente e 5D de diâmetro [5]. Entretanto, a análise do projeto do grupo é pontual e esses domínios apresentados são maiores do que o necessário para o estudo. Como o objeto de estudo se localiza logo atrás da turbina, não faz sentido gastar esforço computacional com regiões distantes. • Malha: A aplicação realizada por um estudo recente sobre avaliação de modelos comerciais em CFD, utilizou o Rigid Body Motion, que leva em consideração o movimento das pás do rotor, obtendo assim um resultado mais preciso do direcionamento do vento após a passagem pelo rotor [6]. Como esse fenômeno também ocorre no projeto do grupo, a aplicação foi similar. A área em volta do anemômetro, do rotor e na camada limite nas superfícies sólidas sofreu maior refinamento, dada a sua importância nos resultados desse trabalho. O refinamento da malha foi avaliado para atingir a sua independência, ou seja, para que ela não interferisse nos resultados das simulações. O formato da malha é poliédrico, por se tratar de uma figura menos complexa e, se comparado à tetraédrica, a primeira exige menos esforço computacional. Além disso, as células são ortogonais ao escoamento de ar, sendo mais apropriado para modelos rotacionais. • Condições de contorno Após estudo prévio a partir da utilização de papers sobre avaliação de modelos comerciais em CFD e modelamento de turbinas eólicas, dentre as condições de contorno analisadas, foi preciso utilizar na entrada do domínio: velocity inlet; na saída do domínio: pressure outlet; e, nas laterais do domínio: simmetry plane. A velocidade de entrada e a rotação estão relacionadas e foram simuladas com valores distintos durante à etapa de após a definição da independência da malha. Ficou estabelecido que a pressão de saída possuía condição ambiente [8]. 14 Vale ressaltar que para o velocity inlet, foi considerada a simplificação de que a velocidade do vento é uniforme. Sabe-se que a velocidade do vento varia conforme a altura em relação ao solo, como evidenciado na Figura 3 [2], mas como não foi possível fazer a verificação e validação desse perfil de velocidade, foi considerada esta simplificação a fim de otimizar o projeto. Figura 3 - Relação velocidade do vento x superfície. Somado a isso, algumas considerações iniciais foram aplicadas: 1. Todas as medições do vento são baseadas em planos horizontais em relação ao solo, portanto bidimensionais. Isso significa que o inflow é zero. 2. A velocidade medida pelo anemômetro é a velocidade absoluta no plano horizontal. 3. A referência da direção do vento é móvel, pois ela é medida em relação à nacele, que pode se mover. • Modelo de turbulência O modelo Large Eddy Simulation (LES) exige requisitos computacionais maiores, se comparado com modelos de duas equações, como SST- kω e realizable-Kɛ [8]. No que diz respeito ao modelo Reynolds Stress Model (RSM), o LES pode precisar de cerca de duas vezes os recursos computacionais para o mesmo cálculo [7]. O modelo SST-Kω apresenta bons resultados para escoamentos simples. Embora o RSM seja considerado um dos mais complexos, ele descreve todas as propriedades médias do escoamento e as tensões de Reynolds sem ajustes individuais. Porém 15 o RSM não é tão bem validado, se comparado ao realizable-Kɛ e SST -kω, pois possui um alto custo nos cálculos, se tornando menos utilizado. Dentre os quatro modelos supracitados, o realizable-Kɛ é o que mais se adequa ao objetivo de estudo do projeto, dado que exige uma carga computacional menor, e consegue predizer com eficácia os efeitos de turbulência sobre a região do anemômetro [7,5]. Em uma análise mais aprofundada, o modelo de turbulência Reynolds Stress Model (RSM) representa o modelo clássico mais completo de turbulência. É considerado um fechamento de turbulência de nível superior. Esse fechamento tem o nome muitas vezes de Second Order Closure. O modelo evita a abordagem de viscosidade turbulenta e parte diretamente para o cálculo dos componentes do tensor de estresse de Reynolds. Esse modelo depende diretamente da equação de transporte de tensão de Reynolds. Ele é capaz de explicar interações complexas em campos de fluxo turbulento e os efeitos direcionais das tensões de Reynolds. Comparando o modelo RSM com outros dois modelos que são baseados em viscosidade turbulenta, o SST-K Ômega e o Realizable Képsilon, é possível observar que existem deficiências desses dois modelos quanto em fluxos turbulentos complexos da vida real que são frequentemente encontrados em aplicações de engenharia. Em escoamentos com alto grau de anisotropia, curvatura significativa da linha de corrente, separação de escoamento, escoamentos com zonas de recirculação ou escoamentos influenciados por efeitos rotacionais médios (que é o caso deste projeto), o desempenho desses modelos é insatisfatório. Os fechamentos dos modelos que se baseiam em viscosidade turbulenta não podem explicar os estados limitantes de turbulência. Em escoamentos turbulentos em decaimento, observou-se que o estado de turbulência se aproxima rapidamente de um estado isotrópico com equipartição da energia cinética turbulenta entre os componentes das tensões de Reynolds. Modelos baseados em viscosidade turbulenta nunca podem replicar esse comportamento de retorno à isotropia, isso porque não podem replicar o comportamento de escoamentos turbulentos no limite de Distorção rápida (Rapid Distortion), onde o escoamento se comporta como um meio elástico (ao invés de viscoso). O modelo de turbulência SST k-ômega é um modelo de viscosidade turbulenta de duas equações. O uso de uma formulação k-ômega nas partes internas da camada limite torna o modelo diretamente utilizável até a parede através da subcamada viscosa. Dessa forma o modelo pode ser usado como um modelo de baixa turbulência sem funções de amortecimento. A formulação SST também muda para um comportamento do modelo de turbulência k- épsilon na corrente livre, o que evita o problema do k-ômega de ser muito sensível às 16 propriedades de turbulência na corrente livre de entrada. Pode ser apontado que o SST k-ômega tem um bom comportamento em gradientes de pressão adversos e em fluxos de separação. O modelo produz níveis de turbulência grandes demais em regiões com significativa deformação normal, regiões de estagnação e regiões com forte aceleração. O modelo de turbulência Realizable Képsilon é um modelo de duas equações, ou seja, que inclui duas equações de transporte extras que representam as propriedades turbulentas do fluxo. Isso garante que seja possível, o modelo levar em conta efeitos históricos como a convecção e a difusão de energia turbulenta. Das equações a primeira é referente à energia turbulenta cinética, k, que é a variável que determina a energia na turbulência. A segunda é com relação a dissipação turbulenta, épsilon, que é a variável que determina a escala da turbulência. O modelo K-épsilon é utilizado majoritariamente para fluxos de camada de cisalhamento livre com gradientes de pressão relativamente pequenos. A precisão do modelo diminui quando existem fluxos contendo grandes gradientes de pressão adversos. Pode se inferir então que o modelo é inadequado para compressores e entradas. • Passo de tempo (time step) A determinação do time step é comumente realizada a partir do tempo que o rotor leva para percorrer a variação de 1° (um grau) [9, 8, 7 e 10] e é avaliada pelo número de Courant- Friedrichs-Lewy (CFL). O CFL é uma condição para a convergência de métodos numéricos instáveis que modelam fenômenos de convecção ou onda e, portanto, desempenha um papel importante nas simulações CFD. Em resumo, o número de CFL pode ser descrito a partir de (1), 𝛥𝑡 (1) 𝐶𝐹𝐿 = 𝑉 𝛥𝑥 que 𝑉 é a magnitude da velocidade, 𝛥𝑡 o time step e 𝛥𝑥 é o tamanho da menor aresta da forma que compõe a malha. Em uma simulação de CFD, o CFL indica o quanto a informação viaja através de uma célula da malha computacional em um passo de tempo. Se o CFL for maior que 1, significa que a informação se propaga por mais de uma célula da malha a cada passo de tempo, tornando a solução potencialmente imprecisa e levando a resultados não físicos ou à divergência da solução em determinados esquemas de integração [12]. Por conta disso, para poder ser aplicado o CFL que mais se adequava ao projeto, foram realizadas simulações teste. 17 O cálculo do time step, foi feito utilizando a velocidade de entrada do vento fornecida pela WEG que é 8 𝑚/𝑠. Como a velocidade de giro do rotor da turbina é 10,7 𝑅𝑃𝑀, convertendo este valor em graus/s, obtém-se 64.2 𝑔𝑟𝑎𝑢𝑠/𝑠. Utilizando a variação de um 1° como parâmetro de cálculo do time step é possível concluir que o rotor precisa de 0,015576 𝑠 para atingi-la. Para validar se este valor pode ser utilizado ou não, foi preciso verificar a condição do CFL. O STAR-CCM+ possui um monitoramento de CFL para mostrar se o valor aplicado está correto e é possível ativar uma configuração automática do CFL ao trocar o modelo da simulação de segregated flow para coupled flow, reduzindo assim o esforço de cálculos. Além disso, o time step também pode ser adaptativo ao ativar a função Adaptative Time-Step na sessão Models da simulação. Utilizou-se inicialmente o valor calculado de 0,015576 s mas a aplicação da calibração dos valores pelo próprio STAR-CCM+ serviu para otimizar todo o processamento necessário nas simulações. • Sub relaxamento O STAR-CCM+ usa o método numérico de Gauss-Seidel para a resolução do sistema de equações matriciais nas simulações. O grupo decidiu utilizar os valores padrão de sub relaxamento do software como parâmetro inicial. Os valores são os seguintes: o Velocidade: 0,8. o Pressão: 0,2. o Turbulência: 0,8. o Viscosidade turbulenta: 1,0. Esses valores geram uma curva com os resíduos da simulação que estão demonstrados na Figura 4. Figura 4 - Gráfico de Resíduos. 18 Pelo gráfico, observa-se que era necessário realizar ajustes nos valores para diminuir a amplitude do resíduo de alguns parâmetros como o da turbulência e, portanto, inúmeras análises foram aplicadas a partir de simulações, pesquisas, alterações e mais simulações. Além dos valores de sub relaxamento, foi possível fazer a alteração do número de maximum inner iterations que tem como padrão 5. Essas iterações internas representam a convergência do método de Gauss-Seidel em um time-step. • Validação da malha Uma malha de baixa qualidade pode resultar em uma redução significativa na precisão e eficiência das simulações em CFD. O STAR-CCM+ fornece a capacidade de verificar a validade da malha e obter estatísticas sobre a qualidade da malha, através dos seguintes parâmetros: ângulo de assimetria da célula, métrica de qualidade da célula, métrica de alteração de volume e indicador de qualidade Chevron. Por isso, fez-se necessário verificar a qualidade dela para poder ter uma acurácia ainda melhor. • DOE O Design of Experiment (DOE) é uma técnica para se planejar experimentos, ou seja, definir quais dados, em que quantidade e em que condições devem ser coletados durante um determinado experimento, buscando, basicamente, satisfazer uma precisão estatística na resposta. Seu emprego no projeto serviu como forma de reduzir a quantidade de simulações, otimizando o tempo de execução do grupo nas análises. 2. Metodologia 2.1 Cálculo da camada limite para geração da malha próxima às superfícies sólidas Para garantir que não houvesse interferência na leitura do anemômetro foi calculado o tamanho da camada limite, como pode ser observado na Figura 5 que será gerada nas pás, na nacele e nos demais elementos sólidos da turbina. 19 Figura 5 - Camada limite [11]. A equação utilizada para a determinação da camada limite [8] pode ser verificada em (3): 𝑥.𝜌.𝑢 𝑅𝑒 = ∞ 𝜇 (2) 0.16𝑥 𝛿 = (3) 𝑅𝑒1/7 O regime de escoamento em uma placa plana deixa de ser laminar e passa para turbulento com 𝑅𝑒 > 106. Entretanto, apesar das superfícies não serem uma placa plana, o próprio movimento do rotor da turbina já faz com que o escoamento seja turbulento. Os valores utilizados no cálculo e o resultado obtido estão descritos na Tabela 3. Tabela 3 - Cálculo da camada limite. Variável Nacele Rotor 𝑥 10 m 3 m 1,18167 kg/m³ 1,18167 kg/m³ 𝜌 𝑢∞ 8 m/s 8 m/s 𝜇 1,825.10−5 1,825.10−5 Re 5,18.106 1,55.106 𝛿 0,18 m 0,06 m 20 Este valor de camada limite foi considerado para a criação da malha. 2.2 Domínio Tendo em vista a semelhança dos estudos que foram realizados nesse projeto e o estado da arte [4,5], definiu-se para o domínio utilizar as seguintes dimensões a partir do eixo do rotor da turbina: 1D+4D (eixo X), 4D (eixo Y) e 4D (eixo Z), onde D é o diâmetro do rotor da turbina. A Figura 6 representa a geometria cilíndrica que foi utilizada como domínio computacional. A entrada e a saída do domínio são representadas pelas bases do cilindro. Já a origem da coordenada está localizada dentro do hub conforme Figura 7. Figura 6 - Geometria cilíndrica. Figura 7 - Origem da coordenada. 21 O sistema de coordenadas possui origem no centro do rotor. O eixo X trata-se da direção axial, portanto a direção de escoamento do vento, possuindo 1D à montante e 4D à jusante. O eixo Y teve um domínio especificado de 4D, e 4D também foi estabelecido radialmente, representando o eixo Z. Isso pode ser mais bem visualizado a partir da Figura 8. Figura 8 - Domínio computacional Há duas partes necessárias para a simulação numérica: o cilindro que representa o domínio computacional e outro que representa as partes móveis. Com base nessas duas partes cilíndricas criadas, foram construídas duas partes virtuais que representam o escoamento de ar. Assim, o cilindro que representa as partes móveis precisou ser subtraído das partes sólidas de seu interior, o rotor e parte do hub2. Da mesma forma, o cilindro que representa o domínio precisou ser subtraído da parte virtual móvel e das partes sólidas. Dessa forma essas duas partes, juntamente com os resfriadores ativo e passivo, foram transformadas nas regiões de análise. As regiões que se referem ao escoamento de ar, como a região móvel e a estática, possuem uma física apropriada ao que era esperado avaliar na simulação. Inicialmente, definiu-se que o domínio de análise era tridimensional. A seguir, informou-se que o escoamento era transiente e que o ar é um gás que escoa a densidade constante devido a não haver efeitos de variação da temperatura significativos e a velocidade ser baixa o suficiente para que o número de Mach fique abaixo de 0,3. Um dos modelos numéricos que foi utilizado nas simulações se refere ao que resolve as equações de Navier-Stokes. Inicialmente foi utilizado o modelo não acoplado, ou seja, um modelo que resolve as equações separadamente e faz a ligação entre elas através de 22 uma aproximação baseada em um algoritmo preditor-corretor. A seguir escolheu-se o modelo de turbulência Realizable Képsilon. 2.3 Criação da malha Na criação da malha alguns fatores foram essenciais a serem considerados. O primeiro se refere a região onde está a turbina e o anemômetro, como pode ser observado na Figura 9, representada por um corte do cilindro. Figura 9 – Vista lateral da malha. A região em destaque se refere à nacele, rotor e hub. O anemômetro está representado pela parte mais escura, na qual a malha é ainda mais refinada. O segundo fator se refere ao rotor. Como ele causa a turbulência do ar, essa parte possui uma malha mais fina. Segundo o cálculo da camada limite existe uma região de escoamento laminar de 0,18 m ao redor da nacele, e de 0,06 m ao redor das pás. A malha próxima à superfície é denominada prism layer, ela possui a mesma espessura da camada limite e foi configurada para ter 15 camadas, stretching de 1,1 e possuir um refinamento de mais de 25 milhões de células. Há duas malhas, uma para a parte rotativa e outra para a parte estática. A malha rotativa usa o surface remesher, automatic surface repair, trimmed cell mesher e o prism layer mesher para a camada limite, ela tem dois custom controls. A malha estática usa as mesmas configurações para criação da malha e possui dois surfaces controls e um volumetric control. Na malha rotativa, o primeiro surface control é para o rotor e especifica um target surface size e minimum surface size de 5% do valor base da malha, ajustado em 1m. O segundo 23 é para rotating domain e especifica um target surface size e minimum surface size de 15% do valor base da malha. Na malha estática, o primeiro surface control é para a entrada, saída e superfície lateral do domínio e especifica um target surface size e minimum surface size de 400% do valor base da malha, ajustado em 1m. O segundo é para rotating domain e especifica um target surface size e minimum surface size de 15% do valor base da malha. O volumetric control é para uma região (block) entre o rotor e os anemômetros, conforme a Figura 10 e especifica um custom size de 5% do valor base da malha. Figura 10 – Custom size. Por padrão as regiões no software são estacionárias, porém para considerar o movimento da superfície é utilizado o Rigid Body Motion. 2.4 Condições de contorno Na entrada do cilindro, que representa o domínio computacional, foi aplicado o velocity inlet, especificado de acordo com as orientações fornecidas pela WEG, sendo sugerido inicialmente uma velocidade de vento de 8 m/s. Na saída foi utilizado pressure outlet. E ao redor do cilindro, foi utilizado simmetry plane. 2.5 Parâmetros para a simulação Para as simulações no STAR-CCM+ foram utilizados os dados enviados pela WEG, exibidos na Tabela 4 , além dos oriundos da revisão bibliográfica. 24 • Dados enviados pela WEG: Tabela 4 - Dados enviados pela WEG. Velocidade do hub [m/s] 6 8 10 Rotação [RPM] 8.3 10.7 11 Inflow [graus] -8, 8 e 0 Densidade do ar [kg/m³] 1,07 e 1,18 Turbulência longitudinal 1 0,10 Turbulência longitudinal 2 0,24 0,20 0,18 Turbulência longitudinal 3 0,19 0,17 0,16 Direção do vento [graus] -15 a 15 A velocidade do hub é a velocidade do vento na entrada do domínio e as turbulências citadas são as intensidades turbulentas. Na Tabela 4 percebe-se que o inflow é o ângulo de ataque vertical do vento no perfil da pá da turbina, como pode ser observado na Figura 11. Figura 11 - Vista superior e lateral da turbina eólica para visualização da direção do vento horizontal e vertical (autoria da WEG). Vale ressaltar, que embora só há dados das densidades do ar, foi considerado também, a temperatura de 25°C. Em relação à quantidade de simulações, foi alinhado com a WEG que se fossem consideradas 3 velocidades do vento, 3 inflow, 2 densidades do ar, 3 turbulências e 3 direções do vento, resultaria em um total de 162 simulações, o que tornaria o processo inviável. Além disso, seriam necessárias de 20 a 30 simulações inicialmente para acertos do modelo e independência de malha. Como estima-se que cada simulação leve em média 1 dia para ser finalizada em um cluster do Insper, seria preciso mais de 162 dias para realização de todas, 25 sendo inviável no que diz respeito ao tempo de elaboração do PFE, adicionadas às dificuldades que o grupo enfrentou e enfrenta com as instabilidades do servidor de licença e a VPN do Insper. Ademais a WEG forneceu um intervalo de tempo em que os dados de vento podem ser usados para a validação no anemômetro. Para a validação da velocidade do vento, já que ele está a uma altura igual à linha de centro do rotor, foi utilizado o campbell_01. Além disso, foi utilizado o campbell_08 que está um pouco abaixo da altura para a validação da direção do vento. Não foi preciso plotar todos os dados como os enviados pela WEG, porque o vento muda de direção e a rotação do rotor muda. A planilha enviada não é de dados com a mesma velocidade do vento. Por fim, foi acordado que o intervalo de direção deveria estar no range de 25°. • Aplicações dos dados da revisão bibliográfica: Durante a revisão bibliográfica, foi possível ter acesso a um vídeo do canal do youtube “Applied Computational Fluid Dynamics” [3] no qual continha uma simulação seguindo as seguintes etapas: 1. Definição do Domínio Computacional. 2. Definição das Regiões. 3. Geração de malha. 4. Definição de Física e Condições de Fronteira. 5. Definição de Monitores, Configurações do Solver e Critérios de Parada. 6. Pós-processamento . Com base nesse vídeo, foram anotados, a partir da tabela no Apêndice A, as principais informações que serviram como parâmetros nas análises futuras do desenvolvimento das simulações. Cada parâmetro a ser considerado, precisou estar adequado com a malha, modelo de turbulência, condições de contorno e domínio computacional. 26 2.6 Simplificações do modelo CAD fornecido pela WEG O modelo CAD enviado pela WEG apresentava uma série de informações que não foram necessárias para a realização da simulação por CFD. Por esse fato, foi preciso realizar inúmeras simplificações para garantir que fosse possível avançar com as simulações no STAR- CCM+. Essas simplificações foram feitas nas pás, nos conectores das pás, no hub, nos anemômetros, nos coolers, na nacele e na torre. Para as simplificações foi necessário a remoção de parafusos que estavam presentes em toda a estrutura, retirada de itens sem relevância para a simulação no anemômetro, retirada dos suportes (mãos francesas) nos coolers além da eliminação de superfícies abertas. Após a simplificação da geometria, foram criadas algumas partes para a turbina eólica com o intuito de agrupar as peças. Problemas apareceram durante a geração da malha para as simulações, como a malha sendo gerada dentro do hub, inclinação do eixo do rotor que estava inadequada, nacele e hub com partes ocas e cooler inadequado para a simulação. Para tanto, inicialmente foi realizado tentativas para fechar as partes ocas do hub e da nacele. Para isso, utilizou o comando extend solid do STAR-CCM+. Foi preciso criar uma superfície que ao ser estendida fechasse o sólido. • Inclinação do eixo do rotor No que diz respeito a inclinação que não estava adequada, como pode ser visualizado na Figura 12, foi aplicado uma inclinação em x de 6° e em seguida, determinou a origem da parte rotativa. 27 Figura 12 - Problemas na inclinação. Para o cooler, foi preciso criar um CAD novo, dado que o que fora criado, estava com muitos problemas. • Alongamento das pás Para acertar o hub, foi preciso manter os conectores das pás e estendê-las até que encontrassem uns nos outros, como pode ser observado na Figura 13. Figura 13 - Hub adequado. 28 • Cooler Além desses problemas, foi observado também uma variação do tamanho fora do cooler, como mostrado na Figura 14 , foi preciso tirar o prism layer do cooler e diminuir o tamanho da célula dentro dele para ter mais células dentro. Figura 14 - Cooler Em relação à nacele o gerador, eles estavam maciços no interior. Entretanto, ao gerar a malha, o resultado foi o evidenciado na Figura 15. Figura 15 - Malha gerada com componentes maciços. 29 A partir disso foi possível perceber que havia um problema com a nacele e com o hub. Após estudo e orientações, foi possível gerar um arquivo com as partes adequadas. Com o arquivo gerado foi preciso criar uma interface entre as duas partes rotacionais. • Vector scene Uma outra etapa importante foi transformar o cooler em meio poroso. Para isso, foi preciso construir o vector scene no STAR-CCM+. Foram utilizados dois tutoriais fornecidos pelo software como guia, sendo eles: Porous Resistance: Isotropic Media e Porous Resistance: Orthotropic Media. O primeiro auxiliou a como criar o meio poroso. Para isso, foi preciso definir as propriedades do material necessárias para os modelos selecionados. Criar interfaces de região, definir condições de contorno, especificar coeficiente de porosidade. Ademais, foi preciso realizar a configuração dos critérios de parada, com base em valores residuais e quantidades monitoradas. E por fim, analisar os resultados da solução usando os recursos de visualização do STAR-CMM+, através da criação de uma cena vetorial, a fim de visualizar o campo vetorial de velocidade à medida que a solução se desenvolve em um plano que divide as regiões fluida e porosa. A região porosa foi definida por seus coeficientes de resistência inercial e viscosa, que são 35 kg/m^4 e 1500 kg/m³, respectivamente. O segundo auxiliou em como configurar um meio poroso não isotrópico. A partir dele, foi possível garantir que o fluido na região porosa não fluísse em nenhuma direção que não fosse a direção do fluxo em massa. Na Figura 16, é possível observar a vista superior. O ponto em branco é a pá. Essa vista leva em consideração o “snap to part”, ou seja, o corte feito no centro do anemômetro. 30 Figura 16 - Vector Scene. Além desses problemas e resoluções, foi preciso alterar o bloco dos anemômetros para uma melhor análise. • Bloco dos anemômetros Depois de pesquisas e orientações, foi possível constatar que o bloco dos anemômetros deveria ir além da parte rotativa. Foi preciso depois realizar o subtract no software, como pode ser observado na Figura 17 , para ter uma boa interface entre eles. Para garantir dessa maneira, que ele chegasse até a parte rotativa e não além. Vale destacar que não precisou subtrair os anemômetros, pois a malha deles sobrepõe essa. 31 Figura 17 - Bloco dos anemômetros. • Malha No processo de refinamento da malha, etapa importante para conseguir um processamento melhor do computador de forma a garantir que os pontos de interesse, tivessem um foco de atenção melhor, houve alguns problemas. Figura 18 - Tentativa de refinamento da malha. Em uma das primeiras tentativas de refinamento, como evidenciado na Figura 18, a malha no centro era mais grossa que na fronteira. E, de acordo com o objetivo do projeto, deveria ser o contrário. Foi criado duas interfaces na parte rotativa. Deveria ter sido criada apenas uma interface. Foi preciso fazer a parte superior do domínio onde está o anemômetro mais refinado e a parte de baixo menos refinada, como evidenciado na Figura 19, pois o escoamento da parte de baixo não influenciava na medição nos anemômetros. Para alterar o crescimento/tamanho da malha, foi utilizado o comando custom control no software. 32 Figura 19 - Malha refinada. O ideal do refinamento é deixar “fino” nos pontos de análise mais importante e “grosso” nos demais. Inicialmente, o refinamento precisava estar onde se dá a medição mais importante no domínio: os anemômetros. Além disso, deveria estar refinado no rotor, pois ele é que altera o sentido de fluxo na entrada do domínio. Abaixo há a Figura 20 que evidencia o refinamento inicial aplicado no anemômetro. Figura 20 - Refinamento no anemômetro A malha frontal precisava ficar independente da presença das pás. Ela precisava ser mais refinada na parte superior da nacele e menos abaixo. Onde não tem pá fica mais grosseira, como evidenciado na Figura 21. 33 Figura 21 - Refinamento na parte superior da nacele superior à parte inferior. A análise de qual malha utilizar é função do comportamento da velocidade e direção do vento no anemômetro. Foi preciso plotar os gráficos para ajudar na análise. 2.11 Independência de malha e boas práticas aplicadas. Foi preciso rodar, inicialmente as seguintes simulações, com base nos dados da Tabela 5 , para avaliar a independência de malha: 34 Tabela 5 - Primeira tentativa de independência da malha. Além desses dados, que contém o subtract da região entre os anemômetros e o rotor menos a parte rotativa, a malha que representa cada anemômetro em si tem 0,02m de base size. O domínio rotativo e a região dos anemômetros (entre os anemômetros e o rotor) devem ter o mesmo tamanho de célula. Os anemômetros em si sempre 0,02m. O objetivo de analisar as simulações com diferentes tamanhos de malha para o domínio, e região rotativa, é para escolher o tamanho de malha com resultado mais adequado. A análise de qual malha utilizar foi em função do comportamento da velocidade e direção do vento, além da comparação dos resíduos gerados pelas simulações. Para tanto, foi preciso plotar gráficos de direção e velocidade do vento no anemômetro, e comparar as curvas dos diferentes tamanhos de malha. • Boas práticas aplicadas. 1° Disponibilizar a simulação que foi considerada boa para o primeiro caso, para os demais e, realizar apenas ajustes de parâmetros no CAD. Vale lembrar que quanto menor a malha, mais tempo resultaria. 2° Avaliar os resíduos e depois verificar se houve diferença significativa entre os modelos de turbulência. 3° Rodar as simulações com os demais parâmetros da WEG, analisando o DoE. 4° Colocar a posição da pá em frente ao anemômetro 1. Para isso, foi necessário um pós-processamento no excel para trabalhar com o arquivo .csv exportado do STAR-CCM+ de cada simulação. Foi preciso ter atenção para a largura da pá, que precisou abranger apenas o tempo em que a pá passa pela frente do anemômetro 1. Dado que a pá tem a variação de 1°, foi 35 preciso saber o ângulo que a entrada e saída de cada pá deveria fazer para ficar na frente do anemômetro 1 e, para isso, saber a largura da pá na altura em que está o anemômetro. 5° Plotar posteriormente as curvas de magnitude de velocidade no plano XY e direção do vento no plano XY, de diferentes simulações, para analisar qual fazia mais sentido ao objetivo do projeto. Nesta etapa, para criar a Field Function que analisa a direção no plano XY foi utilizado: (𝑎𝑡𝑎𝑛(${𝑀𝑎𝑥𝑖𝑚𝑢𝑚 − 𝑎𝑛𝑒𝑚ô𝑚𝑒𝑡𝑟𝑜1 𝑗}/${𝑀𝑎𝑥𝑖𝑚𝑢𝑚 − 𝑎𝑛𝑒𝑚ô𝑚𝑒𝑡𝑟𝑜1 𝑖})) 6° Utilizar a função de autosave no STAR-CCM+ para não perder a simulação em sua totalidade. Alguns dos primeiros gráficos da simulação podem ser observados na Figura 22 e na Figura 23. Figura 22 – Reports Plot. Figura 23 - Direção do vento 1. 36 A partir daqui a análise ficou limitada sem a posição da pá em frente ao anemômetro 1, que nesta análise ocorreu em um total de 9 vezes, considerando todas as rotações. Com o passar do tempo, com um valor definido, a quantidade de vezes foi adequado. Nas demonstrações futuras, a recomendação era que a largura da pá devia abranger apenas o tempo em que a pá passa pela frente do anemômetro. Para colocar as pás, foi preciso acompanhar a simulação para ver em que time step a pá entrava na frente do centro do anemômetro. Para isso, uma aplicação foi observar a malha de frente para o rotor e clicar, pelo STAR-CCM+ na parte do anemômetro 1 para ficar mais claro de visualizar, conforme pode ser observado no zoom na Figura 24. Figura 24 - Zoom do bloco do anemômetro 1. Para saber em qual time step a pá estaria alinhada ao anemômetro, foi preciso observar o ângulo que o anemômetro 1 fazia com a vertical. Além disso, foi preciso saber o ângulo de entrada e de saída que cada pá deveria fazer para ficar na frente do anemômetro 1 e, para isso, saber a largura da pá na altura em que estava o anemômetro. No total são dois times steps por pá por rotação: o time step em que a pá fica na frente e sai. 37 Foram feitas mais simulações. Uma boa prática adotada aqui foi comparar arquivos prontos, uns com os outros, para poder visualizar as diferenças e filtrar os parâmetros mais adequados. Um dos gráficos plotados pode ser observado na Figura 25. Figura 25 -Ângulo de Yaw. Foi levado para plotar esse gráfico a função “média móvel” no software. Entretanto, aqui houve um problema com o azimute. No gráfico da Figura 26, não foi levado a média móvel como parâmetro. Figura 26 - Ângulo de Yaw sem média móvel. Nestes gráficos, cada dado foi retirado com um passo de tempo de 1 segundo. Foi alinhado, nesta etapa, que o intervalo de direção adequado, deveria estar no range de 25°. Além disso, a WEG havia enviado alguns dados, mas foi acordado que nem todos seriam utilizados, pois o vento mudava de direção e a rotação do rotor também alterava. Logo, foi combinado de utilizar os dados do campbell_01 para velocidade e o campbell_08 para direção do vento, que são colunas presentes nesse arquivo de envio. 38 A partir daqui, foram realizadas mais simulações. O resultado de uma delas pode ser observado na Figura 27 e na Figura 28. Figura 27 - Zoom do gráfico de velocidade em m/s. Figura 28 - Gráfico de velocidade em m/s. Nessa forma de visualizar, os campbell estavam no mastro localizado a 300 m à montante da turbina. Assim, se o vento estivesse a 10m/s ele levava 30s para chegar à turbina e um pouco mais até o anemômetro. Foi preciso então considerar a velocidade variável do vento e o tempo que ele leva até o anemômetro para fazer o gráfico comparativo. 39 Após mais simulações, foi possível observar o número de Courant, Friedrichs e Lewy (CFL) estava com valor 50, que representa o default no software. Foi dado sequência a um estudo para analisar a influência do CFL no passo de tempo da simulação. Dentro das análises realizadas, foram feitas simulações com o valor 50 em comparação com 1, 0.7 e com o 0.5, como pode ser observado na Figura 29. Figura 29 - CFLs aplicados distintos. A ordem 2,3 e 1 é a contagem no sentido de rotação das pás. Através desse gráfico, as curvas azul e roxo, 0.5 e 50, foram as que sofrerão oscilações mais acentuadas. A 0.5 levou mais tempo para estabilizar, se comparado com as demais. Entre a 0.7 e a 1, a escolha foi seguir com o CFL igual 1, dado que apresentou um comportamento mais ordenado após a passagem da pá, do que a curva que apresentou 0.7 de CFL. Além disso, foi possível notar que elas possuem efeitos opostos. A análise ficou limitada a Figura 29, dado que como as curvas oscilam muito, o gráfico fica carregado de informações e impacta no entendimento. Logo, foi preciso defasar as curvas, de maneira a ajudar melhor na compreensão e, o resultado pode ser observado a partir da Figura 30. 40 Figura 30 - Curva defasada. Foram utilizados 40 graus, para a defasagem entre às curvas verde e vermelha, de modo a aproximar os picos e vales, a fim de visualizar melhor a diferença. Embora ambas apresentam oscilações similares, o 0,7 tinha uma amplitude mais acentuada, gerando mais instabilidades na análise. A partir dessa constatação, todas as simulações passaram a serem utilizadas com CFL igual a 1, para obter assim uma melhor análise. Essa análise também foi feita a partir dos resíduos, caso a caso, sendo um exemplo evidenciado na Figura 31. Figura 31 - Resíduos. A partir dela também pode constatar que a convergência acontece. Caso a simulação rode por mais tempo, a velocidade de convergência deixa de ser fator fundamental. Se a 41 convergência se dá para o mesmo valor da direção do vento, então qualquer um dos resultados fica bom. Utilizando o CFL igual a 1, foi possível observar mais problemas. Um deles pode ser observado a partir da Figura 32 e da Figura 33. Figura 32 – Reports Plot. Figura 33 - Direção do vento. Foi possível constatar que as curvas não apresentavam um bom resultado, pois a amplitude da direção do vento estava bastante acentuada. Ocorreu problemas na análise do resíduo também, a partir da Figura 34. Figura 34 - Resíduos. 42 A partir dela, foi possível constatar que ocorreu um problema no pico na iteração 3000- 3500. Além disso, os resíduos estavam, em alguns pontos, acima de 1, apresentando uma amplitude indesejada. Para isso, foi modificado o parâmetro implicity unstead para diminuir as células com alto resíduo. Além disso, foi alterado threshold que é uma variável para “turbulent dissipation rate residual”, que estava precisando de um refinamento maior em toda a área em destaque, como pode ser observado na Figura 35. Figura 35 - Thresold. Após melhorias, foram realizadas mais simulações. Nestas, na criação dos gráficos, foram feitas comparações com as médias e com e sem os outliers, que são os pontos que ficam longe da média das velocidades. Nos gráficos de direção, não é tão grande a diferença. Um dos exemplos pode ser observado a partir da Figura 36. 43 Figura 36 - Presença de outlier. Após, foram realizadas mais análises a partir de diferentes simulações. Para escolher a simulação base, que posteriormente foi utilizado para o estudo do DOE, foi utilizado os dados da Tabela 6 com tamanhos de células distintos para cada malha do sistema. Tabela 6 – Parâmetros para simulações. Identificação da Tamanho da célula Tamanho da célula Tamanho da célula simulação do domínio estático domínio rotativo nos anemômetros 0 3,0 0,5 0,02 1 2,0 0,5 0,02 2 2,0 0,4 0,02 3 2,0 0,3 0,02 4 1,0 0,3 0,02 5 1,0 0,2 0,02 6 2,0 0,2 0,02 Foram rodadas mais de 35 simulações. Precisou realizar a comparação das células do domínio e do rotativo. Inicialmente, foram rodados apenas com 10 voltas, com a análise feita a partir da terceira volta, que até o momento, era o que acreditava ser o parâmetro para estabilização dos resultados. Como, pode ser observado, por exemplo na Figura 37. 44 Figura 37 - Presença de instabilidades. É possível perceber que há instabilidades no início até o início da estabilização. Entretanto, posteriormente, optou por realizar 15 voltas em cada simulação e, analisar os resultados, que serão apresentados no tópico 4 deste documento, a partir da 11ª volta. Foram realizadas mais análises, mas agora considerando as curvas com: 1ª (sim.15, pelo gráfico) Com domínio de 5m e 0.2m de rotativo. 2ª (sim.3, pelo gráfico) Com domínio de 10m e 0.3m de rotativo. 3ª (sim.13, pelo gráfico) Com domínio de 15m e 0.5m de rotativo. 4ª (sim.14, pelo gráfico) Com domínio de 20m e 0.7m de rotativo. A sim.15 e a sim.14, se comparada com as demais, não resultaram em um bom comportamento. Pelo gráfico, apresentaram amplitudes indesejadas na parte superior e inferior. A sim.3 e sim.13, foram as mais adequadas. Entretanto, a sim.13 ainda se demonstrou mais controlada. A sim.3 apresentou picos e vales altos, se comparado, como pode ser observado em um zoom (através de um corte) da Figura 38. 45 Figura 38 – Corte do gráfico para evidenciar à relação velocidade versus azimute das diferentes simulações. Além disso, foram feitas análises do resíduo também. Àqueles que ficaram acima de 1, foram descartados, dado que foge do escopo do projeto e, influencia negativamente nos picos e vales. Vale lembrar que em todas essas decisões, foi considerado que o vento entrava de forma perpendicular. Para a decisão de desconsiderar as primeiras voltas e trabalhar com as últimas, também foram realizadas simulações e foi verificado que o valor se mantinha parecido, como por exemplo, pode ser observado (pela curva vermelha tracejada, que representa a média, com os mesmos parâmetros da sim.13) na Figura 39. Figura 39 - Média em relação a sim.13 46 Inicialmente foi rodado com 5 voltas. Após 10 e por fim, 15 voltas para analisar a estabilidade. Considerando todos os parâmetros analisados, àquela que fazia mais sentido em relação ao projeto para ser a simulação base, foi a 13, com 15m no domínio e 0,5m no rotativo. 2.12 Modelos de turbulência Após a definição da independência de malha foi preciso escolher definitivamente o modelo de turbulência a ser utilizado. O grupo testou três opções: RSM, SST- kω e realizable- Kɛ para avaliar qual deles possuiria a melhor performance. Vale citar que até o momento todas as simulações foram feitas utilizando o realizable-Kɛ graças aos papers terem o indicado como melhor correlação com os dados experimentais. 2.13 Escoamento do ar na pá Uma etapa importante foi o entendimento sobre o escoamento do ar na pá da turbina e a influência na leitura do anemômetro que fica localizado na nacele. Nas Figura 40, tem a representação do vento, com uma velocidade V que faz mover a turbina com uma rotação ômega. Figura 40 - Representação da velocidade V e rotação W do vento. Para poder entender o motivo que faz com que a pá rode quando o vento passa por ela, é preciso analisar o corte da pá em uma vista superior, como representado na Figura 41. 47 Figura 41 - Pressão e velocidade. Quando o vento se aproxima da pá, uma parte vai para um lado da pá e a outra parte para a outra parte. Ou seja, o ar se divide. Dessa maneira, ele gera na parte superior, uma região de velocidade menor e pressão maior. E do outro, velocidade maior e pressão menor. Como uma pressão de um lado é maior que a pressão do outro lado, a força de um lado é maior que a força de outro lado. Logo, existirá uma força resultante Fr, como pode ser observada na Figura 42, que vai empurrar a pá. Figura 42 - Força resultante. Nesta representação, a pá gira para baixo, o vento vai da esquerda para a direta, só que após passar pela pá, o vento sai em ângulo, como pode ser observado na Figura 43. 48 Figura 43 - Ângulo formado. O vento sai com um determinado ângulo, mas a velocidade angular é em sentido diferente da velocidade do vento V. O vento vai “bater” na pá em diversas localizações. Inclusive em algumas áreas ele irá passar reto, com um pouco de turbulência que pode inclusive ser oriunda de passagens de outras pás que influenciaram. Na parte superior da representação, o vento que sai do ângulo encontra o vento que entra na outra direção. A velocidade do vento que sai da pá, ao encontrar com a outra velocidade ocorre a mudança de direção. Além disso, o encontro intensifica ainda mais o fenômeno da turbulência, como pode ser observado na representação da Figura 44. Figura 44 - Intensificação da turbulência. Logo, o estudo fez com que fosse concluir que passagem do escoamento pela pá, gera uma diferença de pressão, ocasionando uma força resultante que faz com que a pá se desloca. Isso interfere na direção que será lida no anemômetro. 2.14 Simulações Em uma simulação inicial foi possível perceber que havia muitos erros na configuração. Foi preciso, portanto, realizar simulações menores, com cerca de 10 minutos, para avaliar pouco 49 a pouco os resultados e, assim, realizar os ajustes necessários. Para isso, foi utilizado uma malha com domínio maior, dado que demanda uma carga computacional menor. Além disso, foi retirado os custom controls para não reduzirem muito o tamanho da malha. Foi acordado que inicialmente iria ser descartado a primeira rotação completa da turbina para evitar usar dados que ainda não convergiram adequadamente. Assim, a análise deve ser restrita às duas voltas seguintes. Entretanto, no decorrer do projeto até o fim, foi acordado em descartaríamos as 10 primeiras voltas. 2.15 Design of Experiments (DoE) O primeiro passo foi realizar uma análise para verificar quais ferramentas existentes poderiam suprir a necessidade do grupo e gerar o estudo. Dessa forma o grupo encontrou algumas bibliotecas em Python que fazem a geração dos designs. Dentre elas a PyDoE, DoEgen, pyDoE2 e DoEpy. Além disso foram encontradas ferramentas em R que também geravam o estudo. Antes da escolha da ferramenta foi importante observar os princípios de DoE. Um dos princípios do DoE é a randomização. A randomização é usada para garantir a aleatoriedade do estudo e tenta garantir que que toda distribuição possível de tratamentos tenha a mesma probabilidade de ser escolhida. Juntamente com replicação, outro princípio de DoE que garante a replicação dos diferentes leveis no resultado do estudo, formam os pilares do DoE em conjunto com Local Control. Local Control nada mais é que uma maneira de controlar os leveis a fim de balancear o resultado do estudo aumentando assim a eficiência do estudo. Seguindo o primeiro princípio do DoE foi elaborado um código em Python que dentre um Full Factorial selecionava durante dez milhões de iterações qual dos experimentos mais foi escolhido pela aleatorização. Isso gerou um resultado não tão balanceado que foi descartado Depois disso, houve necessidade de utilizar um Fractional Factorial. Entretanto, o que temos conhecimento, é que a maioria dos estudos de DoE geram experimentos a partir de um range para cada um dos fatores ou leva em consideração dois leveis somente, o que inviabilizou a utilização inicial das ferramentas. Dessa forma foi feito um estudo da ferramenta pyDoE2 e verificou-se a possibilidade de utilizar Generalized Subset Designs que é uma generalização de Fractional Factorial que permite a entrada de um ou mais fatores com diferentes níveis. De tal maneira que foi possível gerar o seguinte estudo presente na Figura 45. 50 Figura 45 - Estudo aplicado. É possível verificar certo controle nas variáveis observando um balanceamento dos leveis. Além disso pode se verificar a replicação dos leveis. Para randomizar o estudo fpo executado os experimentos de forma aleatória. Isso tudo lembrando que foi fixado valores específicos de velocidade do vento no hub, rotação e intensidades turbulentas. Para garantir a eficácia do estudo foi realizado o cálculo do A-efficiency que um medidor de eficácia de DoE. A fórmula para calcular está abaixo: Dessa fórmula o p é número de fatores, Nd é o número de experimentos e X é a matriz normalizada da saída do GSD. O resultado obtido foi de noventa e dois ponto trinta levando em consideração cem pontos possíveis. 51 Para que fossem realizadas todas as simulações do DoE, foi necessário criar nove arquivos base que englobariam todas as combinações de direção do vento e inflow requisitadas pela WEG, no caso -15°, 0° e 15° de direção do vento no plano XY e -8°, 0° e 8° de inflow. Este processo considerou um direcionamento de entrada do vento diferente, e por isso, foi necessário alterar o posicionamento de todo o domínio computacional. Como a face frontal do cilindro foi definida como velocity inlet, ou seja, o local de entrada do vento no sistema, foi preciso que ela estivesse sempre perpendicular à direção do vento, caso contrário a lateral do cilindro também estaria atuando como entrada do vento e isto não estaria correto, portanto, caso o vento possua 15° no plano XY era também preciso rotacionar o domínio computacional no mesmo valor, como evidenciado na Figura 46. Figura 46 - Domínio computacional em vista superior. Apenas o domínio cilíndrico que rotaciona, o resto deve manter-se imóvel, até porque não haveria diferença de angulação caso a turbina também rotacionasse. Fora isso, para toda mudança de orientação do vento e consequentemente do domínio computacional foi necessário refazer a região, alterar seus parâmetros, refazer as interfaces e a malha. Para alterar a direção do vento dentro do software foi preciso descobrir qual deveria ser seu valor vetorial, já que o sistema de coordenadas está dividido em XYZ foi preciso descobrir o valor que cada uma dessas componentes deveria ter para que o vento. Até o momento todas as simulações foram realizadas com 0° de angulação no plano XY e 0° de inflow, logo, todo o 52 direcionamento do vento estava sendo aplicado no eixo X que é o mesmo do eixo da linha de centro da nacele. Para isto foi feito o cálculo vetorial para todos os casos seguindo a seguinte lógica: A partir destas equações foi possível perceber por exemplo que o vento que possui -15° de direção horizontal, -8° de inflow e 10 m/s de magnitude possui na verdade 9.565255025 em X, -2.563002359 em Y e -1.39173101 em Z e deveria ser implementado no software da seguinte forma, como evidenciado na Figura 47: Figura 47 - Propriedades aplicadas. Com o auxílio do excel foi criada uma planilha que informava todos os valores de cada componente da velocidade para cada uma das 39 simulações do DoE. 53 2.16 Matriz de correlação Uma etapa importante do projeto foi a partir da criação e análise da matriz de correlação, como pode ser observado na Tabela 8 e na Tabela 8, para mostrar a relação presente entre as variáveis utilizadas no escopo do projeto. Tabela 7 - Análise da matriz de correlação entre as variáveis de entrada das simulações e a direção do vento da leitura no anemômetro A partir da Tabela 7, é possível observar uma alta correlação entre a direção do vento na entrada do sistema e a sua leitura do anemômetro. Para a análise, foi retirado o inflow por considerar que a direção do vento foi somente calculada com as componentes da direção em x e y, excluindo a componente em z. 54 Tabela 8 - Análise da matriz de correlação entre as variáveis de entrada das simulações e a velocidade do vento lida no anemômetro A partir da Tabela 8, é possível observar uma média correlação entre a velocidade de entrada no hub e a velocidade na leitura do anemômetro. Para a análise, também foi retirado o inflow por considerar que a direção do vento foi somente calculada com as componentes da direção em x e y, excluindo a componente em z. 3. Resultados 3.1 Simplificações do modelo CAD O modelo enviado pela empresa apresenta algumas formas e informações que não serão utilizadas para a simulação de CFD. Dessa forma foram realizadas simplificações no modelo utilizando o software STAR-CCM+. Primeiramente foi realizada a importação do arquivo completo do CAD, que não foi acompanhado do hub. Para isso foi preciso utilizar a opção do STAR-CCM+ "Import CAD model into 3D-CAD", como pode ser visualizado na 8. 55 Figura 48 – Importação do arquivo CAD. Uma vez feito isso iniciou-se o processo de remoções. Foi feita a remoção de parafusos que ligam as pás da turbina com o hub (peça central da turbina onde são ligadas as pás). Pode- se verificar na Figura 49 que os parafusos foram selecionados e removidos em segurança para preservar o resto do modelo CAD. Figura 49 - Remoção dos parafusos externos. 56 É possível observar nas Figura 50, Figura 51 e Figura 52. Figura 50 - Depois da remoção dos parafusos externos. Após isso foram removidos os parafusos internos, as pás que ligavam a junção das pás com o hub de forma interna. Figura 51 - Remoção dos parafusos interno. Além dos parafusos, foram removidos os espaços onde eles ficavam, com o mesmo método dos parafusos, selecionar e deletar. A Figura 52 ilustra a remoção dos espaços dos parafusos das pás. 57 Figura 52 - Remoção dos espaços dos parafusos externos. Foi verificado que o hub era uma estrutura com furos que possibilitavam a geração da malha em seu interior. Foi removido o hub, mantido os conectores das pás e realizada uma extensão dos conectores até que se encontrassem no centro do rotor, como é possível verificar na Figura 53: Figura 53 - Extensão dos conectores das pás. Feito isso, foi possível verificar que existia sobreposição de uma das pás que foi necessário apagar, como evidenciado na Figura 54. 58 Figura 54 - Remoção da sobreposição de uma das pás. Assim como os espaços dos parafusos tiveram que ser apagados os lugares onde ficavam esses parafusos também precisou ser deletado. Com a Figura 55 é possível verificar essa remoção. Figura 55 - Remoção do espaço onde ficavam os parafusos. Teve início as simplificações do cooler. Inicialmente, foram apagadas as mãos francesas que ligam o cooler à nacele, além de apagar os suportes do cooler, porcas e parafusos de acordo com a Figura 56. 59 Figura 56 - Simplificações do cooler. Ao terminar as simplificações continuaram na parte de trás da turbina no anemômetro. O modelo CAD veio com uma serie de porcas e parafusos no anemômetro e uma parte de trás que precisou ser removida. Além disso em torno do anemômetro tinha um raio que foi necessário remoção também. Tudo isso com o mesmo método dos parafusos. A Figura 57 evidencia o que foi removido do anemômetro. Figura 57 - Simplificações do anemômetro. A Figura 58 mostra com mais detalhe algumas das remoções do anemômetro. 60 Figura 58 - Mais remoções realizadas no anemômetro. Para finalização foi feita uma duplicata da haste do anemômetro, como evidenciado na Figura 59 e, utilizada como base e topo. Além disso foi necessário alterar a altura da haste até que ficasse parecida com um disco. Alterando a escala e movendo essa haste modificada até a última base e último topo existente, foi possível continuar com as simplificações. Na Figura 60 é possível verificar o movimento da haste modificada até a base de um dos anemômetros. 61 Figura 59 - Anemômetros modificados. Figura 60 - Movimento da haste modificada. Após isso volta-se às pás. Para consertar uma das pás foi necessário utilizar o método Fill Holes como mostrado na Figura 61. Foi utilizado o mesmo método para as outras duas pás. 62 Figura 61 - Remoção de bordas redundantes. Na nacele foi necessário remover as estruturas de acesso e simplificar retificando a estrutura superior. A Figura 62 e Figura 63 mostram como era e como ficou a nacele, respectivamente. Figura 62 - Nacele sem alterações. 63 Figura 63 - Nacele com alterações. Próximo passo é o corte das pás. Para realizar o corte das pás foram feitos desenhos de retângulos para utilizar o método "Create Revolve Cut" para cada um dos desenhos com o ângulo de 360 graus. Os desenhos começaram da parte de cima das pás e seguindo até o raio de 16 metros das pás. Dessa forma foi possível cortar gradualmente as pás da turbina. Feito isso foi possível começar as simulações. 3.2 Simulações A realização das simulações foi feita baseando-se no vídeo [3] de onde o grupo pode retirar boa parte dos parâmetros para colocar no software STAR-CCM+. Os parâmetros utilizados podem ser encontrados no apêndice A. Após a realização das simulações foi possível verificar velocidades do vento na parte de trás do rotor. Essas velocidades, que podem ser observadas na Figura 64, com o registro dos anemômetros marcando uma velocidade próxima à 8.67 m/s. 64 Figura 64 - Velocidade do vento atrás do rotor. A Figura 65 registra a velocidade do vento na parte da frente do rotor depois das simulações feitas com o software STAR-CCM+. É possível verificar uma velocidade mais alta perto das pás. Figura 65 - Velocidade do vento na frente do rotor. 65 3.3 Validação da malha após a independência de malha A partir do gráfico representado na Figura 66, é possível ver o número de células com um ângulo de assimetria específico para todos os limites no domínio computacional. Figura 66 - Frequency versus skewness angle Para visualizar o ângulo de distorção máximo encontrado no domínio computacional, foi aplicado um Zoom, como pode visualizado na Figura 67. Figura 67 - Zoom de Frequency versus skewness angle. O ângulo de distorção máximo é menor que 85 graus, o que é aceitável para uma simulação de CFD. Abaixo há o passo a passo desenvolvido: 66 1) Métrica de validade de face A métrica de validade de face é uma medida de correção da face normal em relação ao centroide da célula à qual elas estão anexadas. A partir da Figura 68, uma célula de boa qualidade é representada quando as normais estão apontando para fora em relação ao centroide da célula, enquanto em uma célula de validade de face ruim, uma ou mais normais estão apontando para o centroide da célula. Figura 68 - Qualidade da célula. Uma validade de face de 1 é uma indicação de uma célula de boa qualidade, no entanto, os valores menores que 0,5 indicam a malha de volume negativo e devem ser evitados. A partir da Figura 69 e na Figura 70 é possível observar a aplicação da validade da face para o projeto. Figura 69 - Validade da face. 67 Figura 70 - Zoom da validade da face. Todas as células no domínio computacional têm uma validade de face próximo a 1, o que mostra que as normais estão apontando corretamente para fora do centroide da célula. 2) Métrica de qualidade de célula Na métrica de qualidade de célula, geralmente, células planas com faces altamente não ortogonais têm uma qualidade de célula baixa. Para exemplificar a análise há a Figura 71. Figura 71 - Qualidade da malha. Uma célula com qualidade 1 é considerada boa, como ocorre com a célula cúbica. No entanto, uma célula com baixa qualidade tem uma qualidade de célula próxima de zero. Células com qualidade inferior a 10−5 são consideradas ruins. A Figura 72 mostra a qualidade da célula usando os dados do projeto. 68 Figura 72 - Qualidade da célula. Para visualizar a qualidade mínima da célula encontrada no domínio computacional, basta aplicar um Zoom como pode ser observado na Figura 73. Figura 73 - Zoom da qualidade da célula. É possível observar que a qualidade mínima da célula é de cerca de 0,52, o que é aceitável para uma simulação de CFD. 3) Métrica de alteração de volume A métrica de alteração de volume avalia a proporção do volume de uma célula para o maior volume de suas células vizinhas. Esses números mostram se a mudança de volume é boa ou má. As variações acentuadas dos volumes das células nas proximidades de uma célula podem causar sérias imprecisões e instabilidades. As células com baixa qualidade são aquelas com uma alteração de volume de 0,01 ou menos. A análise do projeto está representada na Figura 74. 69 Figura 74 - Alteração de volume. Para visualizar a mudança mínima de volume encontrada no domínio computacional, é preciso ampliar, como mostra a Figura 75. Figura 75 - Zoom da mudança do volume. A mudança mínima de volume é de cerca de 0,04, o que é aceitável para uma simulação de CFD. 4) Indicador de qualidade Chevron O indicador de qualidade Chevron associa um valor de 1 às células nas quais a linha que une os centros das células não passa pela face comum. Tais células são consideradas ruins. O valor zero é atribuído a todas as outras células. Por exemplo, a Figura 76 mostra uma célula normal e uma célula Chevron. 70 Figura 76 - Célula Chevron e célula normal. As células Chevron afetam a precisão e a robustez no cálculo dos fluxos convectivos. A análise da qualidade Chevron do projeto pode ser observada na Figura 77. Figura 77 - Indicador de qualidade Chrevron. Todas as células no domínio computacional têm um indicador de qualidade Chevron igual a 0, o que mostra que não há células Chevron no domínio computacional. 5) Localização de células de má qualidade A Figura 78 representa em destaque as células de má qualidade. Figura 78 - Localização de células de má qualidade. 71 A Figura 79 representa um Zoom para uma melhor verificação. Figura 79 - Zoom para a localização de células de má qualidade. É possível ver as células com os valores do ângulo de assimetria entre 30 e 85. Na Janela de propriedades, como evidenciado na Figura 80, é possível alterar o intervalo do ângulo de assimetria. Figura 80 - Janela de propriedades. 72 6) Remoção de células inválidas O STAR-CCM+ tem a capacidade de identificar e remover as células de baixa qualidade após o processo de geração de malha. Até 6 critérios de malha diferentes podem ser usados ao escolher quais células serão removidas. Os critérios são: 1. Validade mínima de face: minimum face validity. 2. Qualidade mínima da célula: minimum cell validity. 3. Alteração mínima de volume: minimum volume validity. 4. Mínimo de células contíguas: minimum contiguous cells. 5. Área mínima de face conectada: minimum connected face área. 6. Células de volume negativo: negative volume cells. As células removidas são movidas para uma região separada e não são usadas na análise CFD. O STAR-CCM+, adicionando limites de plano de simetria às faces de células vizinhas de células que foram removidas, minimiza o efeito dessas áreas removidas na solução. Aplicando na análise do projeto, a Figura 81 evidencia tais critérios, os valores em “Reset” e a aba “Identify” mostra a qualidade de células. Definindo a qualidade da célula para 0,03, foi possível observar que há 25 células, cuja métrica de qualidade de célula é menor do que 0,03. Figura 81 - Critérios da malha. Se todos os valores na tabela ficar em “valores padrão”, que é quando aciona o “Reset”, nenhuma célula invalida poderá ser identificada, o que é uma indicação de qualidade razoável da malha, como pode ser observado na Figura 82. 73 Figura 82 - Valores padrão. A partir dessa análise, vale destacar que o melhor método para cuidar de células inválidas ou, em geral, de qualquer problema de qualidade de malha é primeiro localizar essas células no domínio computacional e, em seguida, aplicar refinamentos de malha de superfície ou volume adequadas a essas regiões. 3.2 CFL Foi realizado um teste, a partir de simulações, com CFL de 50.0, 1.0, 0.7 e 0.5 para verificar qual o melhor que se adequava ao experimento. O resultado está presente na Figura 83. Figura 83 - Simulações com CFLs diferentes. 74 A partir da Figura 83, é possível concluir que com o CFL 0,5 demora mais tempo para estabilizar e a curva se mostra mais distinta do que as demais. As curvas 1.0 e 0.7, tem efeitos opostos, embora apresente curvas com formatos parecido. A que mais se mostrou adequada com a nossa análise, é a curva com o CFL 1. No gráfico o efeito de direção negativa do vento aparece anteriormente à passagem da pá pelo anemômetro. Algumas vezes ainda fica com um direcionamento positivo. 3.3 Modelos de turbulência Foi realizado simulações para verificar qual o modelo de turbulência mais adequado para ser aplicado, como pode ser evidenciado pela Figura 84 e pela Figura 85. Figura 84 - Simulação com diferentes modelos de turbulência. 75 Figura 85 - Simulação com diferentes modelos de turbulência. A partir das curvas dos três modelos de turbulência, foi possível afirmar que realmente o realizable-Kɛ (curva verde) é o que melhor se comporta em relação ao esperado. Embora o RSM (curva azul) possua um comportamento repetitivo, o que faz sentido a longo prazo, ele ainda não é tão estável quanto o realizable-Kɛ. Já o SST- kω (curva amarela) foi o que se comportou pior, pois seus vales possuem uma amplitude muita elevada tanto na direção quanto na velocidade e picos também fora do comum em sua velocidade, seu comportamento não parece estável já que sua repetição não se observa, a curva parece se comportar diferentemente a cada rotação da turbina. ´É possível verificar, a partir das Figuras 86-88 os resíduos dos modelos de turbulências analisados. 76 Figura 86 - Residuals RSM. Figura 87 – Residuals SST- kω. Figura 88 - Residuals Realizable-Kɛ 77 Pelos gráficos dos resíduos também é possível dizer que o realizable-Kɛ é o melhor entre os três modelos. Os resíduos do RSM possuem um comportamento bem diferente em relação aos outros dois, seu valor de continuidade cresceu conforme o andamento da simulação, o que não é nada bom visto que ele deveria diminuir conforme o progresso da simulação e seu valor está em aproximadamente 10. O comportamento do realizable-Kɛ e do SST- kω são semelhantes, mas ao olhar para o eixo vertical do gráfico percebeu-se que o realizable-Kɛ possui uma ordem de grandeza abaixo, portanto seu resultado é o melhor dos três e consequente ele continuou sendo utilizado nas simulações. 3.4 Resultado das simulações do DoE A partir dos resultados obtidos pelo DOE foi possível realizar várias análises, a primeira foi sobre a influência que a velocidade de entrada do vento possui na leitura do anemômetro, como pode ser observada na Figura 89. Para tal foram utilizadas as simulações WEGDOE10 (6 m/s), WEGDOE23 (8 m/s) e WEGDOE36 (10 m/s). Figura 89 - Análise da influência das velocidades de entrada. Verificou-se que para a velocidade de entrada mais baixa as curvas de velocidade e direção média do anemômetro apresentam um comportamento muito diferente das outras. A curva azul (6 m/s) possui picos e vales muito mais acentuados. Isso ocorre pelo fato de que para 78 velocidades e rotação do rotor menores, a influência da passagem das pás é muito maior. Fica marcado claramente na curva o momento da passagem da pá em frente ao anemômetro. A segunda análise foi referente à influência da direção horizontal de entrada do vento na simulação, como pode ser observada na Figura 90. Foram utilizadas as seguintes simulações: WEGDOE19, WEGDOE23 e WEGDOE40 que foi uma simulação feita exclusivamente para esta análise. As três possuem todos os parâmetros iguais com exceção da direção do vento, que são 15°, 0° e -15° respectivamente. Veja suas curvas na Figura 90. Figura 90 - Análise da Direção horizontal do vento É notável que a direção influencia a leitura do anemômetro, tanto na amplitude, quanto no formato da curva. Foram apresentadas as curvas nas últimas 5 voltas, que totalizam os 1800°, mas para analisar com mais precisão a influência da passagem da pá em frente ao anemômetro foi feita uma média dos 1800° para 120°, que se refere à passagem de uma pá, , como pode ser observada na Figura 91. 79 Figura 91 - Análise da influência da passagem de uma pá em frente do anemômetro É possível perceber que agora a diferença no formato da curva ficou muito mais sutil do que antes, mas a diferença de amplitude é inquestionável. Foi calculada a média de cada uma das curvas e o erro de direção, ou seja, a diferença do valor de direção de entrada do vento com a direção lida pelo anemômetro, em módulo, como pode ser observada nas Figuras 92 e 93. Figura 92 - Média das curvas da análise da direção do vento 80 Figura 93 - Erro da direção do vento lida pelo anemômetro É notável que existe uma correlação entre o erro da leitura e o valor de direção de entrada do vento. Quando o vento está inclinado a favor da direção de rotação da turbina, ou seja, da esquerda para a direita, sentido horário, o valor de erro é maior, no caso aproximadamente 4,5°. Quando o vento está contra a direção de rotação da turbina o erro é bem menor, aproximadamente 1,7. Foi feita uma análise para ver o comportamento da curva quando se varia apenas o inflow, tanto negativamente quanto positivamente. Para tal utilizou-se as seguintes simulações: WEGDOE5, WEGDOE7 e WEGDOE10 para a velocidade de 6 m/s; WEGDOE18, WEGDOE20 e WEGDOE23 para a velocidade de 8 m/s; WEGDOE31, WEGDOE33 e WEGDOE36 para a velocidade de 10 m/s, como evidenciada nas Figura 94-99. 81 Figura 94 - Análise da influência do inflow para velocidade de 6 m/s. Figura 95 - Análise da influência do inflow na direção do vento para a velocidade de 6 m/s. 82 Figura 96 - Análise da influência do inflow para velocidade de 8 m/s. Figura 97 - Análise da influência do inflow na direção do vento para a velocidade de 8 m/s. 83 Figura 98 - Análise da influência do inflow para velocidade de 10 m/s. Figura 99 – Análise da influência do inflow na direção do vento para a velocidade de 8m/s. É possível observar pelos gráficos que há diferenças entre variar o inflow negativamente e positivamente, mas que o comportamento das curvas respectivas não é oposto como poderia ser sugerido. Verifica-se também que para a mesma direção do vento, o comportamento da curva nem sempre é o mesmo, as curvas WEGDOE7, WEGDOE20 e WEGDOE33, cujo inflow é positivo, são praticamente idênticas, enquanto as curvas WEGDOE5, WEGDOE18 e WEDOE31, de inflow negativo, apresentam o comportamento assimilar, levantando um questionamento sobre a qualidade das simulações. 84 Dessa forma, o vento incidir sobre a pá de “cima para baixo”, ou seja, inflow positivo, conforme Figura 100Figura 100, demonstrou ser muito mais estável já que as curvas possuem o comportamento similar. Figura 100 - Vista lateral do domínio computacional com inflow positivo. Uma outra análise foi feita variando a direção do vento, mantendo um inflow inalterado. Para essa comparação foram utilizadas as simulações WEGDOE23, WEGDOE24 e WEGDOE26 com valores de direção do vento de -15°, 0° e 15 respectivamente. Todas com uma velocidade de 8 m/s, e inflow 0°, de acordo com as Figuras 101-104. 85 Figura 101 - Análise da influência da velocidade de 8 m/s variando a velocidade do vento. Figura 102- Análise da influência da velocidade de 8 m/s variando a direção do vento. 86 Figura 103 - Análise da influência da velocidade de 10 m/s variando a direção do vento Figura 104 - Análise da influência da velocidade de 10 m/s variando a direção do vento. Quando se analisa a velocidade das três simulações, na Figura 104, percebe-se que todas possuem curvas semelhantes, variando apenas um pouco da magnitude, tendo a WEGDOE26 em alguns pontos com uma velocidade igual, e até um pouco mais que 8 m/s, indicando um possível mal comportamento do resultado da simulação. A simulação em questão é a que possui uma direção de entrada do vento de 15 °. 87 Já observando a direção do vento, na Figura 102Figura 102, nota-se que a WEGDOE24, e a WEGDOE26 são as simulações com maior amplitude na direção do vento. Podendo indicar um pior resultado da simulação, quando se aplica uma direção de entrada do vento de -15° e 15°, em comparação com a WEGDOE23, que com direção de entrada de 0°, se manteve mais uniforme, e com menor amplitude da direção. Na análise da Figura 103 das velocidades médias das simulações WEGDOE36 (0°), WEGDOE37 (-15°) e WEGDOE39 (15°), é possível observar que as curvas possuem similaridade em seu comportamento. Além disso, é possível verificar que a diferença da direção do vento na entrada gera diferentes velocidades no anemômetro mesmo mantendo uma magnitude de vento na entrada igual. Se a angulação de entrada do vento é positiva a velocidade fica maior do que quando ela é nula ou negativa. Analisando a Figura 104 percebe-se que a variação da direção do vento na entrada do rotor é significante para causar diferenças na leitura da direção do vento no anemômetro. Quando a angulação do vento na entrada do rotor é negativa (curva amarela) nota-se que a direção no anemômetro não ultrapassa 5°. A angulação do vento positiva na entrada do rotor (curva verde) faz com a leitura atinja cerca de 25°, ou seja, o ângulo de entrada horizontal do vento no sistema tem uma grande influência na amplitude da curva. Observa-se também a influência da passagem das pás na frente do anemômetro que causam vales e picos na leitura da direção do vento no anemômetro e sua uniformidade de ocorrência nas três simulações. Para uma análise mais geral que possa ser aplicada ao controlador da turbina eólica da WEG foram feitas todas as 39 simulações e extraídos todos os erros tanto de direção como de velocidade como é possível observar na Figura 105 (erros da direção em graus) e Figura 106 (erros da velocidade em m/s). Figura 105 - Erros da direção em graus. 88 Figura 106 - Erro na leitura do anemômetro x Simulações do DoE. Tanto na velocidade como na direção foi possível observar uma variação de erros próxima de 0 até valores mais altos com uma média nos dois casos próxima de 2 (graus e m/s). Dessa forma é possível verificar que foi possível obter uma superfície de resposta para a calibração da leitura do anemômetro o que potencialmente pode gerar um ajuste no controlador implicando numa maior geração de energia. 3.5 Efeito da passagem da pá em frente ao anemômetro. Além de gráficos e curvas, o grupo analisou visualmente a influência da passagem da pá em frente ao anemômetro. Selecionando um intervalo de tempo referente a 120° que equivale a passagem de uma pá, foram tirados screenshots de todos os times step para que se gerasse um GIF e assim ficasse mais visual o efeito gerado pela pá no vento que a percorre, como pode ser observado na Figura 107. Figura 107 - GIF da simulação A Figura 107 é referente à simulação WEGDOE23, que possui direção horizontal e inflow do vento de entrada igual a zero. Verifica-se que a pá não só altera a direção, mas também a magnitude do vento, logo são plausíveis os resultados apresentados nas curvas do DoE. Nota- se que a região inferior da pá possui alta velocidade e consequentemente baixa pressão, que é 89 o oposto da região superior, o que geraria uma força resultante no sentido de rotação da pá, assim como ocorre naturalmente. Logo, mesmo com o movimento da turbina sendo prescrito, ou seja, não dependente desta força do vento, ainda assim o fenômeno natural é respeitado, isso faz com que possua mais argumentos para defender que a simulação mesmo com suas limitações e simplificações, representa de forma coerente e satisfatória o escoamento do vento pela turbina. 4. Conclusão e trabalhos futuros Este trabalho demonstrou que CFD é uma área de estudo muito extensa de investigação. O objetivo deste projeto foi concluído com sucesso, gerou-se uma superfície de resposta razoavelmente ampla que pode ser utilizada como parâmetro de comparação para a calibração de posicionamento do anemômetro. O estudo desenvolvido durante este semestre possui muito potencial de expansão e aprimoramento. O DOE foi uma forma de otimização da pesquisa devido à limitação de tempo para o projeto. Entretanto limitou a amplitude das análises realizadas, já que para analisar a influência real da mudança de um parâmetro é preciso que apenas ele seja alterado enquanto os outros se mantêm inalterados. A redução do número de simulações inviabilizou isso, mas este fator não invalida os resultados apresentados, apenas indica um ponto de melhora. Portanto, o grupo define este projeto como um ponto de partida para uma investigação mais ampla que pode ser feita pelos engenheiros da WEG ou por outros alunos em uma extensão deste projeto em um outro projeto final de engenharia do Insper. Merece destaque também a metodologia de validação e independência da malha desenvolvida pela equipe. Seus resultados enriqueceram muito a qualidade do estudo já que possibilitaram afirmar com mais propriedade que as simulações estavam adequadas ao cenário real. Mais do que fornecer os dados de calibração do anemômetro, foi elaborada uma documentação precisa de como realizar satisfatoriamente um estudo de simulação CFD que seja condizente com o escoamento natural. 90 5. Referência bibliográfica [1] Isto é WEG: Nosso propósito. WEG headrquarters, 2022. Disponível em: https://www.weg.net/institutional/BR/pt/this-is-weg. Acesso em: 07 mar. 2022. [2] BUHLER, Alexandre Jose. ResearchGate, 2014. Disponível em: https://www.researchgate.net/figure/Figura-2-Perfil-vertical-da-velocidade-do-vento-desde-a- superficie-ate-a-altura-do-vento_fig2_292623126. Acesso em: 07 mar. 2022. [3] Tutorial: CFD simulation of a Wind Turbine. Disponível em: https://www.youtube.com/watch?v=KnfY9xElKYY&t=1312s. Acesso em: 07 mar. 2022. [4] Zhu, Xiaocheng. Sun, Chong. Ouyang, Hua. Du, Zhaohui. Numerical investigation of the effect of towers and nacelles on thenear wake of a horizontal-axis wind turbine model. Publicado em: 10 de ago. 2021. Local: China. [5] Hamlaouia, M.N. Smailia, A. Dobrevb, I. Pereira, M. Fellouahc, H. Khelladi, S. Numerical and experimental investigations of HAWT near wake predictions using Particle. Publicado em: 29 jul. 2021. Local: Algéria, França e Canada. [6] Kragh, Knud A. Hansen, Morten H. Potential of power gain with improved yaw alignment. Publicado em: 18 fev. 2014. Local: Dinamarca. [7] O’Brien, J.M., et al. “An Assessment of Commercial CFD Turbulence Models for near Wake HAWT Modelling.” Journal of Wind Engineering and Industrial Aerodynamics, vol. 176, May 2018, pp. 32–53, 10.1016/j.jweia.2018.03.001. Publicado em: 12 mai. 2020. Local: Irlanda do Norte. [8] Sedaghatizadeh, Nima, et al. “Modelling of Wind Turbine Wake Using Large Eddy Simulation.” Renewable Energy, vol. 115, Jan. 2018, pp. 1166–1176, 10.1016/j.renene.2017.09.017. Local: Australia. [9] Regodeseves, P. García, and C. Santolaria Morros. “Unsteady Numerical Investigation of the Full Geometry of a Horizontal Axis Wind Turbine: Flow through the Rotor and Wake.” Energy, vol. 202, July 2020, p. 117674, 10.1016/j.energy.2020.117674. Accessed 17 Sept. 2020. Local: Espanha. 91 [10] Nilay Sezer-Uzol and Lyle N. Long. “3-D Time-Accurate CFD Simulations of Wind Turbine Rotor Flow Fields” The Pennsylvania State University, University Park, Local: Estados Unidos. [11] IDEAL SIMULATIONS. Courant number in CFD simulations. Disponível em: https://www.idealsimulations.com/resources/courant-number-cfd/. Acesso em: 28 mar. 2022. [12] Planejamento de Experimentos DOE. Disponível em: http://www5.eesc.usp.br/portaldeconhecimentos/index.php/por/Conteudo/Planejamento-de- Experimentos-DOE. Acesso em: 24 de mar. 2022. [13] INTRODUÇÃO A CONVECÇÃO. Disponível em: http://ftp.demec.ufpr.br/disciplinas/TMEC030/Prof_Stephan/01_introducaoConvecao.pdf. Acesso em: 22 de mar. 2022. [14] APPLIED COMPUTATIONAL FLUID DYNAMICS. Tutorial: Checking the Validity and Quality of Mesh (STAR-CCM+). Disponível em: https://www.youtube.com/watch?v=ALGnQFLlrkg&t=121s. Acesso em: 25 fev. 2022. [15] Anônimo. "SST k-omega model". Disponível em: https://www.cfd- online.com/Wiki/SST_k-omega_model. Acesso em: 6 de junho de 2022. [16] Anônimo. "Reynolds stress model (RSM)". Disponível em: https://www.cfd- online.com/Wiki/Reynolds_stress_model_(RSM). Acesso em: 6 de junho de 2022. [17] Anônimo. "Realisable k-epsilon model". Disponível em: https://www.cfd- online.com/Wiki/Realisable_k-epsilon_model. Acesso em: 6 de junho de 2022. [18] Surowiec, Izabella et al. "Generalized Subset Designs in Analytical Chemistry". Disponível em: https://pubs.acs.org/doi/pdf/10.1021/acs.analchem.7b00506. Acesso em: 6 de junho de 2022. [19] Clicumu (github user). "pyDOE2". Disponível em: https://github.com/clicumu/pyDOE2. Acesso em: 6 de junho de 2022. 92 Apêndice A Parâmetro Valor Definition of the Computation Domain To create the cylindrical rotating subdomain, For the start coordinates put where the rotor will be located. [0.0,0.3,0.0]m,m,m whole for the end coordinates put [0.0,-1.25725,0.0]m,m,m. Put 7.5435 m for the radius of this rotating subdomain. Note that this value is equivalent to 1.5 times the radius of the turbine rotor. Outer stationary domain: Conjunto de corner 1: [-18,3, -20,116, -12,9] m,m,m coordenadas para o canto 1 e canto 2 corner 2: [-18,3, 150,87, 11,5] m,m,m Definition of the Regions First all the surfaces are created (hub, left Surfaces blade, left blade_base,right blade, right blade_base) Subtract rotating subdomain from the wind Boolean - Subtract - Wind Tunnel Tunnel Subtract TowerNacelle Boolean - Subtract - TargetPart(Subtract) Subtract rotor from rotating subdomain Rotor - CilindricalRotatingSubdomain-Boolean Add the Rotating Subdomain and Virtual Drag wind tunnel to geometry 1 Assign parts to regions Virtual RotatingSubdomain - Assign parts to regions - Create a region for each part - Create a boundary for each part surface - Do not create interfaces from contacts Mesh Generation Define the interface between the two VirtualWindTunnel - Boundaries - domains RotatingPart1 - Create interface Create Mesh Continua - New - Mesh Continuum Assign Mesh1 to VirtualRotatingSubdomain Right Click on models - Models - Meshing Models - Surface Remesher - Trimmer - Prism Layer Mesher Mesh Alignment Trimmer - Do mesh Alignment Reference Value Reference Value - Maximum Cell Size - 100 Reference Value Reference Value - Number of Prism Layers - 15 93 Reference Value Reference Value - Prism Layer Stretching - 1.1 Mesh VirtualWindTunnel -> Reference 400 Values -> Maximum cell size -> Percentage of base Mesh VirtualWindTunnel -> Reference 0,55 Values -> Surface size -> Relative Minimum Size Mesh VirtualWindTunnel -> Reference 0,55 Values -> Surface size -> Relative Target Size Mesh VirtualWindTunnel -> Reference Slow Values -> Template Growth Rate -> Default Growth Rate Mesh VirtualWindTunnel -> Reference Slow Values -> Template Growth Rate -> Boundary Growth Rate Regions -> VirtualWindTunnel -> Mesh Mesh VirtualWindTunnel Continuum Regions -> VirtualWindTunnel -> Check Boundaries -> Select Bottom, Inlet, Outlet, Sides, Top -> Right click (edit) -> Mesh conditions -> Custom surface size Regions -> VirtualWindTunnel -> 400 Boundaries -> Select Bottom, Inlet, Outlet, Sides, Top -> Right click (edit) -> Mesh conditions -> Custom surface size -> Mesh Values -> Surface Size -> Relative Minimum Size Regions -> VirtualWindTunnel -> 400 Boundaries -> Select Bottom, Inlet, Outlet, Sides, Top -> Right click (edit) -> Mesh conditions -> Custom surface size -> Mesh Values -> Surface Size -> Relative Target Size Regions -> VirtualWindTunnel -> Check Boundaries -> TowerNacelle -> Mesh Conditions -> Custom Surface size Regions -> VirtualWindTunnel -> Relative Minimum Size -> 2.2 === Relative Boundaries -> TowerNacelle -> Mesh Target Size-> 2.2 Values -> Surface size Regions -> VirtualWindTunnel -> Check Boundaries -> Rotating Part -> Mesh Conditions -> Custom Surface size 94 Regions -> VirtualWindTunnel -> Relative Minimum Size -> 11 === Relative Boundaries -> Rotating Part -> Mesh Values Target Size-> 11 -> Surface size Generate Volume Mesh from Here Click Derived Parts (right click) -> New part -> x = 0, 1; y = 0, 0; z = 0,0 === No displayer Section -> Plane Scenes -> MeshVirtualWindTunnel -> Outline -> Uncheck === Surface -> Check Displayers (right click) -> New displayer -> Geometry Scenes -> MeshVirtualWindTunnel -> Regions -> VirtualRotatingSubdomain -> Hub, Displayers -> Geomatry 1 -> Parts (right Left Blade, Right Blade, Left Blade_Base, click) Right_Blase_Base === Regions -> VirtualWindTunnel -> Boundaries -> TowerNacelle Scenes -> MeshVirtualWindTunnel -> Color Mode -> Constant Displayers -> Geomatry 1 Definition of Physics Continua -> Physics 1 -> Models (right click) Time -> Implicity Unsteady === Material -> Gas -> Select Models === Flow -> Segregated flow === Equation of state -> Constant density === Viscous Regime -> Turbulent === Turbulence -> RST === Reynolds Stress Turbulence Models -> Elliptic Blending Continua -> Physics 1 -> Initial conditions -> y = 25 Velocity Tools -> Motions -> Rotation -> Axis y = -1 direction Tools -> Motions -> Rotation -> Rotation 72 rpm Rate Tools -> Coordinate Systems -> Laboratory - Default > Local Coordinate System (right click) -> Cartesian Local Coordinate System -> Cartesian1 Default Rotation -> Manage coordinate system Local Coordinate System -> Cartesian1 Assigning this rotation to the VirtualRotatingSubDomain Definition of boundary conditions Bottom, Top, Side, Side2 Symmetry Plane Inlet Velocity Inlet 95 Physics values -> velocity magnitude 25 m/s Outlet Pressure outlet Physics values -> pressure 0 (absolut pressure) Definition of the Solver Settings Implicity Unsteady -> Time step 0.001157 s Implicity Unsteady -> Temporal 2nd order Discretization Stopping criteria Maximum inner iteration 10 Maximum physical time 3.33s(equivalent to 4 complete rotor rotations) Maximum steps Uncheck Report->new report->Force Rename Force to Thrust Thrust->Axis [0.0,1.0,0.0] Thrust->Smooth values Check Thrust->Parts Left Blade, Right Blade Report->new report->Moment Rename Moment to Torque Moment->Axis [0.0, -1.0,0.0] Moment->Smooth values Check Moment->Parts Left Blade, Right Blade Select Moment and Thrust, right click and Multiple plot create Monitor plot Creating a scene Derived parts->New Part->Isosurface Create Isosurface->Parts Set all in Regions Isosurface->Scalar field Q-criterion Isosurface->Value 5 96 Scene->new scene->scalar Create Scalar-> displayers->outline1->outline Uncheck Scalar-> displayers->outline1->surface Check Scalar-> displayers->outline1->parts select all the wind turbine geometry Scalar->displayers->scalar1->countour style smooth filled Scalar->displayers->scalar1->parts isosurface Scalar->displayers->scalar1->scalar field- velocity ->magnitude >function