مدل سازی تابش خورشید | پیاده سازی در متلب

انرژی خورشیدی بخشی از گرما و نور حاصل از تابش خورشید است که برای کاربردهای مختلف در سطح زمین در دسترس است. از این کاربردها می‌توان به برانگیختن الکترون‌های یک سلول فتوولتائیک، فراهم کردن انرژی برای فرایندهای طبیعی مانند فتوسنتز یا گرم کردن اشیاء اشاره کرد. این انرژی مجانی، پاک و فراوان در طول سال در بیشتر نقاط در دسترس است و به دلیل هزینه سوخت‌های فسیلی و آلودگی هوای ناشی از آن‌ها مورد توجه قرار گرفته است. انرژی موجود در تابش خورشید شامل دو قسمت است: تابش فرازمینی (Extraterrestrial Solar Radiation) که بالاتر از جوّ است و تابش کلی (Global Solar Radiation) که در سطح زیر جوّ قرار دارد.

تابش کلی معمولاً توسط آفتاب‌سنج‌ها (پیلومومترها)، حسگرهای خورشیدی، تابه‌سنج‌ها (اکتینوگراف‌ها) یا گرماسنج‌ها اندازه‌گیری می‌شود. این ابزارهای اندازه‌گیری معمولاً در محل‌های خاصی از مناطق مشخص نصب می‌شوند. با توجه به هزینه بالای این دستگاه‌ها، امکان نصب آن‌ها در بسیاری از محل‌ها وجود ندارد. از سوی دیگر، این دستگاه‌های اندازه‌گیری، تلرانس قابل توجه و دقت کمی دارند؛ در نتیجه ممکن است موجب ثبت اشتباه داده‌ها یا از دست رفتن آن‌ها شوند. مدل‌سازی منبع خورشیدی با استفاده از علم نجوم و اصول هندسه راه‌حل مناسبی برای این رفع مشکل است. همچنین، می‌توان از مقدار تابش اندازه‌گیری‌ شده، برای تشکیل مدل‌های تابش بهره برد. این مدل‌ها روابط ریاضی بین تابش خورشید و متغیرهای هواشناسی مانند دمای محیط، رطوبت و نرخ تابش را توصیف می‌کنند. همچنین می‌توان از این مدل‌ها برای پیش‌بینی تابش در مکان‌هایی استفاده کرد که دستگاه اندازه‌گیری انرژی خورشیدی نصب نشده است.

مدل‌سازی موقعیت خورشید

همان‌طورکه می‌دانیم، زمین در یک مدار بیضوی اطراف خورشید گردش می‌کند. شکل ۱، مدار چرخش زمین به دور خورشید را نشان می‌دهد.

فیلم آموزش فیزیک پایه 1 – جامع و کاربردی در فرادرس

کلیک کنید

مدت زمان هر چرخش زمین به دور خورشید حدود 8766 ساعت، تقریباً برابر با 365٫242 روز است.

نقاط منحصر به‌ فردی در مدار شکل ۱ وجود دارد. تحویل زمستان در 21 دسامبر رخ می‌دهد که در آن، زمین تقریباً 147 میلیون کیلومتر از خورشید فاصله دارد. از سوی دیگر، روز تحویل تابستان که در 21 ژوئن اتفاق می‌افتد، زمین حدود 152 میلیون کیلومتر از خورشید دور است. برای بررسی دقیق‌تر می‌توان گفت زمین در 2 ژانویه به خورشید نزدیک است (147 میلیون کیلومتر) که این نقطه به‌نام «حضیض» (Perihelion) شناخته می‌شود. نقطه‌ای که زمین از خورشید دورتر است (152 میلیون کیلومتر)، «اوج» (Aphelion) نامیده می‌شود و در 3 جولای رخ می‌دهد.

همان‌طورکه در شکل ۲ دیده می‌شود، برای ناظری که در نقطه مشخصی از زمین ایستاده است، موقعیت خورشید با دو زاویه اصلی یعنی زاویه فراز (Altitude angle) $$\alpha$$ و زاویه سمت (Azimuth angle) $$\theta _\mathrm{S}$$ تعیین می‌شود.

در شکل 2، زاویه فراز، زاویه ارتفاع خورشید در آسمان است که نسبت به افق اندازه‌گیری می‌شود. زاویه فراز را می‌توان با

$$\large \begin {align} \mathrm {sin} \, \alpha \, = \mathrm {sin} \, L \, \mathrm {sin} \delta + \mathrm {cos} \, L \, \mathrm {cos} \, \delta \, \mathrm {cos} \, \omega \end {align} \;\;\;\;\; ( 1 )$$

بیان کرد که در آن، $$L$$ «عرض جغرافیایی» (Latitude) محل، $$\delta$$ «زاویه میل» (Angle of declination) و $$\omega$$ «زاویه ساعت» (Hour angle) است.

زاویه میل، زاویه بین بردار زمین-خورشید و صفحه استوا است (شکل ۳ را ببینید) و به‌ صورت زیر محاسبه می‌شود (نتایج به درجه و آرگومان توابع مثلثاتی برحسب رادیان است):

$$\large \begin {align} \sigma _ \mathrm {S} = 2 3 . 4 5 ^ { \circ } \mathrm {sin} \left [ \frac { 2 \pi ( N - 8 1 ) } { 3 6 5 } \right ] \end {align} \;\;\;\;\; ( 2 )$$

زاویه ساعت $$\omega$$، جابه‌جایی زاویه‌ای خورشید از نقطه محلی است و با

$$\large \begin{align} \omega = 15 ^ { \circ } ( \mathrm {AST} - 1 2 \, \mathrm {h} ) \end {align} \;\;\;\;\; ( 3 )$$

به‌ دست می‌آید که در آن، AST «زمان ظاهری یا حقیقی خورشیدی» (True Solar Time) است و با حرکت ظاهری روزانه خورشید تعیین می‌شود. زمان ظاهری بر اساس روز ظاهری محاسبه می‌شود و فاصله زمانی بین دو بازگشت متوالی خورشید به نصف‌النهار محلی است. زمان ظاهری با

$$\large \begin{align} \mathrm {AST} = \mathrm {LMT} + \mathrm {EoT} \pm 4 ^ \circ / ( \mathrm {LSMT} - \mathrm {LOD} ) \end {align} \;\;\;\;\; ( 4 )$$

تعیین می‌شود که در آن، LMT «زمان نصف‌النهار محلی» (Local Meridian Time) و LOD «طول جغرافیایی» (Longitude) و LSMT «زمان نصف‌النهار استاندارد محلی» (Local Standard Meridian Time) و EOT یک «معادله زمانی» است.

LSMT نصف‌النهار مرجع است که برای یک منطقه زمانی خاص استفاده می‌شود و مشابه نصف‌النهار اولیه‌ای است که برای زمان میانگین گرینویچ به‌ کار می‌رود. LSMT را می‌توان با رابطه زیر محاسبه کرد:

$$\large \begin {align} \mathrm {LMST} = 1 5 ^ \circ T _ \mathrm {GMT} \end {align} \;\;\;\;\; ( 5 )$$

EoT اختلاف بین زمان‌های میانگین و ظاهری خورشید است که هر دو در طول جغرافیایی خاصی در یک لحظه واقعی رخ می‌دهند. EoT توسط

$$\large \begin {align} \mathrm {EoT} = 9 . 8 7 \mathrm {sin} ( 2 B ) - 7 . 5 3 \mathrm {cos} B - 1 . 5 \mathrm {sin} B \end {align} \;\;\;\;\; ( 6 )$$

مشخص می‌شود و در آن می‌توان $$B$$ را با رابطه زیر محاسبه کرد:

$$\large \begin {align} B = \dfrac { 2 \pi } { 3 6 5 } ( N - 8 1 ) \end {align} \;\;\;\;\; ( 7 )$$

از طرف دیگر، زاویه سمت که در شکل ۲ دیده می‌شود، جابه‌جایی زاویه‌ای خط مرجع خورشید از محور منبع است. زاویه سمت را می‌توان با استفاده از رابطه زیر به‌ دست آورد:

$$\large \begin {align} \mathrm {sin} \theta = \dfrac { \mathrm { \delta } \mathrm {sin} \omega } { \mathrm {cos} \alpha } \end {align} \;\;\;\;\; ( 8 )$$

مثال 1

برنامه‌ای در متلب بنویسید که زاویه‌های فراز و سمت را در ساعت 13:12 روز دوم جولای برای شهر کوالالامپور محاسبه کند.

حل: بخش‌های اصلی برنامه به‌ صورت زیر است:

• مختصات مکان (طول و عرض جغرافیایی)، شماره روز و زمان میانگین محلی را وارد کنید.
• زاویه میل، معادله زمانی و LMST را محاسبه کنید.
• AST و زاویه ساعت را محاسبه کنید.
• زاویه فراز را به‌دست آورید.
• زاویه سمت را محاسبه کنید.
• نتایج را رسم کنید.

کد متلب این مثال به صورت زیر است:

1%Example 1.1
2%----------------------------------------------------------
3%date 2/7/2015 (N=183)
4%location Kuala Lumpur, Malaysia, L =(3.12), LOD = (101.7)
5L=3.12; %Latitude
6LOD=101.7; %longitude
7N=183; %Day Number
8T_GMT=8; %Time difference with reference to GMT
9LMT_minutes=792; %LMT in minutes
10Ds=23.45*sin((360*(N-81)/365)*(pi/180)); % angle of
11declination
12B=(360*(N-81))/364; %Equation of Time
13EoT=(9.87*sin(2*B*pi/180))- (7.53*cos(B*pi/180))-
14(1.5*sin(B*pi/180)); % Equation of Time
15Lzt= 15* T_GMT; %LMST
16if LOD>=0
17Ts_correction= (-4*(Lzt-LOD))+EoT; %solar time correction
18else
19Ts_correction= (4*(Lzt-LOD))+EoT; %solar time correction
20end
21Ts= LMT_minutes + Ts_correction; %solar time
22Hs=(15 *(Ts - (12*60)))/60; %Hour angle degree
23s i n _ A l p h a = ( s i n ( L * p i / 1 8 0 ) * s i n ( D s * p i / 1 8 0 ) ) +
24(cos(L*pi/180)*cos(Ds*pi/180)* cos(Hs*pi/180)); %altitude
25angle
26Alpha=asind(sin_Alpha) %altitude angle
27Sin_Theta= (cos (Ds*pi/180)*sin (Hs*pi/180))./cos(Alpha_i.
28*pi/180); %Azimuth angle
29Theta=asind(Sin_Theta) %Azimuth angle

پاسخ: $$\mathrm{Alpha}=70.04^\circ$$؛ $$\mathrm{theta}=-1.13^\circ$$.

مثال 2

برنامه نوشته‌ شده مثال 1 را برای محاسبه نمودار زاویه‌های فراز و سمت (هر 5 دقیقه) کل روز خورشیدی دوم جولای شهر کوالالامپور اصلاح کنید.

حل: روز خورشیدی به‌عنوان طول مدت طلوع تا غروب خورشید تعریف می‌شود. بنابراین، باید زاویه فراز و زاویه سمت را برای هر ساعت از طلوع تا غروب خورشید محاسبه کرد. می‌توان زاویه‌های طلوع و غروب خورشید را برابر در نظر گرفت و به‌ صورت زیر به‌ دست آورد:

$$\large \begin{align} \omega _ \mathrm {ss,sr} = \mathrm {cos} ^ { - 1 } ( - \mathrm {tan} L \mathrm {tan} \delta ) . \end {align} \;\;\;\;\; ( 9 )$$

در عین حال، زمان خورشید هر زاویه ساعت را می‌توان با بازنویسی معادله (۳)‌ به‌ صورت زیر محاسبه کرد:

$$\large \begin{align} \dfrac { \omega _ { \mathrm {ss,sr} } } { 1 5 ^ \circ } \pm 1 2 \mathrm { h } = \mathrm {AST} _ \mathrm {sr,ss} . \end {align} \;\;\;\;\; ( 10)$$

اگر بخواهیم زمان طلوع خورشید را محاسبه کنیم، علامت معادله (۱۰) باید منفی باشد و در صورتی ‌که قصد محاسبه زمان غروب آفتاب را داشته باشیم باید آن را مثبت قرار دهیم. پس از آن، بخش اصلی ساختار برنامه را می‌توان به‌ شرح زیر نوشت:

• مختصات مکان (طول و عرض جغرافیایی) و شماره روز را وارد کنید.
• زاویه میل را حساب کنید.
• زاویه‌های ساعت طلوع و غروب خورشید را محاسبه کنید.
• AST طلوع و غروب را محاسبه کنید.
• معادله زمانی و LMST را محاسبه کنید.
• زمان‌های واقعی طلوع و غروب را به‌دست آورید.
• یک حلقه با شروع از طلوع و پایان غروب خورشید با اندازه گام 5 دقیقه تعریف کنید.
• زمان خورشید و زاویه ساعت را در هر گام محاسبه کنید.
• زاویه فراز را در هر گام محاسبه کنید.
• زاویه سمت را در هر گام محاسبه کنید.
• زاویه‌های فراز و سمت محاسبه‌ شده را در آرایه‌ها ذخیره کنید.
• نتایج را رسم کنید.

برنامه متلب این مثال به صورت زیر است:

1%Example 1.2
2%---------------------------------------------------------
3%Date 2/7/2015 (N=183)
4%Location Kuala Lumpur, Malaysia, L =(3.12), LOD = (101.7)
5%Actual solar day time 07:11 to 19:22
6L=3.12; %altitude
7LOD=101.7; %longitude
8N=183; %Day Number
9T_GMT=8; %Time difference with reference to GMT
10Step=5;
11Ds=23.45*sin((360*(N-81)/365)*(pi/180)); % angle of
12declination
13Lzt= 15* T_GMT; %LMST
14B=(360*(N-81))/364; %Equation of Time
15EoT=(9.87*sin(2*B*pi/180))- (7.53*cos(B*pi/180))- (1.5*sin
16(B*pi/180)); %Equation of Time
17if LOD>=0
18Ts_correction= (-4*(Lzt-LOD))+EoT; %solar time correction
19else
20Ts_correction= (4*(Lzt-LOD))+EoT; %solar time correction
21end
22Wsr_ssi=- tan(Ds*pi/180)*tan(L*pi/180);%Sunrise/Sunset
23hour angle
24Wsrsr_ss=acosd(Wsr_ssi);% Sunrise/Sunset hour angle
25ASTsr=abs((((Wsrsr_ss/15)-12)*60));%Sunrise solar time
26ASTss=(((Wsrsr_ss/15)+12)*60);%Sunset solar time
27Tsr=ASTsr+abs(Ts_correction); %Sunrise local time
28Tss=ASTss+abs(Ts_correction); %Sunset local time
29Alpha=[];
30Theta=[];
31for LMT=Tsr:Step:Tss %for loop for the day time
32Ts= LMT + Ts_correction; % solar time at each step
33Hs=(15 *(Ts - (12*60)))/60; % Hour angle degree at each
34step
35s i n _ A l p h a = ( s i n ( L * p i / 1 8 0 ) * s i n ( D s * p i / 1 8 0 ) ) +
36(cos(L*pi/180)*cos(Ds*pi/180)* cos(Hs*pi/180)); %altitude
37angle
38Alpha_i=asind(sin_Alpha) ; %altitude angle
39Alpha=[Alpha;Alpha_i];%store altitude angle in array
40Sin_Theta= (cos (Ds*pi/180)*sin (Hs*pi/180))./cos(Alpha_i.*pi/
41180);%Azimuth angle
42Theta_i=asind(Sin_Theta);%Azimuth angle
43Theta=[Theta;Theta_i];% store azimuth angle in array
44end
45Alpha;
46Theta;
47subplot(2,1,1)%plot results
48plot(Alpha)
49subplot(2,1,2)
50plot (Theta, 'red')

پاسخ: شکل ۴.

مدل‌سازی تابش فرازمینی

اولین گام در مدل‌سازی منبع خورشیدی، تخمین تابش خورشید است. همان‌گونه‌ که می‌دانیم، انرژی تابشی هر جسم را می‌توان به‌ صورت تابعی از دمای آن توصیف کرد. راه متداول برای تخمین انرژی تابشی یک جسم، مقایسه آن با یک جسم سیاه است. جسم سیاه به‌ عنوان یک گسیلنده و یک جذب‌ کننده کامل تعریف شده است. با جذب کامل می‌توان تمام انرژی دریافت‌ شده را با هر بازتابی جذب کرد، در حالی‌ که یک گسیلنده کامل، انرژی را بیش از هر جسم دیگری منتشر می‌کند. قانون پلانک، طول‌ موج منتشر شده جسم سیاه را در یک دمای خاص به‌ صورت زیر توصیف می‌کند:

$$\large E _ \lambda = \dfrac { 3 . 4 7 \times 1 0 ^ 8 } { \lambda ^ 5 \left [ \mathrm { e x p } \dfrac { 1 4 4 00 } { \lambda T } \right ] - 1 }$$

که در آن، $$E_\lambda$$ کلِ تابش در واحد سطح با نرخ انتشار جسم سیاه ($$\mathrm{W}/\mu \mathrm{m^2}$$)، $$T$$ دمای مطلق جسم سیاه ($$\mathrm{K}$$) و $$\lambda$$ طول‌ موج ($$\mathrm{\mu} \mathrm{m}$$) است.

مثال ۳

برنامه‌ای بنویسید که قدرت انتشار طیفی یک جسم سیاه 228 کلوین را در محدوده طول‌ موج‌های 1 تا 60 میکرومتر حساب کند. پس از آن قدرت انتشار را برای طول ‌موج 20 و 30 میکرومتر محاسبه کنید.

حل: بخش اول این مثال را می‌توان به‌ سادگی با استفاده از معادله (۱۱) و محاسبه مقدار آن برای محدوده طول ‌موج مورد نظر به‌ صورت زیر حل کرد:

1%Example 1.3
2%---------------------------------------------------------
3T=288;
4E_lamda=[];
5for lamda=1:1:60;
6E_lamda_i=(3.74*10e8)/(lamda^5*(exp(14400/(lamda*T))-1));
7E_lamda=[E_lamda; E_lamda_i];
8end
9E_lamda
10lamda=1:1:60;
11plot(lamda,E_lamda)

برای محاسبه توان انتشار بین طول‌موج 20 و 30 میکرومتر، می‌توان بخش سایه شکل ۵ را به‌ صورت زیر محاسبه کرد:

$$\large \begin {align} \sum \limits _ { \lambda = 2 0 } ^ { \lambda = 3 0 } { E _ \lambda } = \sum \limits _ { \lambda = 2 0 } ^ { \lambda = 3 0 } { \dfrac { 3 . 4 7 \times 1 0 ^ 8 }{ \lambda ^ 5 \left [ \mathrm { e x p } \dfrac { 1 4 4 0 0 } { \lambda T } \right ] - 1 } } \end {align} \;\;\;\;\; ( 1 2 )$$

که پیاده‌سازی آن در متلب به‌ صورت زیر است:

1T=288;
2E_lamda=[];
3forlamda=20:1:30;
4E_lamda_i= (3.74*10e8)/(lamda^5*(exp(14400/
5(lamda*T))-1));
6E_lamda=[E_lamda; E_lamda_i];
7end
8E_lamda;
9Power=sum(E_lamda)

پاسخ: $$704.0801 \, \mathrm{W/m^2}$$.

دمای درون خورشید حدود 15 میلیون کلوین است، در حالی‌ که دمای سطح آن بسیار سردتر بوده و تقریباً 5778 کلوین (5505 درجه سانتی‌گراد) است. بنابراین، تابشی که از سطح خورشید منتشر می‌شود، یک توزیع طیفی مطابق با پیش‌بینی قانون پلانک برای یک جسم سیاه $$5800\, \mathrm{K}$$ است. مساحت کل زیر منحنی جسم سیاه $$1307 \, \mathrm{W}/\mathrm{m^2}$$ تا $$1393 \, \mathrm{W}/\mathrm{m^2}$$ است که مقدار تابش خورشید فقط در خارج از جو زمین است. این مقدار تابش، «ثابت خورشید» (Solar constant) ($$G_\mathrm{o}$$) نامیده می‌شود، گرچه به‌ دلیل مدار بیضوی زمین، قطر آن و تغییر شرایط در فعالیت خورشید، کاملاً ثابت نیست. مقدار پیشنهادی ثابت خورشیدی از سوی بسیاری از محققان، $$1367 \, \mathrm{W/m^2}$$ است.

مقدار تابش در خارج از جو زمین، به‌ دلیل چرخش زمین حول خورشید تغییر می‌کند. از این‌ رو باید فاصله بین خورشید و زمین را در مدل‌سازی تابش فرازمینی در نظر گرفت. بنابراین، تابش فرازمینی ($$G_\mathrm{ex}$$) با رابطه (۱۳)‌ بیان می‌شود:

$$\large \begin{align} G _ \mathrm { e x } = G _ \mathrm { o } \left ( \dfrac { R _ \mathrm { a v } } { R } \right ) ^ 2 \end {align} \;\;\;\;\; ( 13 )$$

که در آن، $$R_\mathrm{av}$$ فاصله میانگین و $$R$$ فاصله لحظه‌ای بین خورشید و زمین است. فاصله لحظه‌ای بین خورشید و زمین به روزِ سال یا شماره روز بستگی دارد. در واقع، تقریب‌های مختلفی برای عامل ($$R_\mathrm{av}/R$$) وجود دارد. یک تقریب توصیه‌ شده به‌ صورت زیر است:

$$\large \begin {align} \left ( \dfrac { R _ \mathrm { a v } } { R } \right ) ^ 2 = \left ( 1 + 0 . 0 3 3 3 \cos \left ( \dfrac { 2 \pi N } { 3 6 5 } \right ) \right ) \end {align} \;\;\;\;\; ( 1 4 )$$

با جایگذاری معادله (۱۴) در معادله (۱۳)‌، می‌توان تابش فرازمینی واحد زمان را که به واحد متر مربع سطح می‌تابد با رابطه زیر بیان کرد:

$$\large \begin {align} G _ \mathrm { e x } = G _ \mathrm { o } \left ( 1 + 0 . 0 3 3 3 \cos \left ( \dfrac { 2 \pi N } { 3 6 5 } \right ) \right ) \end {align} \;\;\;\;\; ( 1 5 )$$

با‌ این‌ حال، مفهومی که اغلب در مدل‌های تابش خورشید استفاده می‌شود، این است که تابش فرازمینی به یک سطح افقی می‌تابد. یک سطح تخت را در خارج از جو زمین و به موازات سطح زمین در زیر آن در نظر بگیرید. هنگامی‌ که این سطح با خورشید مواجه می‌شود (با یک پرتوی مرکزی، نرمال است)، تابش $$G_\mathrm{ex}$$ روی آن است که حداکثر تابش خورشیدی در آن فاصله است. اگر سطح نسبت به خورشید نرمال نباشد، تابشی خورشیدی که روی آن فرود می‌آید، با کسینوس زاویه بین نرمال سطح و یک پرتوی مرکزی از خورشید کاهش می‌یابد. این مفهوم در شکل ۶ توضیح داده شده است. در این شکل مشاهده می‌شود که نرخ انرژی خورشیدی روی هر دو سطح یکسان است. با این‌ حال، مساحت سطح A بزرگ‌تر از تصویر آن، یعنی سطح فرضی B است که به همین دلیل، نرخ انرژی خورشیدی واحد سطح، روی سطح A کمتر از سطح B است.

بنابراین، تابش خورشیدی فرازمینی ($$G_\mathrm{exH}$$) روی سطح زمین در یک سطح افقی واقع در یک مکان خاص را می‌توان از رابطه (۱۶)‌ محاسبه کرد:

$$\large \begin {align} G _ \mathrm { e x H } = G _ \mathrm { e x } \cos \phi \end {align} \;\;\;\;\; ( 1 6 )$$

که در آن، $$\phi$$ زاویه اوج خورشید است که مستقیماً از بالا نسبت به مرکز هندسی قرص خورشید اندازه‌گیری می‌‌شود. اندازه زاویه اوج برابر با مقدار فراز است و بنابراین معادله (۱۶) را می‌توان به‌ صورت زیر بازنویسی کرد:

$$\large \begin {align} G _ \mathrm { e x H } = G _ \mathrm { o } \left ( 1 + 0 . 0 3 3 3 \cos \left ( \dfrac { 2 \pi N } { 3 6 5 } \right ) \right ) \sin L \sin \delta + \cos L \cos \delta \cos \omega \end {align} \;\;\;\;\; ( 17 )$$

در نهایت، کل انرژی خورشیدی فرازمینی ($$E_\mathrm{ex}$$) بر حسب $$\mathrm{Wh/m^2}$$ به‌ صورت زیر محاسبه می‌شود:

$$\large \begin {align} E _ \mathrm { e x } = \int \limits _ { T _ \mathrm { s r } } ^ { T _ \mathrm { s s } } { G _ \mathrm { e x H } d t } \end {align} \;\;\;\;\; ( 1 8 )$$

مثال ۴

برنامه‌ای بنویسید که پروفیل تابش فرازمینی ساعتی را برای شهر نابلس فلسطین در سی‌ و یکم مارس پیش‌بینی کند.

حل: با توجه به معادله (۱۷)، ابتدا باید مقادیر ساعتی زاویه فراز مکان انتخابی را محاسبه کرد (مثال ۲ را ببنید). سپس می‌توان مقدار تابش ساعتی را با استفاده از معادله (۱۷) به‌ صورت زیر به‌ دست آورد:

1%Example 1.4
2%---------------------------------------------------------
3%Date 31/3/2015 (N=90)
4%Location Nablus, Palestine, L =(32.22), LOD = (35.27)
5L=32.22; %Latitude
6LOD=35.27; %Longitude
7N=90; %day number
8T_GMT=+3; %Time difference with reference to GMT
9Step=60; %step each hour
10Ds=23.45*sin((360*(N-81)/365)*(pi/180)); % angle of
11declination
12B=(360*(N-81))/364; %Equation of time
13EoT=(9.87*sin(2*B*pi/180))- (7.53*cos(B*pi/180))- (1.5*sin
14(B*pi/180)); %Equation of time
15Lzt= 15* T_GMT; %LMST
16if LOD>=0
17Ts_correction= (-4*(Lzt-LOD))+EoT; %solar time correction
18else
19Ts_correction= (4*(Lzt-LOD))+EoT; %solar time correction
20end
21Wsr_ssi=- tan(Ds*pi/180)*tan(L*pi/180); %Sunrise/Sunset
22hour angle
23Wsrsr_ss=acosd(Wsr_ssi); %Sunrise/Sunset hour angle
24ASTsr=abs((((Wsrsr_ss/15)-12)*60)); %Sunrise solar time
25ASTss=(((Wsrsr_ss/15)+12)*60); %Sunset solar time Tsr=ASTsr+
26abs(Ts_correction) %Actual Sunrise time
27Tss=ASTss+abs(Ts_correction) %Actual Sunset time
28sin_Alpha=[];
29for LMT=Tsr:Step:Tss %for loop for the day time
30Ts= LMT + Ts_correction; %solar time at each step
31Hs=(15 *(Ts - (12*60)))/60; %hour angle at each step
32sin_Alpha_i=(sin(L*pi/180)*sin(Ds*pi/180))+
33(cos(L*pi/180)*cos(Ds*pi/180)* cos(Hs*pi/180)); %altitude
34angle
35sin_Alpha=[sin_Alpha;sin_Alpha_i]; % store altitude
36angle results
37end
38LMT=Tsr:Step:Tss
39sin_Alpha;
40Go=1367; %solar constant
41Gext=Go*(1+(0.0333*cos(360*N/365))); %available Gext
42GextH=Gext*sin_Alpha; %Gex on horizontal surface
43plot(LMT,GextH) %plot results

پاسخ: شکل ۷.

مدل‌سازی تابش زمینی روی یک سطح افقی

تابش کل ($$G_\mathrm{T}$$)، تابش موجود در سطح دریا و زیر جو زمین است. تابش کل که روی یک سطح افقی وجود دارد، شامل دو جزء مستقیم یا پرتو (Direct Solar Radiation) و پراکنده (Diffuse Solar Radiation) است. شکل ۸ اجزای تابش خورشیدی را در یک سطح افقی نشان می‌دهد. تابش مستقیم خورشید ($$G_\mathrm{B}$$)، پرتویی است که مستقیماً از خورشید ساطع می‌شود، در حالی‌ که تابش پراکنده خورشید ($$G_\mathrm{D}$$)، تابشی است که ابرها و ذرات دیگر در آسمان پراکنده می‌کنند.

بنابراین می‌توان $$G_\mathrm{T}$$ را به‌ شکل زیر توصیف کرد:

$$\large G _ \mathrm { T } = G _ \mathrm { B } + G _ \mathrm { D } \;\;\;\;\; ( 1 9 )$$

وقتی پرتوهای تابش خورشید از جو عبور می‌کنند، بسیاری از اجزای این پرتو‌ها توسط گازهای آسمان یا مولکول‌های هوا جذب، خنثی و پراکنده می‌شوند. برای یک روز با آسمان صاف، 70 درصد از تابش کلی، تابش مستقیم است. تضعیف این پرتو به ‌دلیل گرد و غبار، آلودگی هوا، بخار آب، ابرها و مه را می‌توان نسبتاً راحت مدل کرد. راه‌حل‌های زیادی برای مدل‌سازی این تضعیف به‌ عنوان تابعی از شماره روز وجود دارد. یکی از این مدل‌ها، مدل اَشری (ASHRAE Model) است که گاهی مدل آسمان صاف نامیده می‌شود. با این مدل، تابش مستقیم خورشید را که به سطح زمین می‌رسد ($$G_\mathrm{B,norm}$$) می‌توان به ‌صورت زیر نوشت:

$$\large G _ \mathrm { B , n o r m } = A ^ { \frac { - K } { \sin \alpha } } \; \; \; \; \; ( 2 0)$$

که در آن، $$A$$ شار فرازمینی ظاهری و $$K$$ یک عامل بدون بعد است که عمق نوری نامیده می‌شود. عامل‌های $$A$$ و $$K$$ را می‌توان به‌ عنوان توابعی از شماره روز بیان کرد:

$$\large A = 1 1 6 0 + 7 5 \sin \left [ \dfrac { 3 6 0} { 36 5 }( N - 2 75 ) \right ] \;\;\;\;\; ( 2 1)$$

$$\large K = 0 . 1 7 4 + 0 . 0 3 5 \sin \left [ \dfrac { 3 6 0 } { 3 65 } ( N - 1 0 0 ) \right ] \;\;\;\;\; ( 22 )$$

اکنون می‌توان تابش مستقیم جذب‌ شده $$G_\mathrm{B}$$ توسط سطح افقی را به‌ شکل زیر توصیف کرد:

$$\large G _ \mathrm { B } = G _ \mathrm { B , n o r m } \sin \alpha \;\;\;\;\; ( 23 )$$

از سوی دیگر، محاسبه تابش پراکنده روی یک کلکتور مسطح افقی، نسبت به محاسبه تابش مستقیم دشوارتر است. تابش ورودی، ممکن است با ذرات جو و بخار آب پراکنده و توسط ابرها منعکس شود. بعضی از تابش‌ها از سطح به آسمان بازتاب و دوباره به زمین پراکنده می‌شوند. در ساده‌ترین مدل‌های تابش پراکنده، فرض می‌شود که تابش با شدت برابر از تمام جهات به یک محل می‌رسد؛ به عبارت دیگر، آسمان به‌ صورت همسان در نظر گرفته می‌شود. بدیهی است که در روزهای مه‌آلود یا بارانی، آسمان در مجاورت خورشید بسیار درخشان‌تر از قبل است و اندازه‌گیری‌ها نیز پدیده مشابهی را در روزهای صاف نشان می‌دهد، اما این عوارض اغلب نادیده گرفته می‌شود. در نهایت، تابش پراکنده خورشید را می‌توان با رابطه (۲۴) تقریب زد:

$$\large G _ \mathrm { D} = 0. 0 9 5 + 0 . 0 4 \sin \left [ \dfrac { 3 6 0 } { 3 6 5 } ( N - 10 0 ) \right ] G _ \mathrm { B , n o r m } \;\;\;\;\; ( 24)$$

مثال ۵

برنامه‌ای در متلب بنویسید که پروفیل تابش پراکنده و کلی یک سطح افقی را در دوم مِی از طلوع تا غروب خورشید برای شهر کویت به‌ صورت ساعتی پیش‌بینی کند.

حل: برنامه مورد نظر را به دو بخش تقسیم می‌کنیم: محاسبه زاویه فراز و سپس محاسبه جزء تابش خورشید. بخش اول در مثال‌های قبلی مانند مثال ۲ نشان داده شد. برای بخش دوم، برنامه معادلات (۱۹) تا (۲۴)‌ به شرح زیر است:

1%Example 1.5
2%---------------------------------------------------------
3%Date 02/05/2015 (N=122)
4%Location Kuwait City, Kuwait, L =(29.36), LOD = (47.97)
5L=29.36; %latitude
6LOD=47.97; %longitude
7N=122; %day number
8T_GMT=+3; %time difference with reference to GMT
9Step=60; %time step
10Ds=23.45*(sind((360*(N-81)/365))); % angle of declination
11%=========================================================
12B=(360*(N-81))/364; %Equation of time
13EoT=(9.87*sin(2*B*pi/180))- (7.53*cos(B*pi/180))- (1.5*sin
14(B*pi/180)); %Equation of time
15Lzt= 15* T_GMT; %LMST
16if LOD>=0
17Ts_correction= (-4*(Lzt-LOD))+EoT; %solar time correction
18else
19Ts_correction= (4*(Lzt-LOD))+EoT; %solar time correction
20end
21Wsr_ssi=- tan(Ds*pi/180)*tan(L*pi/180); %Sunrise/Sunset
22hour angle
23Wsrsr_ss=acosd(Wsr_ssi); %Sunrise/Sunset hour angle
24ASTsr=abs((((Wsrsr_ss/15)-12)*60)); %Sunrise solar time
25ASTss=(((Wsrsr_ss/15)+12)*60); %Sunset solar time
26Tsr=ASTsr+abs(Ts_correction) %Actual Sunrise time
27Tss=ASTss+abs(Ts_correction) %Actual Sunset time
28sin_Alpha=[];
29for LMT=Tsr:Step:Tss %for loop for the day time
30Ts= LMT + Ts_correction; %solar time
31Hs=(15 *(Ts - (12*60)))/60; %Hour angle degree
32sin_Alpha_i=(sin(L*pi/180)*sin(Ds*pi/180))+
33(cos(L*pi/180)*cos(Ds*pi/180)* cos(Hs*pi/180)); %altitude
34angle
35sin_Alpha=[sin_Alpha;sin_Alpha_i]; % Store altitude
36angle
37end
38sin_Alpha
39%============solar radiation calculation===============
40A=1160+(75*sind((360/365)*(N-275))); %extraterrestrial
41solar energy flux
42k= 0.174+ (0.035*sind((360/365)*(N-100))); %k is a
43factor
44C= 0.095+ (0.04*sind((360/365)*(N-100))); %C is a
45factor
46G_B_norm=A*exp(-k./sin_Alpha) % available beam radiation
47in the sky
48G_B=G_B_norm.*sin_Alpha; %collected beam solar radiation
49by the collector on a horizontal surface
50G_D=C*G_B_norm; %diffuse on horizontal surface
51G_T= G_B+G_D
52%--actual global solar radiation data for Kuwait city
53*10e3----
54G_A=[000 0.2431 0.4422 0.5966 0.865 0.976 1.031 1.016
550.936 0.788 0.5904 0.3541 0.1439];
56LMT=Tsr:Step:Tss
57plot(LMT,G_T)
58holdon
59plot(LMT,G_A*1e3, 'red')

پاسخ: شکل ۹.

مدل‌سازی تابش کلی روی یک سطح صاف

اجزای تابش کلی یک کلکتور شیب‌دار در شکل ۱۰ نشان داده شده است.

فیلم آموزش برنامه نویسی متلب برای علوم و مهندسی در فرادرس

کلیک کنید

علاوه بر تابش مستقیم ($$G_{\mathrm{B}, \beta}$$) و تابش پراکنده ($$G_{\mathrm{D}, \beta}$$)، مولفه جدیدی به‌نام بازتاب تابش خورشیدی (Reflected Solar Radiation) $$G_\mathrm{R}$$ اضافه شده است که تابش سطح شیب‌دار را شکل می‌دهد.

این مؤلفه‌ها را می‌توان با رابطه (۲۵) نشان داد:

$$\large G _ \mathrm { T , \beta } = G _ \mathrm { B , \beta } + G _ \mathrm { D , \beta } + G _ \mathrm { R } \;\;\;\;\; ( 25 )$$

معادله (۲۵) را می‌توان با توجه به مؤلفه‌های انرژی خورشیدی در یک سطح افقی به‌ صورت زیر بازنویسی کرد:

$$\large G _ \mathrm { T , \beta } = G _ \mathrm { B } R _ \mathrm { B } + G _ \mathrm { D } R _ \mathrm { D } + G _ \mathrm { T } \rho R _ \mathrm { R } \;\;\;\;\; ( 2 6 )$$

که در آن، $$R_\mathrm{B}$$، $$R_\mathrm{D}$$ و $$R_\mathrm{R}$$ ضرایب هستند و $$\rho$$ «آلبِدو یا سپیدایی زمین» (Ground Albedo) به معنی درصد بازتاب نور از سطح است. همچنین، $$R_\mathrm{B}$$ نسبت بین انرژی خورشیدی کلی در یک سطح افقی و انرژی خورشیدی کلی روی سطح شیب‌دار، $$R_\mathrm{D}$$ نسبت بین انرژی خورشیدی پراکنده در یک سطح افقی و انرژی خورشیدی پراکنده سطح شیب‌دار و $$R_\mathrm{R}$$ عامل انرژی خورشیدی بازتاب‌ شده سطح شیب‌دار است.

از معادله (۲۶) واضح است که کلید یافتن مؤلفه‌های انرژی خورشیدی روی سطح شیب‌دار، تخمین ضرایب $$R_\mathrm{B}$$، $$R_\mathrm{D}$$ و $$R_\mathrm{R}$$ است. متداول‌ترین مدل برای محاسبه $$R_\mathrm{B}$$ «مدل لیو و جردن» (Liu and Jordan Model) است که در آن، $$R_\mathrm{B}$$ به‌ صورت زیر محاسبه می‌شود:

$$\large R _ \mathrm { B} = \dfrac { \cos ( L - \beta ) \cos \delta \sin \omega _ \mathrm { s s } + \omega _ \mathrm { s s } \sin ( L - \beta ) \sin \delta } { \cos L \cos \delta \sin \omega _ \mathrm { s s } + \omega _ \mathrm { s s } \sin L \sin \delta } \;\; \;\;\; ( 2 7 )$$

از آن‌جا که در سطوح در نیم‌کره جنوبی شیب به سمت استوا است، معادله $$R_\mathrm{B}$$ به‌ صورت زیر نشان داده می‌شود:

$$\large R _ \mathrm { B } = \dfrac { \cos ( L + \beta ) \cos \delta \sin \omega _ \mathrm { s s } + \omega _ \mathrm { s s } \sin ( L + \beta ) \sin \delta } { \cos L \cos \delta \sin \omega _ \mathrm { s s } + \omega _ \mathrm { s s } \sin L \sin \delta } \;\;\;\;\; ( 2 8 )$$$$\large R _ \mathrm { B } = \dfrac { \cos ( L + \beta ) \cos \delta \sin \omega _ \mathrm { s s } + \omega _ \mathrm { s s } \sin ( L + \beta ) \sin \delta } { \cos L \cos \delta \sin \omega _ \mathrm { s s } + \omega _ \mathrm { s s } \sin L \sin \delta } \;\;\;\;\; ( 2 8 )$$

در عین‌ حال، بهترین معادله برای $$R_\mathrm{R}$$، رابطه (۲۹) است:

$$\large R _ \mathrm { R } = \dfrac { 1 - \cos \beta } { 2 } \;\;\;\;\; ( 2 9 )$$

از سوی دیگر، برای $$R_\mathrm{D}$$ مدل‌های زیادی ارائه شده که می‌توان آن‌ها را در مدل‌های همسان و ناهمسان دسته‌بندی کرد.

مدل خورشیدی همسان، مبتنی بر این فرض است که تابش همسان بدون در نظر گرفتن جهت اندازه‌گیری، شدت یکسانی دارد و میدان همسان بدون در نظر گرفتن چگونگی ذرات آزمون، عمل می‌کند. این تابش در همه جهات به صورت یکنواخت از یک منبع نقطه‌ای می‌تابد که گاهی «تابنده یا رادیاتور همسان» (Isotropic Radiator) نامیده می‌شود. یکی از مدل‌های رایج تابش پراکنده، مدل لیو و جردن با $$R_\mathrm{D}$$ است:

$$\large R_\mathrm{D}=\dfrac{1+\cos \beta}{2} \;\;\;\;\; ( 30 )$$

$$\large R_\mathrm{D}=\dfrac{1}{3[2+\cos \beta]} \;\;\;\;\; ( 31 )$$

$$\large R_\mathrm{D}=\dfrac{3+\cos (2\beta)}{4} \;\;\;\;\; ( 3 2 )$$

$$\large R_\mathrm{D}=1-\dfrac{\beta}{180} \;\;\;\;\; ( 33 )$$

از سوی دیگر، ناهمسان بودن که به‌ معنای وابستگی به جهت است، در مقابل همسان بودن که به‌ معنای خواص یکسان در همه جهات است قرار می‌گیرد. این را می‌توان به‌ عنوان یک تفاوت اندازه‌گیری در امتداد محورهای مختلف ویژگی‌های فیزیکی مواد (جذب، شاخص شکست، چگالی و غیره) تعریف کرد. بنابراین، مدل‌های خورشیدی ناهمسان، مبتنی بر این فرض هستند که تابش ناهمسان بسته به جهت اندازه‌گیری، شدت متفاوتی دارد و در همه جهات، غیریکنواخت می‌تابد. چند مدل ناهمسان که برای $$R_\mathrm{D}$$ ارائه شده‌اند، به‌ صورت زیر هستند:

$$\large R _ \mathrm { D } = \dfrac { G _ \mathrm { B } } { G _ \mathrm { T } } R _ \mathrm { D } + \left ( 1 + \dfrac { G _ \mathrm { B } } { G _ \mathrm { T } } \right ) \left ( \dfrac { 1 + \cos \beta } { 2 } \right ) \;\;\;\;\; ( 34 )$$

$$\large \begin {align} R _ \mathrm { D } & = 0 . 5 1 R _ \mathrm { B } + \dfrac { 1 + \cos \mathrm { T L T } } { 2 } \nonumber \\ & \, \, \, \, \, \, \, - \dfrac { 1 . 7 4 } { 1 . 2 6 \pi } \left [ \sin \beta = \left ( \beta \dfrac { \pi } { 1 8 0 } \right ) \cos \beta -\pi \sin ^ 2 \left ( \dfrac { \beta } { 2 } \right ) \right ] \end {align} \;\;\;\;\; ( 35 )$$

$$\large R _ \mathrm { D } = \dfrac { G _ \mathrm { B } } { G _ \mathrm { T } } R _ \mathrm { B } + \left ( 1 -\dfrac { G _ \mathrm { B } } { G _ \mathrm { T } } \right ) \left ( \dfrac { 1 + \cos \beta } { 2 } \right ) \left ( 1 + \sqrt { \dfrac { G _ \mathrm { B } } { G _ \mathrm { T } } } \sin ^ 3 \left ( \dfrac { \beta }{ 2 } \right ) \right ) \;\;\;\;\; ( 3 6 )$$

مثال ۶

برنامه‌ای بنویسید که تابش پراکنده و کلی ساعتی را روی یک سطح شیب‌دار برای شهر کویت از طلوع تا غروب خوشید در 31 مارس پیش‌بینی کند.

حل: برنامه متلب این مثال به صورت زیر است:

1%Example 1.6
2%-----31/03/2015 (N=90)
3%Location Kuwait, Kuwait, L =(29.36), LOD = (47.97)
4L=29.36; %latitude
5LOD=47.97; %longitude
6N=90; %day number
7T_GMT=+3; %time difference with reference to GMT
8Step=60; %time step
9Ds=23.45*(sind((360*(N-81)/365))); % angle of
10declination
11%=========================================================
12B=(360*(N-81))/364; %Equation of time
13EoT=(9.87*sin(2*B*pi/180))- (7.53*cos(B*pi/180))- (1.5*sin
14(B*pi/180)); %Equation of time
15Lzt= 15* T_GMT; %LMST
16if LOD>=0
17Ts_correction= (-4*(Lzt-LOD))+EoT; %solar time correction
18else
19Ts_correction= (4*(Lzt-LOD))+EoT; %solar time correction
20end
21Wsr_ssi=- tan(Ds*pi/180)*tan(L*pi/180); %Sunrise/Sunset
22hour angle time
23Wsrsr_ss=acosd(Wsr_ssi) %Sunrise/Sunset hour angle time
24ASTsr=abs((((Wsrsr_ss/15)-12)*60)); %Sunrise solar time
25ASTss=(((Wsrsr_ss/15)+12)*60) %Sunset solar time
26Tsr=ASTsr+abs(Ts_correction); Sunrise time
27Tss=ASTss+abs(Ts_correction); Sunset time
28sin_Alpha=[];
29for LMT=Tsr:Step:Tss-60
30Ts= LMT + Ts_correction; % solar time
31Hs=(15 *(Ts - (12*60)))/60; % Hour angle degree
32sin_Alpha_i=(sin(L*pi/180)*sin(Ds*pi/180))+
33(cos(L*pi/180)*cos(Ds*pi/180)* cos(Hs*pi/180)); %altitude
34angle
35sin_Alpha=[sin_Alpha;sin_Alpha_i];
36end
37sin_Alpha;
38%=======================================================
39A=1160+(75*sind((360/365)*(N-275))); %extraterrestrial
40solar energy flux
41k= 0.174+ (0.035*sind((360/365)*(N-100))); %k is a factor
42C= 0.095+ (0.04*sind((360/365)*(N-100))); %C is a factor
43%-----calculation of solar radiation on horizontal surface
44----------------
45G_B_norm=A*exp(-k./sin_Alpha) % available beam radiation
46in the sky
47G_B=G_B_norm.*sin_Alpha; %collected beam solar radiation
48by the collector on a horizontal surface
49G_D=C*G_B_norm; %diffuse on horizontal surface
50G_T= G_B+G_D;
51%-----calculation of solar radiation on tilt surface
52----------------
53Beta=L; %tilt angle
54Rb=((cos((L-Beta).*(pi/180)).*cos(Ds.*(pi/180)).*sin(Ws
55rsr_ss.*(pi/180))+ ...
56... (Wsrsr_ss.*(pi/180)).*sin((L-Beta).*(pi/180)).*sin(Ds.
57*(pi/180))))./...
58... (((cos(L*(pi/180)).*cos(Ds.*(pi/180)).*sin(Wsrsr_ss.
59*(pi/180)))+...
60... ((Wsrsr_ss.*(pi/180)).*sin(L*(pi/180)).*sin(Ds.*(pi/
61180)))));
62Rd=(1+cos(Beta*(pi/180)))./2;
63Rr= (0.3*(1-cos(Beta*(pi/180))))./2;
64G_B_Beta=(G_B.*Rb);
65G_D_Beta=(G_D.*Rd);
66G_R=(G_T.*Rr);
67G_T_Beta=G_B_Beta+G_D_Beta+G_R;
68LMT=Tsr:Step:Tss-60;
69plot(LMT,G_T)
70holdon
71plot(LMT,G_T_Beta, 'k')

پاسخ: شکل ۱۱

مدل‌سازی تابش براساس اندازه‌گیر‌ی‌های زمین

متغیرهای ورودی که معمولاً برای پیش‌بینی انرژی کلی در نقاط مختلف مورد استفاده قرار می‌گیرند نرخ نور (تابش) خورشید، دمای محیط و رطوبت نسبی هستند.

فیلم آموزش برنامه نویسی متلب MATLAB در فرادرس

کلیک کنید

انرژی کلی در یک سطح افقی، میانگین تابش کلی ضرب در زمان روز خورشیدی است. زمان روز خورشیدی ($$S_\mathrm{o}$$) با رابطه (۳۷) محاسبه می‌شود:

$$\large S _ \mathrm { o } = \dfrac { 2 } { 1 5 } \cos ^ { - 1 } ( - \tan L \tan \delta ) \;\;\;\;\; ( 3 7 )$$

انرژی کلی شدیداً به عاملی به‌ نام نرخ تابش خورشید وابسته است. نرخ تابش بالا، به ‌معنی انرژی زیاد است و برعکس. نرخ خورشید با رابطه (۳۸) داده می‌شود:

$$\large \dfrac { S } { S _ \mathrm { o } } = \dfrac { ( 1 0 - 1 . 2 5 C ) } { 1 0 } , \, \, \, 0 < \dfrac { S } { S _ \mathrm { o } } < 1 \;\;\;\;\; ( 3 8 )$$

که در آن، $$S_\mathrm{c}$$ تعداد ساعات تابان در یک روز آفتابی است و $$C$$ ابعاد میانگین روزانه اندازه‌گیری‌ شده با مقادیر $$0$$ برای آسمان صاف تا $$10$$ برای آسمان کاملاً ابری است.

مدل‌سازی تابش کلی

مدل‌های زیادی برای مدل‌سازی انرژی کلی خورشید روی سطح افقی ارائه شده است. به‌ طور کلی، چهار نوع مدل خطی، غیرخطی، فازی و مبتنی بر هوش مصنوعی برای مدل‌سازی انرژی کلی خورشید وجود دارد. تابش کلی با ساعت خورشیدی رابطه خطی دارد. بنابراین، برای محاسبه انرژی کلی، می‌توان یک مدل خطی براساس ساعت خورشید طراحی کرد. یکی از مدل‌های متداول در معادله (۳۹) آمده است:

$$\large \dfrac { G _ \mathrm { T } }{ G _ \mathrm { e x } } = a + b \dfrac { S } { S _ \mathrm { o } } \;\;\;\;\; ( 3 9 )$$

راه‌های مختلفی برای تعیین ضرایب $$a$$ و $$b$$ به صورت بهینه وجود دارد. اولین مدل تخمین انرژی کلی که به‌ عنوان «مدل آنگستروم» (Angström Model) شناخته می‌شود، از داده‌های مدت تابش خورشید حاصل می‌شود و مقادیر $$a$$ و $$b$$ را به صورت زیر تعریف می‌کند:

$$\large b = \dfrac { \sum _ { i - 1 } ^ { n } { \left [ \left ( \dfrac { G _ \mathrm { T } } { G _ \mathrm { e x } } \right ) _ i -\overline { \left ( \dfrac { S } { S _ \mathrm { o } } \right ) } _ i \right ] \left [ \left ( \dfrac { G _ \mathrm { T } } { G _ \mathrm { e x } } \right ) _ { i - 1 } - \overline { \left ( \dfrac { S } { S _ \mathrm { o } } \right ) } _ { i - 1 } \right ] } } { \sqrt { \sum _ { i - 1 } ^ { n } { \left [ \left ( \dfrac { G _ \mathrm { T } } { G _ \mathrm { e x } } \right ) _ i -\overline { \left ( \dfrac { S } { S _ \mathrm { o } } \right ) } _ { i - 1 } \right ] ^ 2 \left [ \left ( \dfrac { G _ \mathrm { T } } { G _ \mathrm { e x } } \right ) _ { i - 1 } -\overline { \left ( \dfrac { S } { S _ \mathrm { o } } \right ) } _ { i - 1 } \right ] ^ 2 } } } \;\;\;\;\; ( 40)$$

$$\large a = \overline { \left ( \dfrac { G _ \mathrm { T } } { G _ \mathrm { e x } } \right ) } - b \overline { \left ( \dfrac { S } { S _ \mathrm { o } } \right ) } \;\;\;\;\; ( 41 )$$

مثال ۷

با استفاده از نرم‌افزار متلب، یک مدل خطی برای میانگین ماهیانه تابش خورشیدی براساس اطلاعات جدول زیر بنویسید. فرض کنید زمان خورشیدی برابر با زمان محلی است.

حل: برنامه مربوط به این مثال به صورت زیر است:

1%Modeling of PV systems using MATLAB
2%Example 1.7
3fileName = 'PV Modeling Book Data Source.xls';
4sheetName ='Source 1' ;
5L=3.11;
6LOD=101.6;
7Go=1367;
8%----------------------------------------------------
9G_T=xlsread(fileName, sheetName , 'E5:E3640');
10S_So=xlsread(fileName, sheetName , 'I5:I3640');
11N=xlsread(fileName, sheetName , 'C5:C3640');
12LMT= xlsread(fileName, sheetName , 'D5:D3640');
13Ts=LMT; %assumption
14%--------------calculation of Gext-----------------
15Ds=23.45*sin((360*(N-81)/365)*(pi/180));% angle of
16declination
17Hs=15 *(Ts - 12); % Hour angle degree
18sin_Alpha=(sin(L*pi/180).*sin(Ds.*pi/180))+ (cos(L*pi/18
190)*cos(Ds.*pi/180).* cos(Hs.*pi/180)); %altitude angle
20Gext=Go*(1+(0.0333*cos(360*N/365)));
21GextH=Gext.*sin_Alpha;
22G_T_G_ext=G_T./GextH;
23%---------------Modeling of global solar energy----------
24---------------
25N_Liner=1;%order of the model ex nonlinear model (Eq.28),
26N_Liner=2
27P_Liner=polyfit(S_So,G_T_G_ext,N_Liner)
28X_Liner=S_So;
29Y_Liner=0;
30fori=1:N_Liner+1 %constant calculation
31Y_Liner=Y_Liner+P_Liner(i)*X_Liner.^(N_Liner-i+1);
32end
33plot(S_So,G_T_G_ext)
34holdon
35plot(X_Liner,Y_Liner,'-red')

پاسخ:

$$\large \dfrac { G _ \mathrm { T } }{ G _ \mathrm { e x } } = 0 . 2 7 7 2 + 0 . 2 7 7 2 \dfrac { S } { S _ \mathrm { o } } \;\;\;\;\; ( 42 )$$

از سوی دیگر، می‌توان مدلی را با اضافه کردن یک جمله غیرخطی به مدل آنگستروم به صورت زیر پیشنهاد کرد:

$$\large \dfrac { G _ \mathrm { T } } { G _ \mathrm { e x } } = a + b \dfrac { S } { S _ \mathrm { o } } + c \left ( \dfrac { S } { S _ \mathrm { o } } \right ) ^ 2 \;\;\;\;\; ( 4 3 )$$

تعدادی از مطالعات نظری نشان داده است که حساسیت خواص پرتوافکنی ابرها نسبت به ساختار فضایی آن‌ها، بر مدت تابش خورشید و در نتیجه، مقدار انرژی کلی اثر می‌گذارد. بنابراین، در هر فرایند تخمینی، ضروری است که لحظات آماری مرتبه‌ سوم یا بالاتر در نظر گرفته شود. یکی از راه‌های ادغام این لحظات در مدل، گنجاندن توان سوم متغیر مدت خورشید است:

$$\large \dfrac { G _ \mathrm { T } } { G _ \mathrm { e x } } = a + b \dfrac { S } { S _ \mathrm { o } } + c \left ( \dfrac { S } { S _ \mathrm { o } } \right ) ^ 2 + d \left ( \dfrac { S } { S _ \mathrm { o } } \right ) ^ 3 \;\;\;\;\; ( 44 )$$

معمولاً در انرژی و زمان تابش خورشید در طول روز ابهاماتی وجود دارد. برای مقابله با این عدم قطعیت و تخمین مقدار تابش می‌توان از الگوریتم «منطق فازی» (Fuzzy logic) استفاده کرد. مزیت اصلی مدل فازی این است که توانایی توصیف دانش را به‌ شیوه توصیفی انسان و در قالب قواعد ساده، تنها با استفاده از متغیرهای زبانی دارد. گزاره‌های منطقی فازی به‌ صورت عبارت‌های اگر-آنگاه هستند، به‌ عنوان مثال:

• اگر مدت تابش نور آفتاب «طولانی» باشد، آنگاه مقدار انرژی خورشیدی «زیاد» است.
• اگر مدت تابش نور آفتاب «کوتاه» باشد، آنگاه مقدار انرژی خورشیدی «کم» است.

در این دو گزاره، متغیرهای مدت تابش و انرژی، با متغیرهای زبانی مانند «طولانی»، «زیاد»، «کوتاه» و «کم» توصیف می‌شوند. در واقع، این دو گزاره منطقی در یک مدل ساده آنگستروم صدق می‌کنند.

مدل‌سازی تابش پراکنده

بسیاری از مدل‌های خطی رابطه بین $$G_\mathrm{D}/G_\mathrm{T}$$ و شاخص شفافیت $$K_\mathrm{T}$$ را توصیف کرده‌اند که برابر با $$G_\mathrm{T}/G_\mathrm{ex}$$ است. معادله عمومی مدل خطی‌ که انرژی خورشیدی پراکنده را محاسبه می‌کند می‌توان به‌ صورت زیر بیان کرد:

$$\large \dfrac { G _ \mathrm { D } }{ G _ \mathrm { T } } = a + b K _ \mathrm { T } \;\;\;\;\; ( 45 )$$

که در آن، $$a$$ و $$b$$ ضرایب مدل هستند.

می‌توان رابطه مشابهی بین $$G_\mathrm{D}/G_\mathrm{T}$$ و شاخص شفافیت $$K_\mathrm{T}$$ با مدل غیرخطی (۴۶) توصیف کرد:

$$\large \dfrac { G _ \mathrm { D } }{ G _ \mathrm { T } } = a + b K _ \mathrm { T } + c { K _ \mathrm { T } } ^ 2 + d { K _ \mathrm { T } } ^ 3 \;\;\;\;\; ( 4 6 )$$

مثال ۸

مثال ۷ را برای تابش پراکنده تکرار کنید.

حل: برنامه مربوط به این مثال به صورت زیر است:

1%Modeling of PV systems using MATLAB
2%Example 1.7
3fileName1 = 'PV Modeling Book Data Source.xls';
4sheetName = 'Source 1' ;
5L=3.11;
6LOD=101.6;
7Go=1367;
8%----------------------------------------------------
9G_T=xlsread(fileName1, sheetName , 'E5:E3640');
10G_D=xlsread(fileName1, sheetName , 'F5:F3640');
11N=xlsread(fileName1, sheetName , 'C5:C3640');
12LMT= xlsread(fileName1, sheetName , 'D5:D3640');
13Ts=LMT; %assumption
14%--------------calculation of Gext-----------------
15Ds=23.45*sin((360*(N-81)/365)*(pi/180));% angle of
16declination
17Hs=15 *(Ts - 12); Hour angle degree
18sin_Alpha=(sin(L*pi/180).*sin(Ds.*pi/180))+ (cos(L*pi/180)
19*cos(Ds.*pi/180).* cos(Hs.*pi/180)); %altitude angle
20Gext=Go*(1+(0.0333*cos(360*N/365)));
21GextH=Gext.*sin_Alpha;
22G_T_G_ext=G_T./GextH;
23G_D_G_T=G_D./G_T;
24%---------------Modeling of diffuse solar energy------
25-----
26N_Liner=1;
27P_Liner=polyfit(G_T_G_ext,G_D_G_T,N_Liner)
28X_Liner=G_T_G_ext;
29Y_Liner=0;
30fori=1:N_Liner+1
31Y_Liner=Y_Liner+P_Liner(i)*X_Liner.^(N_Liner-i+1);
32end
33plot(G_T_G_ext,G_D_G_T)
34hold on
35plot(X_Liner,Y_Liner,'-red')

پاسخ: شکل ۱۳

روش‌های هوش مصنوعی برای مدل‌سازی تابش خورشیدی

شبکه‌های عصبی مصنوعی (Artificial Neural Networks) یا ANN، سیستم‌های پردازش اطلاعاتی غیرالگوریتمی و موازی هستند که رابطه بین ورودی‌ها و خروجی‌ها را با نمونه‌هایی مانند داده‌های ثبت‌ شده قبلی یاد می‌گیرند. در این شبکه‌ها، نرون‌ها از طریق تعداد زیادی از پیوندهای وزن‌دار که می‌توانند سیگنال‌ها را عبور دهند به هم وصل شده‌اند. یک نرون، ورودی‌ها را با اتصالات ورودی خود دریافت می‌کند، آن‌ها را ترکیب می‌کند، در حالت کلی یک عملیات غیرخطی انجام می‌دهد و نتایج نهایی را به خروجی می‌برد. در اینجا از برنامه متلب برای آموزش و تشکیل ANN برای مدل‌سازی انرژی کلی خورشید استفاده می‌کنیم. همچنین از یک «شبکه پیش‌خور پرسپترون چندلایه» (Feedforward Multilayer Perceptron Network) یا FFMLP استفاده شده است، زیرا رایج‌ترین شبکه عصبی است که از نمونه‌ها یاد می‌گیرد و برای این کار مناسب است. نمودار ساختار پایه‌ای شبکه عصبی، در شکل ۱۴ نشان داده شده است. این شبکه سه لایه دارد: ورودی، مخفی و خروجی. هر لایه از طریق اتصالاتی با اندازه مختلف به نام وزن ارتباط دارد.

از چهار متغیر جغرافیایی و اقلیمی به‌ عنوان ویژگی‌های ورودی FFMLP استفاده می‌کنیم. این متغیرها عبارتند از شماره روز، طول جغرافیایی، عرض جغرافیایی و نسبت ساعت تابش روزانه (مدت تابش اندازه‌گیری‌ شده در طول روز بر حداکثر مدت تابش ممکن روز). یک گره خروجی به‌ عنوان پیش‌بینی شاخص شفافیت روزانه به‌ عنوان خروجی در نظر گرفته شده است. تابع تبدیل نرون‌ها یک تابع سیگموئید لگاریتمی $$f(z_i)$$ است:

$$\large f ( z _ i ) = \dfrac { 1 } { 1 + e ^ { - z _ i } } \;\;\;\;\; ( 47 )$$

$$\large z _ i = \sum \limits _ { j = 1 } ^ { 4 } { w _ { i j } x _ j + \beta _ i } \;\;\;\;\; ( 48 )$$

که در آن، $$z_i$$ مجموع وزن‌دار ورودی‌ها، $$x_j$$ سیگنال وارد شده به نرون $$j$$اُم از لایه ورودی، $$w_{ij}$$ وزن پیوند نرون $$j$$ به نرون $$i$$ در لایه مخفی و $$\beta _i$$ بایاس نرون $$i$$اُم است.

شبکه‌های عصبی برای حل مسئله آموزش داده می‌شوند، نه اینکه برای انجام کاری برنامه‌ریزی شوند. یادگیری از طریق آموزش حاصل می‌شود. به‌ عبارت دیگر، آموزش، فرایندی است که در آن، شبکه‌ها یاد می‌گیرند خطای بین مقدار خروجی و مقدار مرجع را به حداقل برسانند. در FFMLP، خطا با استفاده از الگوریتم آموزش اصطلاحاً پس‌انتشار، از لایه خروجی به لایه مخفی و سپس لایه ورودی به عقب منتشر می‌شود و با پس‌انتشار خطا، وزن پیوندهای بین نرون‌ها به‌روز می‌شود. در صورتی‌ که تعداد نرون‌های مخفی کافی باشد، با یک شبکه چندلایه می‌توان هر تابع چندمتغیره پیوسته را با هر درجه دقت به‌ صورت ریاضی تقریب زد. یک مشکل احتمالی که ممکن است رخ دهد، این است که شبکه به‌ جای یادگیری و تعمیم ساختار اصلی داده‌ها، جزئیات غیرضروری موارد خاص آن‌ها را یاد بگیرد.

همان توپولوژی شبکه عصبی که در مدل‌سازی انرژی کلی مورد استفاده قرار گرفت، در مدل‌سازی انرژی خورشیدی پراکنده نیز استفاده می‌شود. شکل ۱۵، شبکه عصبی پیش‌بین انرژی پراکنده نشان می‌دهد. این شبکه، چهار ورودی شامل عرض جغرافیایی، طول جغرافیایی، شاخص شفافیت و شماره روز و یک متغیر خروجی دارد که انرژی پراکنده خورشیدی است. شکل ۱۶، این شبکه عصبی را نشان می‌دهد.

مثال ۹

یک شبکه عصبی FFMLP ارائه دهید که تابش کلی ساعتی و تابش پراکنده را همان‌گونه که در شکل ۱۶ نشان داده شده پیش‌بینی کند.

حل: برنامه مربوط به این مثال به صورت زیر است:

1%Modeling of PV systems using MATLAB
2%Example 1.9
3fileName = 'PV Modeling Book Data Source.xls';
4sheetName = 'Source 1' ;
5G_T=xlsread(fileName, sheetName , 'E5:E2997');
6G_D=xlsread(fileName, sheetName , 'F5:F2997');
7Hum=xlsread(fileName, sheetName , 'H5:H2997');
8T=xlsread(fileName, sheetName , 'J5:J2997');
9S=xlsread(fileName, sheetName , 'I5:I2997');
10M=xlsread(fileName, sheetName , 'A5:A2997');
11D=xlsread(fileName, sheetName , 'B5:B2997');
12H=xlsread(fileName, sheetName , 'D5:D2997');
13G_T_Test=xlsread(fileName, sheetName , 'E2998:E3640');
14G_D_Test=xlsread(fileName, sheetName , 'F2998:F3640');
15Hum_Test=xlsread(fileName, sheetName , 'H2998:H3640');
16T_Test=xlsread(fileName, sheetName , 'J2998:J3640');
17S_Test=xlsread(fileName, sheetName , 'I2998:I3640');
18M_Test=xlsread(fileName, sheetName , 'A2998:A3640');
19D_Test=xlsread(fileName, sheetName , 'B2998:B3640');
20H_Test=xlsread(fileName, sheetName , 'D2998:D3640');
21%------------
22inputs = [M,D,H,T,H,S];
23I=inputs';
24targets= [G_T, G_D];
25T=targets';
26%-------ann Model development and training
27net = newff(I,T,5);
28Y = sim(net,I);
29net.trainParam.epochs = 100;
30net = train(net,I,T);
31%------------testing the developed model--------
32test=[M_Test,D_Test,H_Test,T_Test,H_Test,S_Test];
33Test1=test';
34G_Mi = sim(net,Test1);
35G_M= G_Mi';
36G_Tp=[];
37G_Dp=[];
38for i=1:1:length(G_M)
39G_Tp=[G_Tp;G_M(i,1)];
40G_Dp=[G_Dp;G_M(i,2)];
41end
42G_Tp;
43G_Dp;
44subplot(2,1,1)
45plot (G_T_Test)
46hold on
47plot(G_Tp,'red')
48subplot(2,1,2)
49plot (G_D_Test)
50hold on
51plot(G_Dp,'red')