태양광 발전 시스템의 효율을 극대화하기 위해서는 지표면의 고도, 경사도(Tilt) 및 방위각(Azimuth)에 따른 태양 복사 에너지의 변화를 정확하게 예측하는 것이 필수적입니다. 본 문서에서는 MATLAB을 활용하여 임의의 평면에서 수신되는 태양 복사량을 계산하는 물리적 모델의 구현 방법을 설명합니다.

1. 태양 복사 모델의 물리적 구성

경사면에서 수신되는 전체 일사량(\(I_{total}\))은 크게 세 가지 성분으로 구분됩니다.

• 직달 일사량(Direct Beam Radiation, \(I_b\)): 태양으로부터 가려짐 없이 직접 지표면에 도달하는 성분입니다.
• 산란 일사량(Diffuse Sky Radiation, \(I_d\)): 대기 중의 입자에 의해 산란되어 도달하는 성분입니다.
• 반사 일사량(Ground Reflected Radiation, \(I_r\)): 주변 지면에서 반사되어 경사면에 도달하는 성분입니다.

2. 태양 위치 및 좌표 계산 알고리즘

먼저 특정 위도와 경도, 시간에서 태양의 고도각(Altitude)과 방위각(Azimuth)을 산출해야 합니다. 아래 함수는 적위(Declination)와 시각(Hour Angle)을 기반으로 태양 위치를 계산합니다.

function [sunElev, sunAzim] = getSunPosition(lat, lon, targetTime)
    % lat: 위도, lon: 경도, targetTime: datetime 객체
    
    dayOfYear = day(targetTime, 'dayofyear');
    % 태양 적위 계산 (Cooper의 식)
    declination = 23.45 * sin(deg2rad(360 * (284 + dayOfYear) / 365));
    
    % 지방시 및 시각 계산
    localHour = hour(targetTime) + minute(targetTime)/60;
    timeOffset = (lon - 120) * 4 / 60; % 기준 자오선(한국/중국 등) 기준 보정
    solarTime = localHour + timeOffset;
    hourAngle = 15 * (solarTime - 12);
    
    % 고도각 계산
    sinElev = sin(deg2rad(lat)) * sin(deg2rad(declination)) + ...
              cos(deg2rad(lat)) * cos(deg2rad(declination)) * cos(deg2rad(hourAngle));
    sunElev = rad2deg(asin(sinElev));
    
    % 방위각 계산
    cosAzim = (sin(deg2rad(declination)) - sinElev * sin(deg2rad(lat))) / ...
              (cos(deg2rad(sunElev)) * cos(deg2rad(lat)));
    sunAzim = rad2deg(acos(max(min(cosAzim, 1), -1)));
    
    if hourAngle > 0
        sunAzim = 360 - sunAzim;
    end
end

3. 경사면 복사량 산출 모듈

임의의 경사각(\(\beta\))과 면의 방위각(\(\alpha\))을 가진 평면에 도달하는 복사 강도를 계산합니다. 여기서는 Hottel의 맑은 하늘 모델을 기반으로 대기 투과율을 고려합니다.

function totalIrradiance = computeSlopeIrradiance(lat, lon, timestamp, tilt, surfaceAzim, Gsc)
    % Gsc: 태양 상수 (약 1367 W/m^2)
    
    [alt, azim] = getSunPosition(lat, lon, timestamp);
    
    if alt <= 0
        totalIrradiance = 0;
        return;
    end
    
    % 대기 투과율 모델 (간략화된 Hottel 모델)
    transmissivity = 0.7^(1/sin(deg2rad(alt)));
    
    % 입사각(Angle of Incidence) 계산
    cosTheta = sin(deg2rad(alt)) * cos(deg2rad(tilt)) + ...
               cos(deg2rad(alt)) * sin(deg2rad(tilt)) * cos(deg2rad(azim - surfaceAzim));
    
    % 1. 직달 성분
    Ib = Gsc * transmissivity * max(0, cosTheta);
    
    % 2. 산란 성분 (Isotropic Sky 모델)
    Id = Gsc * 0.3 * (1 - transmissivity) * sin(deg2rad(alt)) * (1 + cos(deg2rad(tilt))) / 2;
    
    % 3. 지면 반사 성분 (반사율 0.2 가정)
    Ir = Gsc * sin(deg2rad(alt)) * 0.2 * (1 - cos(deg2rad(tilt))) / 2;
    
    totalIrradiance = Ib + Id + Ir;
end

4. 응용 예제: 최적 경사각 탐색

특정 위치에서 특정 날짜의 누적 에너지 수신량을 최대화하는 최적 경사각을 찾는 시뮬레이션입니다.

% 설정 파라미터
latitude = 37.56; % 서울 기준
longitude = 126.97;
targetDate = datetime(2025, 5, 20);
solarConstant = 1367;

tilts = 0:1:90;
energyYield = zeros(size(tilts));

for i = 1:length(tilts)
    hourlyEnergy = 0;
    % 일출부터 일몰까지 10분 간격 적분
    for h = 6:1/6:18
        currentTime = targetDate + hours(h);
        val = computeSlopeIrradiance(latitude, longitude, currentTime, tilts(i), 180, solarConstant);
        hourlyEnergy = hourlyEnergy + val * (1/6);
    end
    energyYield(i) = hourlyEnergy;
end

[maxEnergy, idx] = max(energyYield);
fprintf('최적 경사각: %d도, 예상 일일 복사량: %.2f Wh/m^2\n', tilts(idx), maxEnergy);

plot(tilts, energyYield, 'LineWidth', 2);
grid on;
xlabel('Tilt Angle (deg)');
ylabel('Daily Energy (Wh/m^2)');
title('경사각에 따른 일일 수신 에너지 변화');

5. 지형 차폐 및 고도 데이터 활용

복잡한 산악 지형에서는 주변 산의 그림자에 의한 차폐 효과를 고려해야 합니다. DEM(Digital Elevation Model) 데이터를 사용하면 각 격자점에서의 가시선(Line of Sight) 분석을 통해 그림자 여부를 판단할 수 있습니다. viewshed 함수나 벡터의 내적을 활용하여 태양 벡터가 지형 선에 의해 가려지는지 체크함으로써 모델의 정밀도를 높일 수 있습니다.

또한 대규모 연산을 수행할 경우 MATLAB의 parfor를 이용한 병렬 처리나 GPU 배열을 활용하여 수만 개의 격자점에 대한 복사량 분포 지도를 빠르게 생성할 수 있습니다.