1. Implementação do Método dos Painéis

O programa, cujo o código implementado no MATLAB segue descrito abaixo,

efetua o cálculo do coeficiente de arrasto, o momento de arfagem sobre o centro

aerodinâmico, além de calcular a distribuição do coeficiente de pressão, a distribuição

de velocidade e a distribuição de intensidade dos vórtices ao redor do aerofólio. O

contorno do aerofólio é aproximado por painéis de vórtices cujas intensidades variam

linearmente.

As variáveis de entrada e saída para o programa são:

Tabela 1 – Relação de Entradas e Saídas do código desenvolvido.

Entrada

XB Vetor que contém as coordenadas x das seções do aerofólio

YB Vetor que contém as coordenadas y das seções do aerofólio

ALPHA Ângulo de Ataque em graus

Saída

cl Coeficiente de arrasto da seção

cmac Coeficiente do momento de arfagem sobre o centro aerodinâmico

xp Localização em x do ponto central do painel

yp Localização em y do ponto central do painel

cp Vetor contendo os coeficientes de pressão de cada painel no ponto (xp,yp)

V

Vetor contendo a velocidade dividida pela velocidade do fluxo livre para

cada painel no ponto (xp,yp)

gamma

Vetor contendo a circulação normalizada de cada painel. Os valores de

gamma são a força do vórtice dadas por (XB,YB)

Para implementação do programa, deve-se lembrar que os vetores XB e YB

precisam possuir as mesmas dimensões e de que o numero de painéis é dada pela

dimensão desses vetores subtraída de uma unidade. Deve-se contar também que o

primeiro e o ultimo elemento devem estar localizados no bordo de fuga do aerofólio.

Para efeitos de teste do código que segue abaixo utilizaremos o seguinte

aerofólio para efetuar os cálculos:

X Y1 0

0.933 -0.0050.75 -0.0170.5 -0.033

0.25 -0.0420.067 -0.033

0 00.067 0.0450.25 0.0760.5 0.072

0.75 0.0440.933 0.013

1 0

Figura 1 – Dados do aerofólio utilizado para teste do código.

Na tela de comando do MATLAB, devem então ser digitados os seguintes

parâmetros para chamar a função.

>XB=[1 0.933 0.75 0.5 0.25 0.067 0 0.067 0.25 0.5 0.75 0.933 1];

>YB=[0 -0.005 -0.017 -0.033 -0.042 -0.033 0 0.045 0.076 0.072 0.044 0.013 0];

>ALPHA=8;

>[cl,cmax,xp,yp,cp,V,gamma]=airfoil(XB,YB,ALPHA)

Inicia-se abaixo o código escrito para o MATLAB para o método dos painéis.

Vale salientar que será utilizado o comando do MATLAB para comentários, o %, para

realizar as explicações no decorrer do código. O código deverá ser salvo como

airfoil.m

function [cl,cmac,xp,yp,cp,V,gamma]=airfoil(XB,YB,ALPHA)%Caracterizando o Aerofólio:%Verifica os tamanhos dos vetores if size(XB,1)~=1 | size(YB,1) ~=1

error ('Erro na especificação do aerofolio : XB e YB precisam ter as mesmas dimensões')end %Identifica o valor minimo para x e sua posição [XMIN,ILE]=min(XB); %Achando o valor mínimo de x e sua posiçãoMP1=size(XB,2); %Lendo o número de pontos dadosM=MP1-1; %Número de PainéisALPHA1=ALPHA*pi/180; %Calculando angulo em radianos %Verificando os valores de y para que o o aerofólio seja fechado e que%a condição de Kutta possa ser testada (y final e y inicial no mesmo ponto) if (YB(1)~=0|YB(MP1)~=0)error ('Erro na Especificação do Aerofólio : Y(0) e o Y(end) precisam ser 0')end %Verificando se quando X apresenta seu valor mínimo, y está no eixo x if (XB(1)~=XB(MP1)|YB(ILE)~=0)error ('Erro na Especificação do Aerofólio')end %Plotando o Aerofólioxmin=min(XB);xmax=max(XB);ymin=10*min(YB);ymax=10*max(YB);plot(XB,YB);axis([xmin,xmax,ymin,ymax])title('Aerofólio', 'Color', 'b')input('Pressione qualquer tecla para continuar'); %Passando os valores de linhas para colunas XB=XB'; YB=YB'; %Calculando a corda c=XB(1)-XMIN; XAC=XMIN+c/4; %Criando os vetores com dimensões adequadas e reservando espaço na memória%para que os cálculos possam ser efetuados RHS=zeros(MP1,1);gamma=zeros(MP1,1); AN=zeros(MP1,MP1); AT=zeros(M,MP1);THETA=zeros(M,1); X=zeros(M,1); Y=zeros(M,1); S=zeros(M,1);SX=zeros(M,1); SY=zeros(M,1); SINE=zeros(M,1);

COSINE=zeros(M,1); %Calculando as coordenadas centrais do painel X=(XB([1:M])+XB([2:MP1]))/2; Y=(YB([1:M])+YB([2:MP1]))/2; %Calculando as distâncias relativas a x e y para 2 pontos consecutivos SX=XB([2:MP1])-XB([1:M]) SY=YB([2:MP1])-YB([1:M]); %Calculando a distância entre dois pontos consecutivos S=sqrt(SX.^2+SY.^2); %Calculando o angulo para 4 quadrantes e os respectivos senos e cosenos THETA=atan2(SY,SX); SINE=sin(THETA);COSINE=cos(THETA); for I=1:Mfor J=1:Mif I==JCN1(I,J)=-1.0;CN2(I,J)=1.0;CT1(I,J)=pi/2;CT2(I,J)=pi/2;elseXTEMP=X(I)-XB(J);YTEMP=Y(I)-YB(J);A=-XTEMP*COSINE(J)-YTEMP*SINE(J);B=XTEMP^2+YTEMP^2;TTEMP=THETA(I)-THETA(J);C=sin(TTEMP);D=cos(TTEMP);E=XTEMP*SINE(J)-YTEMP*COSINE(J);F=log(1+S(J)*(S(J)+2*A)/B);G=atan2(E*S(J),B+A*S(J));TTEMP=TTEMP-THETA(J);P=XTEMP*sin(TTEMP)+YTEMP*cos(TTEMP);Q=XTEMP*cos(TTEMP)-YTEMP*sin(TTEMP);CN2(I,J)=D+Q*F/(2*S(J))-(A*C+D*E)*G/S(J);CN1(I,J)=D*F/2+C*G-CN2(I,J);CT2(I,J)=C+P*F/(2*S(J))+(A*D-C*E)*G/S(J);CT1(I,J)=C*F/2-D*G-CT2(I,J);endendend % Calculando os Coeficientes de InfluênciaAN=[CN1(:,1),CN1(:,[2:M])+CN2(:,[1:M-1]),CN2(:,M);...1,zeros(1,M-1),1];AT=[CT1(:,1),CT1(:,[2:M])+CT2(:,[1:M-1]),CT2(:,M)];RHS=[sin(THETA-ALPHA1);0]; % Calculando a circulação

gamma=AN\RHS; % Calculando a velocidade V=cos(THETA-ALPHA1)+AT*gamma; % Calculando o coeficiente de Pressão cp=1-V.^2; % Calculando a circulação total G=S'*(gamma([1:M])+gamma([2:MP1]))/2; % Calculando o coeficiente de arrasto cl=4*pi*G/c; % Calculando o momento de arfagem sobre o centro aerodinâmico cmac=S'*(cp.*((X-XAC).*COSINE+Y.*SINE)/c^2); xp=X; yp=Y; % Organizando dados para plotarCpplot=zeros(M,2);Cfp=[X,cp]; for i=1:M; P=(M/2)+i; if i>(M/2); P=i-(M/2); end Cpplot(i,:)=Cfp(P,:);end %Plotando plot(Cpplot(:,1),-Cpplot(:,2))xlabel('Corda', 'Color', 'b')ylabel('-cp', 'Color', 'b')title('Distribuição dos coeficientes de pressão', 'Color', 'b')