% Euler_template V1.0 % Final Project % % THP 04/2K+1 % Linearized Euler Code for A Thin Circular Arc Airfoil % % Flow is Assumed to be rho = 1, u = M_inf , v = 0 at Inflow % Linearization is about rho = 1, u = M_inf, v = 0 % Thin Airfoil BC % % Euler equations are written as % ( Pressure = rho is used i.e. Constant Temperature) % % DQ DQ DQ % -- + A -- + B -- = 0 % Dt Dx Dy % % where: | 0 1 0 | | rho | % A = | -M_inf^2 2 M_inf 0 |, Q = | rho u | % | 0 0 M_inf | | rho v | % % | 0 0 1 | % B = | 0 0 M_inf | % | 1 0 0 | % M_inf : Free Stream Mach number % nx : Number of points on airfoil (x=0,1) % % Inputs: % cfl : dt = min(dx,dy)*cfl/(M_inf+1 + 1) % where Spectral radius of A is (M_inf+1) and B is (1) % nmax : Number of time steps to go % clear; M_inf = 0.8; % remove comment % $$$ M_inf = input('Enter M_inf = '); % Grid Generation xmin = -1.0; xmax = 3.0; ymin = 0.0; ymax = 2.0; nx = 10; % $$$ nx = input('Enter number of points on body = '); dx = 1.0/(nx-1); dy = dx; x = [xmin:dx:xmax]'; y = [ymin:dy:ymax]'; jmax = size(x,1); kmax = size(y,1); [X,Y] = meshgrid(x,y); jle = nx; jte = 2*nx-1; fprintf(' X of Leading edge point = %g X of trailing edge point = %g \n',x(jle),x(jte)); % Surface Definition tau = 0.1; % tau = input('Enter tau = '); ys = zeros(jmax,1); s = [jle:jte]; ys(s) = tau*0.5*x(s).*(1.0-x(s)); % Initialize Q's Q(:,:,1) = ones(jmax,kmax,1); % rho Q(:,:,2) = M_inf*ones(jmax,kmax,1); % rho u Q(:,:,3) = zeros(jmax,kmax,1); % rho v % Jacobian Matrices A = zeros(3,3); B = zeros(3,3); A(1,2) = 1.0; A(2,1) = -M_inf^2 + 1; A(2,2) = 2.0*M_inf; A(3,3) = M_inf; B(1,3) = 1; B(2,3) = M_inf; B(3,1) = 1; % Thin airfoil BC v = M_inf D(ys)/Dx Q(s,1,3) = M_inf*(ys(s+1) - ys(s-1))*0.5/dx; % Indices J = [1:jmax]; K = [1:kmax]; JM = [2:jmax-1]; KM = [2:kmax-1]; % Inputs cfl = input('Enter cfl = '); nmax = input('Enter nmax = '); dt = min(dx,dy)*cfl/(2.0+M_inf); fprintf(' AT CFL = %g dt = %g \n',cfl,dt); for nstep = 1:nmax % rhs Central differencing Qx_c(JM,KM,:) = 0.5*(Q(JM+1,KM,:) - Q(JM-1,KM,:))/dx; Qy_c(JM,KM,:) = 0.5*(Q(JM,KM+1,:) - Q(JM,KM-1,:))/dy; for m = 1:3 R(JM,KM,m) = 0.0; for n = 1:3 R(JM,KM,m) = R(JM,KM,m) + ... A(m,n)*Qx_c(JM,KM,n) + ... B(m,n)*Qy_c(JM,KM,n) ; end end % $$$ You can add these next 4 lines and run the code with cfl = 0.1 % $$$ for 1000 steps to see what a reasonable solution would be. % $$$ I added this in to give yuou at a sample working code. % $$$ Remove % line# % line1 epse = 0.1; % line2 Qx_d(JM,KM,:) = (Q(JM+1,KM,:) -2*Q(JM,KM,:) + Q(JM-1,KM,:))/dx; % line3 Qy_d(JM,KM,:) = (Q(JM,KM+1,:) -2*Q(JM,KM,:) + Q(JM,KM-1,:))/dy; % line4 R(JM,KM,:) = R(JM,KM,:)-epse*(Qx_d(JM,KM,:) + Qy_d(JM,KM,:)); resid(nstep) = ... (norm(R(:,:,1)) + norm(R(:,:,2)) + norm(R(:,:,3)))/(3*jmax*kmax); if nstep == 1 resid_norm = resid(1); end resid(nstep) = resid(nstep)/resid_norm; fprintf(' nstep = %g Resid = %g \n',nstep,resid(nstep)); Q(JM,KM,:) = Q(JM,KM,:) - dt*R(JM,KM,:); % Euler Explicit Step % BC % At j = 1 if M-inf < 1 Extrap rho if M_inf < 1.0 Q(1, KM,1) = Q(2, KM,1); end % At jmax DQ/Dx = 0 using backward diff. Q(jmax,KM,1) = Q(jmax-1, KM,1); Q(jmax,KM,2) = Q(jmax-1, KM,2); Q(jmax,KM,3) = Q(jmax-1, KM,3); % At k=1 D rho/Dy = 0, Du/Dy = 0, fix v RHO = Q(JM,1,1); Q(JM,1,1) = Q(JM,2,1); Q(JM,1,2) = Q(JM,2,2)./Q(JM,2,1).*Q(JM,1,1); % fixing v, not rho v so this little rescaling is necessary Q(JM,1,3) = Q (JM,1,3)./RHO.*Q(JM,1,1); % At kmax fix all do nothing end % Fix up corners for plot Q(1,1,:) = Q(2,1,:); Q(jmax,1,:) = Q(jmax-1,1,:); Q(1,kmax,:) = Q(2,kmax,:); Q(jmax,kmax,:) = Q(jmax-1,kmax,:); figure; subplot(2,2,1); surf(X,Y,Q(:,:,1)'); xlabel('x'); ylabel('y'); zlabel('rho'); view(-25,50); subplot(2,2,2); U = Q(:,:,2)./Q(:,:,1); V = Q(:,:,3)./Q(:,:,1); surf(X,Y,U'); xlabel('x'); ylabel('y'); zlabel('u'); view(-25,50); subplot(2,2,3); P = Q(:,2,1); cp = (P - 1.0)/(0.5*M_inf^2); plot(x,-cp); xlabel('x'); ylabel('-cp'); subplot(2,2,4); semilogy([1:nmax],resid); xlabel('n'); ylabel('resid'); title(['Euler Explicit',' CFL = ',num2str(cfl),... ' M_{\infty} = ',num2str(M_inf)]);