hi,
I am a beginner to matlab. I have written a function code for the implementation of Crank Nicolson Method within which i hv called a script named Compressibility as given below. But weneva i run this code matlab shows its busy…for a long time. When i press ctrl+C the cmmand window shows the execution was stopped by user at line 44..evrytime i press ctrl+c. the codes r given below. How can I run the program n get my desired output circumventing this problem??? Plzz plzz help!!
function [ P_out ] = Cranknicolson(H,L,D,T,dt,dx,I,N,Q_init,P_in,rho)%
%%parameters needed:
g = 9.81; %acc. due to gravity(m/s2)
%height of receiving end(m)=H
%Length of gas duct(m)=L
f = 0.003; %friction factor
%diameter of gas duct(m)=D
R = 392; %m2/(s2.K)
%Temperature(K)=T
%density of gas(kg/m3)= rho
%Initial flow=Q_init
%Input pressure=P_in
Z=1; %initial guess of Z value
compressibility;c = sqrt(Z*R*T);% C is the iso-thermal speed of sound
r = sqrt(L.^2+H.^2);Y = H/r; %value of sine(theta)
P_out = zeros(N,1);A = zeros((I-1),(I-1));C = zeros((I-1),(I-1));G = zeros((I+1),N);G1 = zeros((I-1),1);G(1,:) = P_in.^2;alpha = (16*f*Q_init)/(D.^3*c.^2*pi);beta = (2*g*Y)/c.^2;a = dt/(alpha*dx.^2);b = (alpha*beta*dx)/2;A1 = a-b;A2 = -2*a-2;A3 = a+b;A4 = b-a;A5 = 2*a-2;A6 = -a-b;for i=1:(I-1) A(i,i)=A2; C(i,i)=A5; while i>1 A(i-1,i)=A3; A(i,i-1)=A1; C(i-1,i)=A6; C(i,i-1)=A4; endendfor i=2:(I+1) x=(i-1)*dx; G(i,1) = (G(1,1)+(xi/sigma))*exp(-sigma*x)-(xi/sigma);endG1(1:(I-1),1)=G(2:I,1);for n=1:(N-1) G(I+1,n+1)=(G(I,n)+(xi/sigma))*exp(-sigma*dx)-(xi/sigma); B = [ (A4*G(1,n) - A1*G(1,n+1));zeros(I-3,1);(A6*G(I+1,n) - A3*G(I+1,n+1)) ]; G2 = A\((C*G1)+B); G1(1:I-1,1) = G2(1:I-1,1); G(2:I,n+1)=G2(1:I-1,1);endP_out(1:N,1)=sqrt(G(I+1,1:N));end**Compressibility**:%To find out the value of Z(compressibility factor)
%%parameters needed:
Qn0 = 50; %initial flow(m3/s)
P_0 = 50e+5; %input pressure(pascal)
%Other values
r = sqrt(L.^2+H.^2);Y = H/r;%Z0 = 1; %initial value
c = sqrt(Z0*R*T);sigma = (2*g*Y)/c.^2;xi = (f/D)*(((8*rho*c*Qn0)/(pi*D.^2)).^2); P_L = sqrt((P_0.^2 + (xi/sigma))*(exp(-sigma*L))-(xi/sigma)); P_avg = 0.5*(P_0+P_L); Z = (1 - (P_avg/39e6));while abs(Z-Z0)>1e-6 Z0=Z; c = sqrt(Z0*R*T);sigma = (2*g*Y)/c.^2;xi = (f/D)*(((8*rho*c*Qn0)/(pi*D.^2)).^2); P_L = sqrt((P_0.^2 + (xi/sigma))*exp(-sigma*L)-(xi/sigma)); P_avg = 0.5*(P_0+P_L); Z = (1 - (P_avg/39e+6));endfprintf('The value of Z is %f\n',Z)
Best Answer