شبیه سازی معادلات توزیع دما در راکتور هسته ای در MATLAB

[پیرجو]

شبیه سازی معادلات توزیع دما در راکتور هسته ای در MATLAB

global n delta rspan options r T  
n = 100;
delta = 1.e-3;
rspan=linspace(delta,1,n+1);
options = odeset('RelTol',1e-12,'AbsTol',1e-12,'NormControl','on','MaxOrder',5);
T0=newton_mfile(@cp6_BV,2,1.e-10,0); 
r=[0;r];  
T=[[T0,0];T]; 
[Tdirect,rdirect] = cp6_direct(n);
figure(1)
hold off
Texact= 9/32-r.^2/4-r.^4/32;
plot(r,Texact,'r-')
hold on
grid on
plot(r,T(:,1),'b:');
plot(rdirect,Tdirect,'g--');
legend('Exact','Shooting','Direct')

%----------------------------------------------------

%newton_mfile

function x_root=newton_mfile(funhandle,x1,eps,reg)
del = max([1.e-10,eps^2]);
x_guess(1) = [x1];
f_guess(1) = feval(funhandle,x_guess(1));
if abs(f_guess) < eps
   x_root=x_guess;
   disp('excellent guess');
   return
end 
iter = 1;
check = logical(1);
tau=0.1;
while check
    iter = iter +1;
    x_guess(iter)=0;
    f_guess(iter)=0;
    f_h = feval(funhandle,x_guess(iter-1)+del);
    df=(f_h-f_guess(iter-1))/del;    
    dx(iter) = f_guess(iter-1)/df; 
    if iter >=  3
        tau=tau*dx(iter-1)/(dx(iter-1)-dx(iter));
    end
    if reg == 0
        x_guess(iter)=x_guess(iter-1)-dx(iter);
    else
        x_guess(iter)=x_guess(iter-1)-tau*dx(iter);
    end
    f_guess(iter) = feval(funhandle,x_guess(iter));
    if reg == 0 & abs(tau) < 1.0e-5 | abs(f_guess(iter)) < eps 
        x_root=x_guess(iter);
        check = logical(0);
    end
end
    
%-------------------------------------------------------

%cp6_BV

function f = cp6_BV(z)
global n delta rspan options r T
Tdelta  = z-1/4*delta^2-delta^4/32; 
dTdelta = -1/2*delta-delta^3/8;
[r,T] = ode15s(@rhs_cp6,rspan,[Tdelta dTdelta],options); 
f = T(end,1);  
%------------------------------------------------------------

%cp6_direct

function [T,r] = cp6_direct(N)

A = zeros(N+1);
S = zeros(N+1,1);
r = linspace(0,1,N+1)';
S(2:N) = -1 - r(2:N).^2/2;  
h = 1.0/N;
for i = 2:N 
    A(i,i-1) = 1/h^2-1/(2*h*r(i));
    A(i,i)   = -2/h^2;
    A(i,i+1) = 1/h^2+1/(2*h*r(i));
end
A(1,1:2) = [-1,1]; 
A(N+1,N+1) = 1;
T = inv(A)*S;