수평원기둥 파력발전장치의 시간영역 운동해석을 통한 연간 에너지생산량 추정

1. 서 론

입사하는 횡파에 의해 횡 운동하는 파력발전장치의 로터의 형상은 Salter’s duck과 같이 비대칭 단면 형상을 갖도록 설계하거나(Salter[1974]; Salter et al.[1975]), 원형 대칭 단면에 비대칭 형상과 동일한 효과를 주기 위하여 회전축을 중심에서 벗어나게 설계한다(Lucas et al.[2009]).

Cho et al.[2018a; 2018b]은 파랑중 횡 운동하는 수평원기둥 로터를 파력발전장치의 1차 변환장치로 채택하고 로터의 회전축을 원기둥 중심으로부터 이격시켰을 때 로터의 횡 운동 특성을 수치 계산과 모형실험을 통하여 살펴보았다. 로터의 운동성능에 큰 영향을 주는 설계변수는 잠긴 깊이와 회전축의 위치이며, 이 값들을 변화시켜가면서 최대 추출효율을 주는 최적의 값들을 도출하였다. Ko et al.[2018a]은 수평원기둥 로터에 대하여 점성과 비선형 효과를 함께 고려하기 위해 OpenFOAM을 이용하여 2차원 단면에 대한 수치해석을 수행하였다. 수치계산 결과를 검증하기 위하여 자유감쇠실험 및 규칙파 모형실험 결과와 비교한 결과, 입사파의 파고가 비교적 낮을 때는 수치결과와 모형실험결과가 서로 잘 일치하였다. 또한 검증된 수치코드를 이용하여 파고 변화 및 회전축 위치 변화에 따른 로터의 횡 운동의 비선형적 특성에 대하여 살펴보았다.

파력발전장치를 통한 연간 에너지 생산량(Annual Energy Production, AEP) 평가는 현재 개발된 수치프로그램을 가지고 가능하며 본 연구의 궁극적인 목표중의 하나이다. 연간 에너지 생산량은 LCOE(Levelized cost of energy) 계산에 이용되며 이는 개발된 파력발전장치의 상용화 가능성을 평가하는 중요한 지표가 된다. 파력발전장치의 연간에너지 생산량은 설치 해역의 파랑통계 특성, 장치의 발전성능, 추출파워 매트릭스, 1·2차 변환장치의 에너지 손실, 송전에 의한 손실 등으로 결정된다. 파랑통계 특성은 특정 유의파고(Significant wave height)와 피크주기(Peak period)의 발생 빈도로 나타나는 파랑산포도(Wave scatter diagram)로 표현된다. 불규칙한 해역에서의 파력발전장치의 발전성능은 특정한 파랑조건에서의 추출파워에 파랑산포도에 나타난 발생 빈도수를 곱한 기계파워 매트릭스로 표현된다. 여기에 2차 변환장치의 PTO(Power take-off) 변환효율을 곱하면 시간평균 전력파워를 구할 수 있다.

최종적으로 송전효율과 정격출력에서 크게 벗어난 최대 전력 출력을 제한하는 운용효율을 고려하면 연간 에너지 생산량을 추정할 수 있다.

본 연구는 수평원기둥 로터의 시간영역 운동해석을 실시하여 불규칙 파랑중 로터의 횡 운동 변위와 속도 그리고 추출파워를 구하였다. 이때 로터의 횡 운동에 의해 시시각각 변하는 횡 방향 정유체 모멘트(복원모멘트)를 부분 비선형, 완전 비선형으로 구분하여 계산하였다. 또한 회전축의 위치에 따라 달라지는 점성에 의한 감쇠력을 고려하였다. 점성에 의한 감쇠력을 얻기 위하여 회전축의 위치를 바꿔가면서 CFD 코드인 OpenFOAM을 사용하여 자유감쇠 수치실험을 실시하여 회전각도의 함수로 점성감쇠계수를 나타내었다. Cummins[1962]가 제안한 시간영역 운동방정식을 푸는 방법으로 OrcaFlex를 사용하였고, 이때 필요한 주파수영역에서의 부가 관성모멘트, 정유체 모멘트, 방사 감쇠계수, 파기진 모멘트는 WAMIT을 사용하여 계산하였다. 개발된 수치 프로그램을 이용하여 최대 파워를 보이는 최적의 회전축 위치를 정한 후 설치 대상 해역인 제주도 차귀도에서 얻은 파랑산포도에 나타난 파랑분포와 빈도수를 고려하여 기계파워 매트릭스를 구하였다. 여기에 2차 변환 효율을 곱하여 시간평균 전력파워를 구하고 이로부터 연간 에너지 생산량을 추정하였다.

2. 시간영역 해석

2.1 해석 모델

본 연구에서 다룰 수평원기둥 파력발전장치의 설치 대상 해역은 제주도 차귀도에 위치한 파력발전장치 실해역 시험장 5번 정박지로 가정하였고, 수심은 80 m이다.

수평원기둥 파력발전장치의 추출효율을 극대화하는 많은 설계변수들이 있다. Cho et al.[2018a]의 연구 결과를 인용하여 수평원기둥로터의 잠긴 깊이(d)를 1.6 m로 하였고, 원기둥의 중심으로부터 반지름 방향으로 1.5 m (R = 0.75r) 떨어진 곳에 회전축을 위치시켰다(Fig. 1). 그 밖의 수평원기둥 로터의 상세 제원은 Table 1과 같다(Cho et al.[2018a]).

Fig. 1.

Schematic diagram of a horizontal cylinder rotor in beam sea.

Table 1.

Specification of the horizontal cylinder rotor

로터의 회전축 위치를 나타내는 또 다른 설계변수인 설치각도(α)는 로터의 횡 운동 및 추출 파워에 큰 영향을 미친다. 따라서 최대 추출파워를 주는 최적의 회전축의 각도를 선정하는 것이 필요하다. 본 연구에서는 회전축 각도를 0도에서 360도까지 10도 간격으로 변화시켜가며 설치 대상해역의 대표 파랑스펙트럼(유의파고, 피크주기)을 갖는 불규칙파에 대해 수평원기둥 로터의 횡 운동 변위와 속도를 시간영역 해석법을 통하여 구하고, 이로부터 시간평균 추출파워를 계산하여 최적의 회전축 설치각도의 범위를 선정하였다. 이 때 계산시간을 줄이기 위하여 정유체 모멘트 계산시 완전 비선형(Fully nonlinear) 모델을 이용하지 않고 부분 비선형(Partially nonlinear) 모델을 사용하였다. 1차로 선정된 설치각도 범위내에서 완전 비선형 모델을 사용하여 로터의 최적의 회전축의 설치각도를 정하였다. 이렇게 선정된 최적의 위치에 회전축을 위치시킨 후 설치 대상 해역의 파랑산포도에 나타난 각 유의파고와 피크주기에 따라 수평원기둥 로터의 횡 운동 시간영역 해석을 반복 실시하였다. 각 파랑 조건에 대한 해석 결과를 이용하여 시간평균 파워를 계산하여 파워 매트릭스를 작성하고 파랑산포도의 파랑 발생 빈도수를 고려하여 시간평균 전력파워를 구하고, 운용 효율과 송전 효율을 곱하여 최종적으로 연간 에너지 생산량을 계산하였다. 수평원기둥의 횡 운동해석 시, 2차 변환장치인 유압장치에 무리가 가지 않도록 로터의 운동이 운동 구속범위인 -60도에서 +60도를 벗어나지 않게 수치적으로 제한하였다. 이를 위하여 로터의 운동 변위가 -60도와 +60도를 벗어날 경우에 다시 운동 구속 범위내로 돌아올 수 있도록 큰 복원력 값을 인위적으로 주었다. 이때 강제로 적용한 복원력이 로터의 운동에 큰 영향을 주지 않도록 값을 적절히 조정하였다.

2.2 시간영역 운동방정식

수평원기둥 로터는 회전축을 중심으로 횡 운동만 하며 다른 운동모드는 구속되어 있고 회전축에는 정적인 하중이 작용하지 않는다고 가정하였다. 시간영역에서 로터의 1 자유도 횡 운동방정식은 식 (1)과 같다(Cummins[1962]).

$\[ \begin{array}{l}\left(I_{55}+a_{55}\left(\infty \right)\right)\ddot{\xi }+\left(b_{vis}+{\tilde{b}}_{PTO}\left({\omega }_N\right)\right)\dot{\xi }\\ +\underset{0}{\overset{t}{\int }}{L}_{55}\left(\tau \right)\dot{\xi }\left(t-\tau \right)d\tau +K_{55}\xi =M_E\left(t\right)\end{array} \]$

(1)

여기서 I₅₅는 로터의 관성 모멘트(Moment of Inertia), α₅₅(∞)는 ω→∞일 때의 횡 방향 부가 관성 모멘트, b_vis는 점성과 마찰에 의한 감쇠계수, $$ {\tilde{b}}_{PTO}\left({\omega }_N\right)$$는 에너지 추출에 따른PTO 감쇠계수이다. L₅₅(τ)는 시간기억함수(time memory function), K₅₅는 정유체 모멘트 계수, M_E는 파기진 모멘트이다. $$ \xi , \dot{\xi }, \ddot{\xi }$$는 각각 로터의 횡 방향 변위, 속도, 가속도를 나타낸다.

2.2.1 파기진 모멘트

수평원기둥 로터가 불규칙파중에 놓여 있을 때 시간영역에서의 횡 방향 파기진 모멘트는 조화함수의 중첩원리를 사용하면 파랑스펙트럼과 주파수영역에서 구한 단위 입사파의 진폭당 횡 방향 파기진 모멘트(X₅(ω))로부터 얻어진다.

$\[ \begin{array}{l}{M}_E\left(t\right)=\sum _{n=1}^{N_w}\left|X_5\left({\omega }_n\right)\right|A_n\cos \left({\omega }_{n}t+x_n+2\pi {\epsilon }_n\right),\\ \begin{array}{cc}{A}_n+\sqrt{2S_{\zeta }\left({\omega }_n\right)\mathrm{\Delta }{\omega }_n},& \mathrm{\Delta }{\omega }_n=\left(1+0.2{\epsilon }_n\right)\mathrm{\Delta }\omega \end{array}\\ \begin{array}{cc}{\omega }_n={\omega }_{n-1}+\frac{1}{2}\left(\mathrm{\Delta }{\omega }_n+\mathrm{\Delta }{\omega }_{n-1}\right),& j=2,\dots ,N_w\end{array}\end{array} \]$

(2)

여기서 ω₁ = 0.1 rad/s, Δω = 3.0/N_w rad/s, N_w = 300, ε_n는 범위[0,1] 내의 난수를 뜻한다. 식 (2)에서 파랑스펙트럼을 주파수에 대하여 등 간격으로 나누지 않은 이유는 일정시간 후에 파기진 모멘트가 반복하는 것을 피하기 위함이다. 여기서 S_ζ(ω)는 파랑 스펙트럼으로 JONSWAP 스펙트럼을 사용하였다. JONSWAP 스펙트럼은 식 (3)과 같다.

$\[ \begin{array}{l}{S}_{\xi }\left(\omega \right)=\frac{\alpha g^2}{{\omega }^5}\exp \left[-1.25{\left(\frac{\omega }{{\omega }_P}\right)}^{-4}\right]{\gamma }^{\exp \left[-\frac{{\left(\omega -{\omega }_P\right)}^2}{2{\sigma }^2{\omega }_P^2}\right]},\\ \mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h} \alpha \approx 5.061{\left(\frac{{\omega }_P}{2\pi }\right)}^4{H_s}^2\left(1-0.287\log \gamma \right)\end{array} \]$

(3)

여기서 H_s는 유의파고, $$ {\omega }_P\left(=\frac{2\pi }{T_P}\right)$$는 피크주파수. γ는 Peak enhancement factor, g는 중력가속도이다. ω < ω_P일 때 σ = 0.07이며, ω ≥ ω_P일 때 σ = 0.09이다. 식 (2)에서 파기진 모멘트가 로터에 갑자기 작용하는 것을 막기 위하여 램프 함수(Ramp function) S(t) = 3(t/t_s) − 2(t/t_s)을 도입하였다. 여기서 t_s = 5T_P는 램프 함수가 적용된 시간이다.

2.2.2 PTO감쇠력

식 (1)의 PTO 감쇠력은 2차 변환장치인 PTO 장치 설치에 따른 감쇠력이다. 최적의 PTO 계수는 주파수의 함수로 식 (4)와 같이 쓸 수 있으나, 여러 개의 주파수들이 섞여 있는 불규칙파에서 PTO 감쇠계수를 주파수에 따라 매번 바꾸는 것이 불가능하므로 횡 운동이 가장 크게 일어나는 공진주파수에서의 값($$ {\tilde{b}}_{PTO}\left({\omega }_N\right)=b_{55}\left({\omega }_N\right)+b_{vis}$$)을 사용하였다 (Cho et al.[2018a]; Cho et al.[2018b]).

$\[ {\tilde{b}}_{PTO}\left(\omega \right)=\frac{\sqrt{{\left({\omega }_N^2-{\omega }^2\right)}^2{\left(I_{55}+a_{55}\left(\omega \right)\right)}^2+{\omega }^2{\left(b_{55}\left(\omega \right)+b_{vis}\right)}^2}}{\omega } \]$

(4)

여기서 공진주파수 $$ {\omega }_N=\sqrt{\frac{K_{55}}{I_{55}+a_{55}\left({\omega }_N\right)}}$$이다.

2.2.3 점성 감쇠력

점성 감쇠력은 포텐셜이론으로 해석이 불가능하므로 모형실험 혹은 점성 수치해석을 이용하여 구한다. Ko et al.[2018a]은 점성수치해석 프로그램인 OpenFOAM을 이용하여 수평원기둥 로터에 대한 자유감쇠 수치실험을 실시하였다. 수치해석 결과를 검증하기 위하여 자유감쇠 모형실험 결과와 비교하였다. 본 연구에서도 OpenFOAM을 이용하여 수평원기둥 로터의 회전축 각도 변화에 따라 자유감쇠 수치실험을 실시하여 점성 감쇠력을 구하여 곡선맞춤(Curve fitting)을 통하여 점성 감쇠력을 회전축 각도의 함수로 나타내 시간영역 운동해석에 활용하였다.

OpenFOAM을 이용하여 로터의 횡 방향 자유감쇠실험을 실시하면 Fig. 2와 같은 횡 운동 시계열 결과를 얻을 수 있다. 연속적인 4개의 파봉점 및 파곡점의 변위 값(z_ai, i = 1,2,3,4)으로부터 감쇠비(Damping ratio, κ)를 식 (5)와 같이 로그함수를 통해 구할 수 있다. 이렇게 얻은 감쇠비는 방사 감쇠력과 점성 감쇠력을 모두 포함한 전체 감쇠력이다. 여기서 방사 감쇠력을 빼면 순수한 점성 감쇠력을 구할 수 있다.

Fig. 2.

Determination of logarithmic decrement (Journée and Massie[2001]).

$\[ \kappa =\frac{1}{2\pi }\ln \left\{\frac{z_{a_1}+z_{a_2}}{z_{a_3}+z_{a_4}}\right\} \]$

(5)

점성 감쇠력은 회전축의 위치에 따라 변하므로 회전축의 각도를 0도부터 360도까지 60도 간격으로 6개의 회전축 각도(0도, 60도, 120도, 180도, 240도, 300도)에 대하여 각각 로터의 자유감쇠 수치실험을 실시하고, 그 결과를 곡선 맞춤하였다. Fig. 3은 OpenFOAM 계산에 사용된 수평원기둥 로터의 주변 격자 체계를 보여주고 있다. 계산의 효율성을 높이기 위하여 로터 주변과 자유수면 근처의 격자의 간격을 더 줄였다(Ko et al.[2018a]; Ko et al.[2018b]). Fig. 4는 특정 회전축 각도에서 로터의 자유감쇠 수치실험의 시계열 횡 운동 변위를 보여주고 있다. 회전축의 각도에 따라 로터의 감쇠 정도가 달라짐을 확인할 수 있다.

Fig. 3.

2D grid system near WEC and free surface.

Fig. 4.

Time histories of roll motion with different COR position (R=0.75r).

Fig. 5는 자유감쇠 수치실험을 통하여 얻은 6개의 회전축 각도에서 얻은 감쇠비를 가지고 5차 다항식을 이용한 곡선 맞춤을 실시한 결과이며, 곡선 맞춤식은 식 (6)과 같다.

Fig. 5.

Curve fitting of damping ratio.

$\[ \kappa =c_5{\alpha }^5+c_4{\alpha }^4+c_3{\alpha }^3+c_2{\alpha }^2+c_1\alpha +c_0 \]$

(6)

여기서

c₅ = -3.31 × 10^-13, c₄ = 5.38 × 10^-10,
c₃ = -2.58 × 10^-7, c₂ = 4.69 × 10^-5,
c₁ = -3.01 × 10^-3, c₀ = 0.153.

2.2.4 정유체 모멘트

횡 방향 정유체 모멘트 계수는 선형이론 사용 시 상수이지만 비선형 정유체 모멘트인 경우에는 로터의 횡 운동에 따른 잠긴 형상(Submerged shape)의 변화와 로터의 횡 운동에 의해 회전축으로부터 로터의 무게중심과 부력중심이 이동하면서 발생하는 모멘트암(Moment arm)의 변화를 고려해야 한다. 식 (7)의 오른쪽에서 첫번째 항이 로터의 잠긴 형상 변화를, 두번째 항이 모멘트암의 영향을 나타낸다.

$\[ K_{55}=\rho g{\iint }_{A_{WP}}{x}^{2}ds+\rho gV\left(z_B-z_G\right) \]$

(7)

여기서 ρ는 해수의 밀도, A_WP는 수선면적, V는 잠긴 체적, z_G와 z_B는 각각 무게와 부력 중심의 z좌표이다.

일반적으로 운동 변위가 작은 부유체인 경우 정유체 모멘트는 상수인 복원력 계수와 횡 운동 변위의 곱으로 표현된다. 그러나 현재 해석 모델인 수평원기둥 로터는 수면 위에 노출된 편심된 회전축을 가지고 있어 작은 횡 운동 변위에도 잠긴 형상의 변화가 크기 때문에 선형 정유체 모멘트로 가정하여 해석할 경우, 로터의 횡 운동을 정확하게 예측할 수 없다. 따라서 본 연구에서는 정유체 모멘트를 Table 2와 같이 3가지 모델로 구분하여 시간영역 해석에 각각 적용하였다. 여기서 부분 비선형모델은 횡 운동에 따른 모멘트암의 변화로 발생하는 비선형 모멘트를 고려한 모델이다. 완전 비선형 모델은 매시간 변하는 로터의 횡 운동 변위에 따른 로터의 잠긴 부분의 형상과 모멘트암의 변화를 모두 반영하여 계산된 비선형 모멘트를 나타낸 모델이다.

Table 2.

Hydrostatic moment model

3. 시간영역 운동해석 결과

식 (1)에 주어진 시간영역 운동방정식을 풀기 위하여 상용 프로그램인 OrcaFlex를 사용하였다. OrcaFlex는 로터의 관성모멘트, 정유체 모멘트, 부력 중심, 주파수영역 해석을 통해 얻은 횡 방향 동유체 모멘트, 파기진 모멘트를 입력 값으로 하여 로터의 횡 방향 운동 변위, 속도, 가속도를 계산한다. 주파수영역 해석은 경계요소법에 기반을 둔 상용코드인 WAMIT을 이용하였다.

로터의 초기 과도 응답을 줄이기 위하여 램프 함수를 이용하여 400초 동안 점진적으로 외력을 증가시켰다. 식 (1)을 풀어 구한 로터의 횡 운동 속도 시계열 값을 가지고 시간평균 추출파워(P_m)를 구하면 식 (8)과 같다.

여기서 T는 총 발전 시간이다. 본 계산에서는 총 발전 시간을 30분으로 하여 시간평균 추출파워를 구하였다.

3.1 회전축 위치

수평원기둥 로터의 최적의 회전축의 위치를 정하기 위하여 설치 해역의 대표 파랑스펙트럼을 대상으로 부분 비선형 모델을 가지고 시간영역 운동해석을 실시한 뒤, 높은 추출파워를 보이는 회전축 각도 범위를 선정하고, 이 각도 범위내에서 완전 비선형 모델을 적용하여 가장 높은 추출파워를 보이는 최적의 회전축의 각도를 정하였다. Fig. 6은 설치 대상 해역인 제주도 차귀도 인근 해상의 대표 파랑스펙트럼인 유의파고 2.0m, 피크주기 6.65 s, γ = 1.3인 JONSWAP 파랑스펙트럼을 보여주고 있다. Fig. 7은 주어진 JONSWAP 파랑스펙트럼을 이용하여 생성된 불규칙파의 시계열 데이터이다.

Fig. 6.

JONSWAP wave spectrum with HS = 2.0m, TP = 6.65 s and γ = 1.3.

Fig. 7.

Wave elevation generated by given JONSWAP wave spectrum.

Fig. 8은 선형과 부분 비선형 모델을 적용하여 로터의 시간영역 운동해석을 실시하였을 때의 로터의 추출파워 결과이다. 전반적으로 부분 비선형 모델을 적용했을 때의 해석 결과가 선형 모델을 적용했을 때의 결과보다 다소 작게 나타났지만 그 차이는 매우 작았다. 즉, 파랑중 로터가 횡 운동하면서 나타나는 모멘트암에 의한 복원모멘트의 변화는 무시할 수 있다. 추출파워는 선형 모델에서는 회전축이 320도일 때 10.37 kW, 310도에 위치하였을 때 10.32 kW로 가장 높게 나타났으며, 부분 비선형 모델에서도 320도에 위치했을 때 10.25 kW, 310도에 위치했을 때 10.11 kW로 가장 높게 나타났다. 두 모델의 결과를 종합하면 로터의 회전축이 300도~340도 사이에 위치했을 때, 높은 추출파워를 보임을 확인하였다. 장치의 설치 및 유지보수 측면에서도 회전축이 수면 가까이 위치하는 것보다 로터의 상부 위치인 이 구간에 위치하는 것이 바람직하다.

Fig. 8.

Extracted power with linear (solid line with circle) and partially nonlinear (dashed line with triangle) hydrostatic stiffness.

완전 비선형 모델을 적용한 시간영역 운동해석은 선형과 부분 비선형 모델 적용하였을 때 높은 추출파워를 보이는 회전축 각도 300도~340도 범위를 대상으로 이뤄졌다. 완전 비선형 모델을 적용하기 위하여 각 회전축의 위치에서 허용 운동 범위인 ±60도 내에서 10도 간격으로 복원모멘트를 WAMIT을 사용하여 미리 구한 뒤 Fig. 9와 같이 곡선 맞춤을 통하여 곡선 맞춤식을 구하고 이를 이용하여 실시간으로 변하는 로터의 횡 운동 변위에 맞는 복원모멘트를 찾아 시간영역 운동방정식에 대입하여 횡 운동 변위를 계산하였다.

Fig. 9.

Roll displacement–hydrostatic moment curve according to the position of COR.

Fig. 10은 회전축의 각도가 330도에 위치한 로터의 횡 운동 해석 시 정유체 모멘트를 부분 비선형 모델과 완전 비선형 모델로 각 각 적용했을 때의 횡 운동 변위와 파워의 시계열 결과를 비교한 것이다. Fig. 10(a)의 횡 운동 변위 결과를 보면 부분 비선형 모델을 적용하였을 때 로터의 횡 운동이 2차 변환장치의 안정성을 고려한 횡 운동 허용 구간인 +60도~-60도를 넘어가지는 않으나, 완전 비선형 모델을 적용하였을 때 횡 운동의 변위가 상대적으로 크게 일어나 일부 시간에서는 허용 운동 범위를 넘어가는 것을 볼 수 있다. 추출파워 결과도 완전 비선형 모델을 적용하였을 때 더 큰 값을 보이고 있다(Fig. 10(b)).

Fig. 10.

Time histories of rotor simulation in the case of α = 330 deg.

Table 3은 우수한 추출효율을 보이는 5개의 회전축의 각도에 대하여 부분 비선형 모델과 완전 비선형 모델을 각각 적용하였을 때 계산된 시간평균 추출파워를 비교한 결과이다. 완전 비선형 모델을 적용하였을 때, 가장 높은 추출파워를 보이는 회전축의 각도는 330도로 부분 비선형 모델을 적용했을 때 나타난 320도와는 다른 결과를 보였다. 부분 비선형 모델을 적용했을 때와 완전 비선형 모델을 적용했을 때의 정량적인 크기도 차이를 보였다. 완전 비선형 모델을 적용한 추출파워가 약 40% 높게 나타났다. 본 수평원기둥 로터와 같이 회전축이 수면보다 위에 위치하고, 회전 운동을 하면서 잠긴 형상의 변화가 큰 모델인 경우에는 완전 비선형 모델을 사용해야 보다 정확한 결과를 얻을 수 있음을 알 수 있었다. 이상의 결과를 종합하면 흘수가 1.6 m인 로터의 회전축 위치가 R = 1.5m과 α = 330도일 때 가장 높은 추출파워를 보였다. 이후 모든 계산에서 이 최적값을 사용하였다.

Table 3.

Comparison of the time-averaged extracted power

3.2 연간 에너지 생산량

수평원기둥 로터를 이용한 파력발전장치의 연간 에너지 생산량 계산을 위하여 로터의 최적의 위치에 회전축을 고정시킨 후에 로터의 횡 운동 변위와 속도를 구하고 이로부터 추출파워를 계산하였다. 개발된 해석 프로그램을 제주도 차귀도 해역에서 얻은 파랑산포도에 적용하여 파워 매트릭스를 계산하고 이로부터 연간 에너지 생산량을 추정하였다. Table 4는 연간 에너지 생산량을 계산하기에 앞서 입력 자료인 제주도 차귀도 파랑 산포도를 보여주고 있다. 이 파랑 통계자료는 제주도 차귀도 해역 근방에서 1979년 1월 1일부터 2008년 9월 1일까지 관측과 수치계산을 통하여 얻은 자료이다. 파랑 산포도 내부에는 발생 빈도수가 백분율로 표시되어 있으며, 편의상 소수점 둘째 자리까지 표시하였다. 비어 있는 셀들은 측정 기간동안 한번도 발생하지 않는 파랑 조건들이다. 입사파의 파향 변화는 고려하지 않았으며, 모든 파는 로터의 정면에 수직방향으로 입사한다고 가정하였다.

Table 4.

Wave scatter diagram (Measured from 1979.01.01.10:00 to 2008.09.01.09:00)

Table 5는 파랑 산포도의 각 셀에 해당하는 유의파고와 피크주기에 대하여 로터의 횡 운동 시간영역 해석을 통해 계산된 기계파워(Mechanical power, P_mi)를 Table 4의 빈도수(f_i)와 곱하여 얻은 기계파워 매트릭스이다.

Table 5.

Mechanical power matrix (kW, Pmi × fi)

Table 5의 기계파워 매트릭스의 결과를 모두 합하여 식 (10)과 같이 2차 변환효율(η₁ = 0.8)을 곱하면 시간평균 전력파워(P_ave)를 구할 수 있다. 계산된 로터의 시간평균 전력파워는 약 3.46 kW이다.

$\[ {\overline{P}}_{ave}=\sum _{i=1}^N{\overline{P}}_{e_i}\times f_i={\eta }_1\times \sum _{i=1}^N{\overline{P}}_{m_i}\times f_i \]$

(10)

계산된 시간평균 전력파워를 가지고 연간 에너지 생산량을 식 (11)과 같이 계산하였다.

$\[ AEP={\overline{P}}_{ave}\times T_h\times {\eta }_2\times {\eta }_3 \]$

(11)

여기서 T_h는 8760시간(1년)이며, η₂는 운용효율(Operating rate)로 0.95, η₃은 송전효율(Transmission rate)로 0.98을 사용하였다(Neary et al.[2014]). 회전축이 R = 1.5 m, α = 330도에 위치한 수평원기둥 로터를 이용하여 발전되는 연간 에너지 생산량은 28.26 MWh이다.

4. 결 론

본 연구에서는 수평원기둥 로터의 추출파워에 큰 영향을 미치는 편심 된 회전축의 최적 위치를 정하고, 최적 회전축 위치에서 로터의 추출파워를 계산하고, 설치 대상 해역인 제주도 차귀도 해역에서 얻은 파랑 산포도에 적용하여 최종적으로 연간 에너지 생산량을 추정하였다. 회전축의 최적 위치를 선정하기 위하여 설치 대상 해역의 대표 파랑 스펙트럼(유의파고 2.0 m, 피크주기 6.55s)을 이용하여 회전축을 0도에서 360도까지 10도 간격으로 각도를 변화시켜가며 추출파워를 계산하였다. 이때 2차 변환장치인 유압장치의 안정성을 고려하여 로터의 운동을 -60도에서 +60도까지 제한하였다. 시간영역 해석 시, 복원모멘트를 선형 모델, 부분 비선형 모델, 완전 비선형 모델로 나눠 해석하였고 각각에 대하여 추출파워의 결과를 비교하였다. 선형 모델과 부분 비선형 모델의 추출파워의 결과는 큰 차이가 없었으며 가장 높은 추출파워를 보이는 회전축 각도는 320도였다. 그러나 완전 비선형 모델을 적용한 경우에는 다른 두 모델의 결과와 다소 다른 결과를 보였다. 가장 높은 추출파워를 보이는 회전축 각도는 330도였으며 선형과 부분 비선형 모델보다 큰 추출파워를 보였다. 이는 수면 위에 위치한 회전축을 중심으로 로터가 운동하면서 횡 운동에 따른 잠긴 형상의 변화가 크게 나타났기 때문으로 사료된다. 완전 비선형 모델 적용을 통하여 얻은 회전축의 최적 각도 330도를 연간 에너지 생산량을 추정할 때 사용하였다. 설치 대상 해역의 파랑 산포도의 각 셀의 파랑조건들에 대하여 추출파워를 구하여 로터의 기계 파워 매트릭스를 구하고 이로부터 시간평균 전력파워를 계산하였다. 최종적으로 2차 변환장치 효율, 운용 효율, 송전 효율을 모두 고려하여 얻은 연간 에너지생산량은 28.26 MWh로 나타났다.

본 연구에서는 입사파의 진행방향은 로터에 수직 방향이고 로터의 폭이 로터의 직경에 비하여 긴 2차원 문제로 가정하여 간단하게 해석하였다. 향후 입사파의 진행방향을 고려한 다방향 불규칙파에 로터가 놓여 있거나 다수의 로터가 여러 줄로 배열 된 경우에는 3차원 수치 해석모델 개발이 불가피하며 이에 대한 연구가 필요하다. 본 연구에서는 연간에너지생산량을 추정할 때, 일방향 불규칙 입사파로 가정하여 계산된 추출파워를 이용하였으나 실제 해양환경과 같이 다방향 불규칙 입사파를 고려하여 계산된 추출파워를 이용해야 한다. 또한 2차 변환장치의 종류에 따른 2차 변환효율, 로터와 회전축 사이의 기계적 마찰에 의한 감쇠 등을 반영한 해석이 진행되어야 보다 실제적인 연간 에너지 생산량을 구할 수 있다. 이 밖에 송전 효율 및 정격용량과 관련이 있는 운용 효율에 대한 보다 구체적인 값도 필요하다.

Acknowledgments

본 연구는 산업통상자원부(MOTIE)와 한국에너지기술평가원(KETEP)의 지원을 받아 수행한 연구 과제입니다(No. 20163010071690).

References

• Cho, I.H., Koh, C.H., and Bae, Y.H., (2018a), Performance Analysis of a Horizontal Cylinder Wave Energy Converter with Off-centered Rotational Axis, J. Korean Soc. Mar. Environ. Energy, 21(1), p10-22. [https://doi.org/10.7846/jkosmee.2018.21.1.10]
• Cho, I.H., Kim, J.R., and Bae, Y.H., (2018b), Time-domain Analysis of Horizontal Cylinder Wave Energy Converter with Off-centered Rotational Axis, J. Korean Soc. Mar. Environ. Energy, 21(2), p97-106. [https://doi.org/10.7846/jkosmee.2018.21.2.97]
• Cummins, W.E., (1962), The Impulse Response Function and Ship Motions, International Symposium on Ship Theory, Hamburg, Germany, p8.
• Journée, J.M.J., and Massie, W.W., (2001), Offshore hydromechanics, first ed., Delft University of Technology, The Netherlands.
• Ko, H.S., Cho, I.H., and Bae, Y.H., (2018a), Non-Linear Dynamic Analysis of an Eccentrical Cylinder-Type Wave Energy Converter using OpenFOAM, J. Comput. Fluids Eng., 23(3), p18-26. [https://doi.org/10.6112/kscfe.2018.23.3.018]
• Ko, H.S., Bae, Y.H., and Cho, I.H., (2018b), Dynamic Analysis of a Floating Structure Using OpenFOAM, J. Comput. Fluids Eng., 23(1), p101-112. [https://doi.org/10.6112/kscfe.2018.23.1.101]
• Lucas, J., Salter, S.H., Cruz, J., Taylor, J., and Bryden, I., (2009), Performance optimization of a modified duck through optimal mass distribution, in proc. of the 8th European Wave and Tidal Energy Conference, Uppsala, Sweden, p10.
• Neary, V.S., Previsic, M., Jepsen, R.A., Lawson, M.J., Yu, Y.-H., Copping, A.E., Fontaine, A.A., and Hallett, K.C., (2014), Methodology for Design and Economic Analysis of Marine Energy Conversion (MEC) Technologies, Sandia National Laboratories, Albuquerque, NM.
• Salter, S.H., (1974), Wave Power, Nature, 249, p720-724.
• Salter, S.H., Jeffrey, D.C., and Taylor, J., (1975, First year interim report on Edinburgh wave power project: study of mechanism to extract power from sea waves, Technical report, Wave-power project University of Edinburgh, September).