-
Notifications
You must be signed in to change notification settings - Fork 0
/
OF.m
197 lines (178 loc) · 8.21 KB
/
OF.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
function [HST_max,TOT_max,AEQ,HST,TOT]=OF(PUL,AMB)
%% This function calculates thermal mode for OF power transformer per IEC 60076-7
% Input data:
% PUL - profile of transformer loading in pu
% AMB - profile of ambient temperature, degC
% Output data:
% HST_max - maximal hot-spot temerature of winding, degC
% TOT_max - maximal top-oil temperature, degC
% HST - profile of hot spot temperature, degC
% TOT - profile of top oil temperature, degC
% AEQ - ageing equivalent, pu relatve to normal ageing (1pu)
%% Constants
% Thermal characteristics of OF power transformer per IEC60076-7
delta_theta_or = 56; % Top-oil (in tank) temperature rise in steady state at rated losses (no-load losses + load losses),K
delta_theta_hr = 22; % Hot-spot-to-top-oil (in tank) gradient at rated current, K
tao_0 = 90; % Average oil time constant, min
tao_w = 7; % Winding time constant, min
R = 6; % Ratio of load losses at rated current to no-load losses
x = 1; % Exponential power of total losses versus top-oil (in tank) temperature rise (oil exponent)
y = 1.3; % Exponential power of current versus winding temperature rise (winding exponent)
k11 = 1; % Thermal model constant
k21 = 1.3; % Thermal model constant
k22 = 1; % Thermal model constant
%% Solving the difference equations
NPUL=length(PUL); % Finding the length of load data
NAMB=length(AMB); % Finding the length of ambient temperature data
% Checking that input data is in minute format
if NPUL==1440 && NAMB==1440
% do nothing
elseif NPUL==24 && NAMB==24
PUL=Convert2minute(PUL); % Convert loading data to minute format
AMB=Convert2minute(AMB); % Convert amb. temperature data to minute format
elseif NPUL==1440 && NAMB==24
AMB=Convert2minute(AMB); % Convert amb. temperature data to minute format
elseif NPUL==24 && NAMB==1440
PUL=Convert2minute(PUL); % Convert loading data to minute format
else
error('Check the length of input data')
end
% Change the variables name
K=PUL;
theta_a=AMB;
Dt=1;% time step 1 minute
% Although the system may not strictly be in the steady state at the start of a calculation period,
% this is usually the best one can assume, and it has little effect on the result
K_0=K(1);
theta_a_0=theta_a(1);
theta_0 = ((1+K_0.^2.*R)./(1+R)).^x.*delta_theta_or+theta_a_0; % top-oil temperature
delta_theta_h1 = k21*K_0.^y*delta_theta_hr; % Hot-spot-to-top-oil (in tank) gradient at start
delta_theta_h2 = (k21-1)*K_0.^y*delta_theta_hr; % Hot-spot-to-top-oil (in tank) gradient at start
Loss_of_life = 0;
% Create an array of hot-spot temperature and top-oil temperature
HST=NaN(length(K),1);
TOT=NaN(length(K),1);
% Solving difference (not differentiate!) equations: iterative approach (see
% Annex C in IEC 60076-7 for equations)
for i=1:1:length(K)
D_theta_0 = (Dt./(k11.*tao_0)).*((((1+K(i).^2.*R)./(1+R)).^x).*(delta_theta_or)-(theta_0-theta_a(i)));
theta_0 = theta_0+D_theta_0;
D_delta_theta_h1 = Dt./(k22.*tao_w).*(k21.*delta_theta_hr.*K(i).^y-delta_theta_h1);
delta_theta_h1 = delta_theta_h1+D_delta_theta_h1;
D_delta_theta_h2 = Dt./(1./k22.*tao_0).*((k21-1).*delta_theta_hr.*K(i).^y-delta_theta_h2);
delta_theta_h2 = delta_theta_h2+D_delta_theta_h2;
delta_theta_h = delta_theta_h1-delta_theta_h2;
HST(i,:) = theta_0+delta_theta_h; % hot spot temperature
TOT(i,:)=theta_0; % top oil temperature
end
% Calculating ageing
DL=NaN(length(K),1);
for i=1:1:length(HST)
% DL(i,:) = (exp((15000./(110+273)-15000./(HST(i)+273)))).*Dt;
DL(i,:) = (2^((HST(i)-98)/6)).*Dt;
end
Loss_of_life = Loss_of_life+DL;
ASUM=sum(Loss_of_life);
% Last outputs
AEQ=ASUM/length(K);
HST_max=max(HST);
TOT_max=max(TOT);
end% end of function
% % Äàííûé ñêðèïò îïðåäåëÿåò èçíîñ èçîëÿöèè è òåìïåðàòóðó ÍÍÒ
% % Òåïëîâàÿ ìîäåëü ñèëîâîãî òðàíñôîðìàòîðà OF IEC 60076-7
% % Ââîä ïàðàìåòðîâ òðàíñôîðìàòîðà èç òàáëèöû â ïðèëîæåíèè Å
% delta_theta_or = 56; % input('ââåäèòå delta_theta_or '); % Ðîñò òåìïåðàòóðû â âåðõíèõ ñëîÿõ â íîðìàëüíîì ðåæèìå, Ê
% delta_theta_hr = 22; % input ('ââåäèòå delta_theta_hr '); % Ïðåâûøåíèå òåìïåðàòóðû ÍÍÒ íàä òåìïåðàòóðîé â âåðõíèõ ñëîÿõ ìàñëà, Ê
% tao_0 = 90; % input ('ââåäèòå tao_0 '); % Ïîñòîÿííàÿ âðåìåíè ìàñëà, ìèí
% tao_w = 7; % input ('ââåäèòå tao_w '); % Ïîñòîÿííàÿ âðåìåíè îáìîòîê, ìèí
% R = 6; % input ('ââåäèòå R '); % ïîñòîÿííàÿ îòíîøåíèå ïîòåðü ïðè íîìèíàëüíîé íàãðóçêå ê ïîòåðÿì íà õîëîñòîì õîäó, ìèí
% x = 1; % input ('ââåäèòå x ');
% y = 1.3; % input ('ââåäèòå y ');
% k11 = 1; % input ('ââåäèòå k11 ');
% k21 = 1.3; % input ('ââåäèòå k21 ');
% k22 = 1; % input ('ââåäèòå k22 ');
% %t=5;
%
% %Ââîä èñõîäíûõ äàííûõ: ãðàôèê íàãðóçêè, òåìïåðàòóðû
% % AMB=linspace(33,33,1440)';
% NPUL=length(PUL); % Finding the length of load data
% NAMB=length(AMB); % Finding the length of ambient temperature data
%
% % Checking that input data is in minute format
% if NPUL==1440 && NAMB==1440
% % do nothing
% elseif NPUL==24 && NAMB==24
% PUL=PUL_to_1min(PUL,60); % Convert loading data to minute format
% AMB=PUL_to_1min(AMB,60); % Convert amb. temperature data to minute format
% elseif NPUL==1440 && NAMB==24
% AMB=PUL_to_1min(AMB,60); % Convert amb. temperature data to minute format
% elseif NPUL==24 && NAMB==1440
% PUL=PUL_to_1min(PUL,60); % Convert loading data to minute format
% elseif NPUL==48 && NAMB==48
% PUL=PUL_to_1min(PUL,30); % Convert loading data to minute format
% AMB=PUL_to_1min(AMB,30); % Convert amb. temperature data to minute format
% elseif NPUL==96 && NAMB==96
% PUL=PUL_to_1min(PUL,15); % Convert loading data to minute format
% AMB=PUL_to_1min(AMB,15); % Convert amb. temperature data to minute format
% elseif NPUL==288 && NAMB==288
% PUL=PUL_to_1min(PUL,5); % Convert loading data to minute format
% AMB=PUL_to_1min(AMB,5); % Convert amb. temperature data to minute formatelse
% else
% error('Check the length of input data')
% end
% K=PUL;
% theta_a=AMB;
% % load('initial_data.mat','TIM')
% Dt=1;
%
%
% % Ðàñ÷åò íà÷àëüíûõ óñëîâèé
% K_0=K(1);
% theta_a_0=theta_a(1);
% theta_0 = ((1+K_0.^2.*R)./(1+R)).^x.*delta_theta_or+theta_a_0; % íà÷àëüíîå çíà÷åíèå òåìïåðàòóðû â âåðõíèõ ñëîÿõ ìàñëà â áàêå
% delta_theta_h1 = k21*K_0.^y*delta_theta_hr; % íà÷àëüíîå çíà÷åíèå ïðåâûøåíèÿ òåìïåðàòóðû ÍÍÒ íàä òåìïåðàòóðîé âåðõíèõ ñëîåâ ìàñëà â áàêå
% delta_theta_h2 = (k21-1)*K_0.^y*delta_theta_hr; % íà÷àëüíîå çíà÷åíèå ïðåâûøåíèÿ òåìïåðàòóðû ÍÍÒ íàä òåìïåðàòóðîé âåðõíèõ ñëîåâ ìàñëà â áàêå
% Loss_of_life = 0;
%
% % Create an array of hot-spot temperature and top-oil temperature
% HST=NaN(length(K),1);
% TOT=NaN(length(K),1);
%
% % Ðåøåíèå ðàçíîñòíûõ óðàâíåíèé
% for i=1:1:length(K)
%
% D_theta_0 = (Dt./(k11.*tao_0)).*((((1+K(i).^2.*R)./(1+R)).^x).*(delta_theta_or)-(theta_0-theta_a(i)));
% theta_0 = theta_0+D_theta_0;
%
% D_delta_theta_h1 = Dt./(k22.*tao_w).*(k21.*delta_theta_hr.*K(i).^y-delta_theta_h1);
% delta_theta_h1 = delta_theta_h1+D_delta_theta_h1;
%
% D_delta_theta_h2 = Dt./(1./k22.*tao_0).*((k21-1).*delta_theta_hr.*K(i).^y-delta_theta_h2);
% delta_theta_h2 = delta_theta_h2+D_delta_theta_h2;
%
% delta_theta_h = delta_theta_h1-delta_theta_h2;
%
% HST(i,:) = theta_0+delta_theta_h; % Òåìïåðàòóðà ÍÍÒ â òðàíñôîðìàòîðå
% TOT(i,:)=theta_0;
% end
% % Calculating ageing
% DL=NaN(length(K),1);
% % Ïîòåðÿ ñðîêà ñëóæáû
% for i=1:1:length(HST)
% % DL(i,:) = (exp((15000./(110+273)-15000./(HST(i)+273)))).*Dt;
% DL(i,:) = (2^((HST(i)-98)/6)).*Dt;
% end
% Loss_of_life = Loss_of_life+DL;
% ASUM=sum(Loss_of_life);
% % Current_ageing=0;
% % for i=1:length(DL)
% % Current_ageing(i)=Current_ageing(end)+DL(i);
% % end
% AEQ=ASUM/length(K);
% HST_max=max(HST);
% TOT_max=max(TOT);
% HST_1=HST(1);
% HST_end=HST(end);
% end
%
%