I am trying to calculate and plot subsonic coefficients of pressure symbolically, but it is returning a piecewise function that I can't easily plot, what am I doing wrong?

8 views (last 30 days)
I am trying to calculate and plot subsonic coefficients of pressure symbolically, but it is returning a piecewise function that I can't easily plot, am I doing something?
clear all; clc; close all
syms x ksi
%N = 20;
gamma = 1.4;
limits = [0 1];
M = 1.2;
%alpha = deg2rad([-5:1:20]);
%X = linspace(0,1,N);
%Geometry
c = 1;
tau = 1/10; %thickness ratio
ZksiU = 2*tau*c*(ksi/c-(ksi/c)^2);
ZksiL = -0.5*tau*c*(ksi/c-(ksi/c)^2);
fun = diff(ZksiU,ksi,1)/(x-ksi);
Cpi = -(2/pi)*int(fun,[limits(1) limits(2)]) %incompressible pressure coefficient
CpSub = -(2/M^2*(gamma+1))*((1-M^2)-((1-M^2)^3/2+(3/4)*M^2*(gamma+1)*Cpi)^(2/3)) %Subsonic pressure coefficient
returns:
fun = ((2*ksi)/5 - 1/5)/(ksi - x)
Cpi = piecewise(ksi in Dom::Interval([0], [1]), (5734161139222659*((2*ksi)/5 - 1/5)*int(1/(x - ksi), x, 0, 1))/9007199254740992, ~in(ksi, 'real') | ksi < 0 | 1 < ksi, -(5734161139222659*((2*ksi)/5 - 1/5)*(log(-ksi) - log(1 - ksi)))/9007199254740992)
CpSub = piecewise(ksi in Dom::Interval([0], [1]), (10*((464467052277035379*int(-1/(ksi - x), x, 0, 1)*((2*ksi)/5 - 1/5))/281474976710656000 - 1331/31250)^(2/3))/3 + 22/15, ~in(ksi, 'real') | ksi < 0 | 1 < ksi, (10*(- (464467052277035379*((2*ksi)/5 - 1/5)*(log(-ksi) - log(1 - ksi)))/281474976710656000 - 1331/31250)^(2/3))/3 + 22/15)

Accepted Answer

Walter Roberson
Walter Roberson on 18 Oct 2021
syms x ksi
%N = 20;
gamma = 1.4;
limits = [0 1];
M = 1.2;
%alpha = deg2rad([-5:1:20]);
%X = linspace(0,1,N);
%Geometry
c = 1;
tau = 1/10; %thickness ratio
ZksiU = 2*tau*c*(ksi/c-(ksi/c)^2)
ZksiU = 
ZksiL = -0.5*tau*c*(ksi/c-(ksi/c)^2);
fun = diff(ZksiU,ksi,1)/(x-ksi)
fun = 
Cpi = -(2/pi)*int(fun,[limits(1) limits(2)]) %incompressible pressure coefficient
Cpi = 
CpSub = -(2/M^2*(gamma+1))*((1-M^2)-((1-M^2)^3/2+(3/4)*M^2*(gamma+1)*Cpi)^(2/3)) %Subsonic pressure coefficient
CpSub = 
KSI = linspace(-1, 2, 1000);
CP = double(subs(CpSub, ksi, KSI));
plot(KSI, real(CP), 'k-.', KSI, imag(CP), 'r-+')
vpa(subs(CpSub, 0.5))
ans = 
NaN
Why is it getting NaN for a fair bit of the region? Well, you are asking to integrate 1/(x-ksi) with ksi in the range 0 to 1, and x in the range 0 to 1. Except for the case of ksi = 0 exactly or ksi = 1 exactly, you are trying to integrate through a fundamental discontinuity, equivalent to integrating through 1/X with X being 0 at some point in the range.
Hence, the integral is only defined for ksi outside the 0 to 1 range.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!