원자로 수치해석 과제 #1

: Matrix Eigenvalue Problem of 1D, 1G Neutron Diffusion Equation

2008-11833. 반영석

1. 1-D, 1-G neutron diffusion equation의 해석적인 해를 구하고 eigenvalue [eqn: lambda]의 참값을 구하여라.

[eqn: - d over dx D(x) d over dx d phi(x) + SIGMA_a (x) phi(x) = lambda nu SIGMA_f (x)phi(x)]   , [eqn: x in [0,H]]

위와 같이 나타내어지는 neutron diffusion equation을 풀어본다. 여기서 [eqn: SIGMA_a], [eqn: nu SIGMA_f], [eqn: D(x)], [eqn: H]값은 모두 주어진 상수이다. 위의 방정식은 다음과 같이 모양을 바꿀 수 있다.

[eqn: {d ^{2} phi } over {dx ^{2}} +[-{SIGMA_a - lambdanuSIGMA_f } over {D}]phi= {d^2 phi} over{ d x^2}+B^2 phi =0]

단, [eqn: B^2 ==-{SIGMA_a - lambdanuSIGMA_f } over {D}]이다. 위와 같은 꼴의 방정식을 Helmholtz Equation이라고 하는데 [eqn: B^2]값만 구하면 우리가 원하는 eigenvalue [eqn: lambda]를 구할 수 있다. 위의 방정식의 해는 삼각함수가 만족시키는데 [eqn: phi(x)=C_1 cos Bx + C_2sin Bx]라고 두고 경계조건 [eqn: phi_0]와 [eqn: phi_N]가 0임을 이용하면 곧 [eqn: phi(x)= C* sin Bx], [eqn: B= pi over H]임을 알 수 있다. 이로부터 역으로 [eqn: lambda= 1 over {nu SIGMA_f} ( SIGMA_a + D B^2 )]에서 원하는 값인 eigenvalue [eqn: lambda]를 얻는다.

[eqn: lambda= 1 over {nu SIGMA_f} ( SIGMA_a + D ({pi over H })^2 )]

실제 이 값은 9.71180e-001 이다.

point scheme을 적용할 때, discretization에 의해 다음과 같이도 쓸 수 있다.

[eqn: phi _{k-1} -2 phi _{k} + phi _{k+1} +h ^{2} B ^{2} phi _{k} =0]

또는 [eqn: phi _{k-1} -2 phi _{k} + phi _{k+1} +tilde B ^{2} phi _{k} =0] (여기서 [eqn: h^2 B^2 = tilde B ^2])

여기서 [eqn: phi_k = e^mk]라고 둔다면 [eqn: cosh m = 1 -{ tilde{ B } }^2 over 2]이 나오는데 여기서 좌변 [eqn: cosh m]은 항상 1 이상이고 우변은 항상 1보다 작으므로 해가 될 수 없다. 그러나 [eqn: phi_k = e^imk]라고 둔다면 [eqn: cosm= 1 - {tildeB}^2 over 2]가 되어 [eqn: phi_k]는 [eqn: C_1 cosmk + C_2 sin mk]꼴로 나타나게 된다. Boundary condition인 양 쪽 끝에서 flux가 zero가 되는, 즉 [eqn: phi_0]와 [eqn: phi_N]가 0임을 이용하면 [eqn: C_1 =0][eqn: m= {l pi} over N]라는 식이 나오는데 또 여기서 flux는 양의 값만 물리적인 의미가 있으므로 [eqn: phi]는 오직 fundamental mode ([eqn: l=1]) 에서만 유효하다. 위에서와 같은 방법으로 따라가면 error는 [eqn: h^2]에 비례함도 알 수 있다.

2. point scheme을 이용하여 N*N의 행렬들을 구성하여라.

Point scheme을 이용하면 위의 1D, 1G neutron diffusion equation을 수치해법적인 방법으로 풀 수 있다. [eqn: M PHI = lambda F PHI] 꼴의 행렬식을 풀게 되는데 이를 위해서는 (N+1) by (N+1)의tri-diagonal한 matrix [eqn: M]과 같은 크기인 diagonal matrix [eqn: F]가 필요하다. 둘의 구조는 다음과 같다.

표  M의 구조

[eqn: SIGMA_1^D + alpha_L over h + 1 over 2 SIGMA_a,1]	[eqn: -SIGMA_1^D]
[eqn: -SIGMA_1^D]	[eqn: SIGMA_1^D + SIGMA_2^D + barSIGMA_a,1]	[eqn: -SIGMA_2^D]
	[eqn: -SIGMA_2^D]	...
			...	[eqn: -SIGMA_N^D]
			[eqn: -SIGMA_N^D]	[eqn: SIGMA_N^D + alpha_R over h + 1 over 2 SIGMA_a,N]

(단, [eqn: SIGMA_n^D = D_n over h^2]이다.)

표  F의 구조

[eqn: 1 over 2 nu SIGMA_f,1]
	[eqn: bar {nu SIGMA_f,1]
		...
			[eqn: bar {nu SIGMA_f,N-1]
				[eqn: 1 over 2 nu SIGMA_f,N]

여기서 빈칸은 0을 뜻한다. 또 매질의 균질성에 의해 [eqn: SIGMA_1 ^D = SIGMA_2^D = ``cdots]이고 [eqn: bar{nu SIGMA_f,n]도 모두 동일하다. 이에 근거하여 행렬을 구성해 보았다.

function [M]= con_M(N)

% constant

sig_a=0.1;

D=2;

H=100;

% N : number of meshes

M=zeros(N+1,N+1);

h=H/N; % mesh size

sig_D=D/(h^2);

for i=2:N+1

M(i,i)=2*sig_D+sig_a;

M(i,i-1)=-sig_D;

M(i-1,i)=-sig_D;

end

M(1,1)=10^20;

M(N+1,N+1)=10^20;

function [F]=con_F(N)

% constant

sig_f=0.105; %nu sig_f

F=sig_f*eye(N+1,N+1);

F(1,1)=F(1,1)/2;

F(N+1,N+1)=F(N+1,N+1)/2;

function [lam]=solve(N)

M=con_M(N);

F=con_F(N);

A=inv(F)*M;

%F^-1 M 의 minimum eigenvalue

lam=min(eig(A));

function [LH,LE,A]=draw()

for i=1:8

A(i)=solve(i*50);

E(i)=abs(A(i)-real_eigen);

end

N=50:50:400;

H=100./N;

LE=log(E);

LH=log(H);

plot(LH,LE)

행렬 [eqn: M]의 (1,1)과 (N+1, N+1)성분에는 [eqn: alpha]값이 들어가는데 zero-flux조건에 의해 그 항이 무한대로 발산하므로 다른 항을 무시하고 [eqn: 10^20]을 입력하였다.

3. N이 50, 100, 150, 200 ... 일 때의 eigenvalue [eqn: lambda]의 최소값을 구하여라. mesh크기에 따라서 오차율이 어떻게 변하는지 보고 수업시간에 해석적으로 거동을 살핀 것과 비교하여라.

# of meshes	minimum eigenvalue	error
50	9.71174E-01	6.18389E-06
100	9.71179E-01	1.54613E-06
150	9.71180E-01	6.87179E-07
200	9.71180E-01	3.86541E-07
250	9.71180E-01	2.47387E-07
300	9.71180E-01	1.71796E-07
350	9.71180E-01	1.26218E-07
400	9.71180E-01	9.66350E-08

위의 표는 2.에서의 draw라고 정의한 함수로 구한 mesh의 개수에 따른 eigenvalue와 위에서 해석적으로 구한 eigenvalue의 참값 [eqn: lambda*]와의 차이 즉, error를 기록한 표이다.

위의 그래프는 [eqn: ln(h) - ln( DELTA lambda)] 그래프이다. 위에서 살펴보았듯이 Error([eqn: DELTA lambda])는 mesh size의 제곱에 비례함을 알 수 있었다.