-
Notifications
You must be signed in to change notification settings - Fork 1
/
cable_thermal_model_IEC_60853_2.m
1011 lines (851 loc) · 44.5 KB
/
cable_thermal_model_IEC_60853_2.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
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
function [T_max,Temperature,I_adm]=cable_thermal_model_IEC_60853_2(Iech_load)
%% ************************************************************************
%--------------------------------------------------------------------------
% Input of cable_thermal_model_IEC_60853_2.m:
% Iech_load - hourly current profile in A;
% Output of cable_thermal_model_IEC_60853_2.m:
% Temperature - hourly profile of cable temperature in degC from t=00:00 to
% t_end=length(Iech_load)e.g.[25x1] 25 since 1 value for t=0;
% T_max - maximal temperature of cable conductor, [1x1];
% I_adm - steady-state current, calculated in accordance with IEC287
% [1x1];
%% ***********************************************************************
%% ************************************************************************
%% Electrothermal modelling of submarine cable
%% ************************************************************************
% Global variables (used in different functions)
global S_ame_pertes Lcable Rho_0_E_Al Alpha_E_Al Tinit Lambda1 Lambda2 Alpha_E_Pb
global T1 T2 T3 T4 TA TB Ta Tb Qc Qi Qs Qf Qa Qj QA QB Wd t a b I_D_1
global pVW Depth De Tamb Rho_sol Diff_sol R_DC20 R_AC90_eval Ys Yp Alpha_E_AG
global Theta_inf_dielectric optimization memoring theta_VECTOR_ISGT days
global thermal_memory theta_VECTOR_global thermal_memory_preload TEMPERATURE_D_1 preloading
%% ************************************************************************
%% Initial temperature of cable (dielectric losses are added later!)
%% ************************************************************************
Tamb=273.15+18; % Ambient temperature in Kelvin (K)
% p10 IEC853-2: it has been assumed that the voltage had been applied for
% a sufficiently long time for the conductor temperature rise due to
% dielectric loss to have reached a steady state. The total temperature
% rise of the conductor above ambient is then the sum of the steady state
% temperature rise due to the dielectric loss (as given by IEC 60287) and
% the transient temperature via riations due to change of current.
% Tinit=Tamb; % Initial temperature of cable in Kelvin (K)
% Tinit is updated in the end of the script!!!
% Here, Tinit = Tamb but it is also neccesary to add the temperature rise
% due to dielectric losses. The dielectirc losses and temperature rise
% could be calculated by IEC 60287. But to do that , thermal parameters of
% cable should be firstly calculated. (see next parts of the script)
% Temperature rise from dielectric losses is be added in the end of script
%% ************************************************************************
%% System constants
%% ************************************************************************
U_0=225e3/sqrt(3) ; % Phase to neutral voltage in V, 130 kV
f=50; % Frequency in Hz
w=2*pi*f; % Network pulsation
%% ************************************************************************
%% Cable characteristics and geometry
%% ************************************************************************
Lcable=1 ; % cable length (m)
%--------------------------------------------------------------------------
% Core / conductor [for 1 core of three-core cable]
%--------------------------------------------------------------------------
D_ame =38.2e-3; % core diameter (m)
rame=D_ame/2; % core radius (m)
% S_ame_pertes=pi*rame^2; % actual conductor cross-section (m2)
S_ame_pertes=1000e-6; % nominal conductor cross-section (m2)
% Note: We use 1000 mm2 instead of 1146 mm2 because the rest is air
% (which does not transmit current). This is also done in the COMSOL model.
%--------------------------------------------------------------------------
% Semiconductor 1 (Semiconducting compound)
%--------------------------------------------------------------------------
Thick_SC1=1.5e-3; % thickness of semiconductor (m)
D_SC1_int=D_ame; % internal diameter (m)
D_SC1_ext=D_ame+2*Thick_SC1; % external diameter (m)
%--------------------------------------------------------------------------
% Insulation (XLPE)
%--------------------------------------------------------------------------
Thick_Isol=22e-3; % thickness of insulation (m)
D_Isol_int=D_SC1_ext; % internal diameter (m)
D_Isol_ext=D_SC1_ext+2*Thick_Isol; % external diameter (m)
%--------------------------------------------------------------------------
% Semiconductor 2 (Semiconducting compound)
%--------------------------------------------------------------------------
Thick_SC2=1.5e-3; % thickness of semiconductor(m)
D_SC2_int=D_Isol_ext; % internal diameter (m)
D_SC2_ext=D_Isol_ext+2*Thick_SC2; % external diameter (m)
%--------------------------------------------------------------------------
% Sealing mats (Semiconducting water-swellable tapes)
%--------------------------------------------------------------------------
Thick_Matelas=1e-3; % thickness of sealing mats(m)
D_Matelas_int=D_SC2_ext; % internal diameter (m)
D_Matelas_ext=D_SC2_ext+2*Thick_Matelas; % external diameter (m)
%--------------------------------------------------------------------------
% Metallic screen (Lead alloy)
%--------------------------------------------------------------------------
Thick_Pb_Screen=2.5e-3; % thickness of sreen(m)
D_Pb_Screen_int=D_Matelas_ext; % internal diameter (m)
D_Pb_Screen_ext=D_Matelas_ext+2*Thick_Pb_Screen; % external diameter (m)
D_moy_screen=D_Matelas_ext+Thick_Pb_Screen; % average diameter (m)
%--------------------------------------------------------------------------
% Internal protection sheath (Semiconducting PE)
%--------------------------------------------------------------------------
Thick_Gaine_PH=3e-3; % sheath thickness (m)
D_Gaine_PH_int=D_Pb_Screen_ext; % internal diameter (m)
D_Gaine_PH_ext=D_Gaine_PH_int+2*Thick_Gaine_PH; % external diameter (m)
%--------------------------------------------------------------------------
% Filler (Plastic shaped elements)
%--------------------------------------------------------------------------
Dext_filler = 228e-3; % D_Armure_Matelas_int (below)
%--------------------------------------------------------------------------
% Armour mats (PP yarn)
%--------------------------------------------------------------------------
Thick_Armure_Matelas=0.5e-3; % thickness of armour mats (m)
D_Armure_Matelas_int=228e-3; % internal diameter (m)
D_Armure_Matelas_ext=D_Armure_Matelas_int...% external diameter (m)
+2*Thick_Armure_Matelas;
%--------------------------------------------------------------------------
% Armour (Steel)
%--------------------------------------------------------------------------
Ep_arm=6e-3; % armour thickness (m)
Dfil_Arm=Ep_arm; % diameter of round wire (m)
Dint_arm=D_Armure_Matelas_ext; % internal diameter (m)
Dext_arm=Dint_arm+2*Ep_arm; % external diameter (m)
Dmoy_Arm=Dint_arm+1*Ep_arm; % average diameter (m)
Nfil_Arm=99; % wire number in armour (-)
%--------------------------------------------------------------------------
% Outer serving (PP yarn)
%--------------------------------------------------------------------------
Ep_Gaine_Ext=5e-3; % thickness of outer serving (m)
Dgaine_ext_int=Dext_arm; % internal diameter (m)
Dgaine_ext_ext=Dext_arm+2*Ep_Gaine_Ext; % external diameter (m)
%% ************************************************************************
%% Physical parameters of the materials (electrical and thermal)
%% ************************************************************************
%--------------------------------------------------------------------------
% Aluminium
%--------------------------------------------------------------------------
Rho_0_E_Al=2.8264e-8; % electrical resistivity of Alu at 20°C % (p30 Table1
% 60287-1-1)
Alpha_E_Al=4.03e-3; % temperature coefficient(p30 Table1 IEC 60287-1-1)
Cth_Al=2.5e6; % volumetric specific heat (J/m3.K)(p77 IEC 60853-2)
%--------------------------------------------------------------------------
% Semiconductor1 (Semiconducting compound)
%--------------------------------------------------------------------------
% Note: p19 IEC 60287-2-1:for thermal calculations semi-conducting layers
%(including metallized carbon paper tapes) are considered as part of
% the insulation.
%--------------------------------------------------------------------------
% Insulation - XLPE
%--------------------------------------------------------------------------
Rho_th_isol=3.5; % thermal resistivity (K.m/W) p75 IEC 60853-2
Cth_Isol=2.4e6; % volumetric specific heat(J/m3.K) p75 IEC 60853-2
eps_r=2.5; % relative permittivity (p64 CEI 60287-1-1)for
% cables greater than 18/30 (36) kV cables (unfilled)
tan_delta=0.001; % loss factors (p64 CEI 60287-1-1)for
% cables greater than 18/30 (36) kV cables(unfilled)*
%*Safe values at maximum permissible temperature, applicable to the highest
% voltages normally specified for each type of cable.
%--------------------------------------------------------------------------
% SemiConductor2 (Semiconducting compound)
%--------------------------------------------------------------------------
% Note: p19 IEC 60287-2-1:for thermal calculations semi-conducting layers
%(including metallized carbon paper tapes) are considered as part of
% the insulation.
%--------------------------------------------------------------------------
% Water-swellable tape
%--------------------------------------------------------------------------
% Note: p19 IEC 60287-2-1:for thermal calculations semi-conducting layers
%(including metallized carbon paper tapes) are considered as part of
% the insulation.
% Note: Swelling tape is assumed to be a part of the insulation since its
% thickness is small, and it is assumed to have the same thermal
% resistivity. (p6 Nielsen, Applied Sciences 2019)
%--------------------------------------------------------------------------
% Screen (lead alloy)
%--------------------------------------------------------------------------
Cth_ecran_Pb=1.45e6; % 1.7e6; % volumetric specific heat(J/m3.K)
Rho_0_E_Pb=21.4e-8; % electrical resistivity of lead alloy at20°C
% (p30 Table1 60287-1-1)
Alpha_E_Pb=4.0e-3; % temperature coefficient of lead alloy at20°C
% (p30 Table1 IEC 60287-1-1)
% Note: Metallic layers are neglected, as their thermal resistance is
% negligible compared with that of polycomposite materials.
% (p6 Nielsen, Applied Sciences 2019)
% Another note: the thermal resistance of the metallic parts in the cable,
% even though not equal to zero, is so small that it is usually neglected
% in rating computations (p34 Anders Rating of electric power cables 1997)
% (IEC does not provide instructions about thermal resistances of
% metallic parts?)
%--------------------------------------------------------------------------
% Internal sheath (Semiconducting PE)
%--------------------------------------------------------------------------
% Note: p19 IEC 60287-2-1:for thermal calculations semi-conducting layers
%(including metallized carbon paper tapes) are considered as part of
% the insulation.
Rho_th_gaine_PH=3.5; % thermal resistivity (K.m/W)
% p75 IEC 60853-2
Cth_Gaine_PH=2.4e6;%1.7e6; % volumetric specific heat(J/m3.K)
% p75 IEC 60853-2
%--------------------------------------------------------------------------
% Armour matelas (polypropylene yarn)
%--------------------------------------------------------------------------
Rho_th_armure_matelas=5.5; % thermal resistivity (K.m/W)
% p55 IEC 60287-2-1
Cth_armure_matelas=1.8e6; % volumetric specific heat(J/m3.K)
% not given in IEC. taken from p8 Nielsen Applied Sciences 2019
%--------------------------------------------------------------------------
% Filler (Plastic shaped elements)
%--------------------------------------------------------------------------
% Kappa_filler=Kappa_Gaine_tri_int*0.5; % 0.5:article EWTEC 2017 (EMODI)
%CHB -> Plutot revue Elsevier renewable energy
% Cth_filler=Cth_armure_matelas*0.5;
%--------------------------------------------------------------------------
% Armour (steel or stainless steel)
%--------------------------------------------------------------------------
% Cp_th_acier=1.006*500;
% Rho_m_acier=7900;
Cth_Armure=3.8e6; % volumetric specific heat(J/m3.K)
% p77 TableE2 IEC 60853-2
% Note: Metallic layers are neglected, as their thermal resistance is
% negligible compared with that of poly-composite materials.
% (p6 Nielsen, Applied Sciences 2019)
% Another note: the thermal resistance of the metallic parts in the cable,
% even though not equal to zero, is so small that it is usually neglected
% in rating computations (p34 Anders Rating of electric power cables 1997)
stainless_steel=1; % stainless steel is 1; steel is 0
% the resistance of the armour at 20°C (?/m). p30 IEC 60287-1-1
if stainless_steel==1
Rho_E_AG=70e-8; % Stainless Steel
else
Rho_E_AG=13.8e-8; % Steel
end
% Note that Temperature coefficient of resistance is negligible for
% stainless steel; thus,no correction is necessary in Sub-clause 4.4.2.
% (see p77 IEC853-2)
% temperature coefficient of steel at20°C (p30 Table1 IEC 60287-1-1)
if stainless_steel==1
Alpha_E_AG=0; % Stainless Steel
else
Alpha_E_AG=4.5e-3; % Steel
end
%--------------------------------------------------------------------------
% Outer serving (polypropylene yarn)
%--------------------------------------------------------------------------
Rho_th_Gaine_ext=5.5; % thermal resistivity (K.m/W)
% Rho_th_Gaine_ext=10; % thermal resistivity (K.m/W) Nielsen 2019
% p55 IEC 60287-2-1
Cth_Gaine_ext=1.8e6; % volumetric specific heat(J/m3.K)
% p77 TableE2 IEC 60853-2
%--------------------------------------------------------------------------
% Soil
%--------------------------------------------------------------------------
% Kappa_Sol=1/0.7;
% Cp_sol=2e6;
Rho_sol=0.8; % thermal resistivity (K.m/W) p73 Table D1 IEC 60853-2
Diff_sol=0.6e-6; % diffusivity of soil m2/s p73 Table D1 IEC 60853-2
De=Dgaine_ext_ext; % external diameter of cable
Depth=1.75; % depth of cable burying m
Depth=1.87; % Value used in simulations, m
%% ************************************************************************
%% Calculation of equivalent thermal restistance (T1,T2,T3) for a model of
%% three-core cable
%% ************************************************************************
% 2 Calculation of thermal resistances
% 2.1 Thermal resistance of the constituent parts of a cable, T1, T2 and T3
% Important note (p19 IEC 60287-2-1):
% Where screening layers are present, for thermal calculations metallic
% tapes are considered to be part of the conductor or sheath while
% semi-conducting layers (including metallized carbon paper tapes) are
% considered as part of the insulation. The appropriate component
% dimensions must be modified accordingly.
% 4.2.1.2 Three-core cables (% p25 IEC 60853-2)
% The three-core cable is replaced by an equivalent single-core
% construction % dissipating the same total conductor losses.
% The equivalent single-core conductor has a diameter: dc but firstly
% p22 IEC 60287-2-1
%--------------------------------------------------------------------------
% T1,the thermal resistance between one conductor and the sheath:
% 2.1.1.5 SL and SA type cables
% The thermal resistance T1 is calculated in the same way as for
% single-core cables (see next 2.1.1.1).
% 2.1.1.1 Single-core cables (p19 IEC 60287-2-1)
% The thermal resistance between one conductor and the sheath T1 is given
% by:
T1 = (Rho_th_isol/(2*pi))*log(1+2*((D_Matelas_ext/2)-rame)/(2*rame));
%--------------------------------------------------------------------------
% T2, The thermal resistance of fillers and bedding under the armour:
% 2.1.2.2 SL and SA type cables (p.27 IEC 60287-2-1)
% Firstly, it is neccesary to calculate G the geometric factor given in
% figure 6 (p71 IEC 60287-2-1)
% Thickness between sheaths and armour
Ep_Gaine_Armure = Thick_Gaine_PH+Thick_Armure_Matelas;
% Ep_Gaine_Armure = Thick_Armure_Matelas;
% dans Dorison, au début du livre (p36, tome 2) facteur X p36 Dorison : contradiction avec la norme (p70 IEC 60287-2-1)
% Verification avec le dessin 6 (p70 IEC 60287-2-1)
% Thickness of material between sheaths and armour expressed as a fraction
% of the outer diamter of the sheath
% facteur_X = Ep_Gaine_Armure/D_Pb_Screen_ext;
facteur_X = Ep_Gaine_Armure/D_Armure_Matelas_int;
% facteur_X = Ep_Gaine_Armure/D_Gaine_PH_ext;
if (0<facteur_X && facteur_X<=0.03)
facteur_G = 2*pi*(0.00022619+2.11429*facteur_X-20.4762*facteur_X^2); % upper curve p51 CEI 60287-2-1
% facteur_G = 2*pi*(0.000202380+2.03214*facteur_X-21.6667*facteur_X^2); % Nielsen 2019
elseif (0.03<facteur_X && facteur_X<=0.15)
facteur_G = 2*pi*(0.0142108+1.17533*facteur_X-4.49737*facteur_X^2+... % upper curvep 51 CEI 60287-2-1
10.6352*facteur_X^3);
% facteur_G = 2*pi*(0.0126529+1.101*facteur_X-4.56104*facteur_X^2+... % Nielsen 2019
% 11.5093*facteur_X^3);
end
% Finally, the thermal resistance of fillers and bedding under the armour:
T2 = (Rho_th_gaine_PH/(6*pi))*facteur_G;
% T2 = (Rho_th_armure_matelas/(6*pi))*facteur_G;
%--------------------------------------------------------------------------
% T3,the thermal resistance of outer covering (serving):
% 2.1.3 Thermal resistance of outer covering (serving) T3(p27 IEC60287-2-1)
% The external servings are generally in the form of concentric layers and
% the thermal resistance T3 is given by:
T3 = (Rho_th_Gaine_ext/(2*pi))*log(1+2*(Ep_Gaine_Ext)/(Dext_arm));
%--------------------------------------------------------------------------
% T4, the external thermal resistance:
% 2.2 External thermal resistance T4 (p29 IEC 60287-2-1)
% 2.2.2 Single isolated buried cable T4(p31 IEC 60287-2-1)
u = 2*Depth/De;
T4 = (Rho_sol/(2*pi))*log(u+sqrt(u^2-1));
%% ************************************************************************
%% Calculation of the equivalent thermal capacitances for a model of
%% three-core cable.
%% ************************************************************************
% 4.2.1.2 Three-core cables (p25 IEC 60853-2):
% The three-core cable is replaced by an equivalent single-core
% construction dissipating the same total conductor losses:
% external diameter of dielectricc. The same value of diameter over
% dielectric (under the sheath) as for the three-core cable
Di=D_Matelas_ext;
% Important note (p25 IEC 60853-2):
% Where screening layers are present: for thermal calculations metallic
% tapes are considered to be part of the conductor or sheath while
% semi-conducting layers (including metallized carbon paper tapes) are
% considered as part of the insulation. The appropriate component
% dimensions should be modified accordingly.
% T1 is one-third of the value for one of the cores of the three-core cable
% as given in IEC Publication 287,
% The equivalent single-core conductor has a diameter (dc):
dc=(Di*exp(-(2*pi*T1/3)/Rho_th_isol));
pVW = (1/(2*log(Di/dc)))-(1/((Di/dc)^2-1));
% Alternatives:
% % (Dorison, dès % la page 136 (pour câble tri-écranté et p156).
% % Coefficient de répartition p (long) et p* (court)
%
% % p156: "r1 et r2 sont les rayons de l'âme et sur isolant" ->
% % interprétation T1: rayon intérieur du diélectrique et T2 rayon extérieur
% % du diélectrique
% %
% r2_sur_r1=(D_CU_Screen_int/2)/rame;% internal diameter écran cuivre
% %(contenant âme+ 2 écrans semi-conducteurs + diélectrique)/rayon de l'âme
%
% %pVW=1/((r2_sur_r1)^2-1)+1/(2*log(r2_sur_r1)); % p
%
% % p (Dorison, p156) (corrigé pendant la révision)
% pVW=-1/((r2_sur_r1)^2-1)+1/(2*log(r2_sur_r1));
%
% % Alternative calculation of pVW (p10 Nielsen, Applied Sciences 2019)
% Dc_star=D_Isol_ext*e((-2*pi*T1)/(Rho_th_isol*3));
% pVW=1/(2*log(D_Isol_ext/Dc_star))-1/((D_Isol_ext/Dc_star)^2-1);
%--------------------------------------------------------------------------
% (p25 IEC 60853-2):
% Thermal capacitances are calculated on the following assumptions:
% a) The actual conductors are considered to be completely inside the
% diameter of the equivalent single-core conductor, the remainder of the
% equivalent conductor being occupied by insulation.
% (p27 IEC 60853-2):
% b) The space between the equivalent single-core conductor and the sheath
% is considered to be completely occupied by insulation
% (p27 IEC 60853-2):
% The factor p is then calculated using the dimensions of the equivalent
% single-core cable and is applied to the thermal capacitance of the
% insulation based on assumption b) above.
% thermal capacitance of conductor (or equivalent single-core conductor of
% a three-core cable, see Sub-clause 4.2.1.2)(p25 IEC 60853-2):
Qc = pi*rame^2*Cth_Al;
Qc = S_ame_pertes*Cth_Al;
% Qc = pi*(dc/2)^2*Cth_Al;
% Qc is for equivalent single-core cable ie dc is used instead of rame
% total thermal capacitance of dielectric per conductor (or equivalent
% single-core conductor of a three-core cable, see Sub-clause 4.2.1.2)
Qi = Cth_Isol*(pi/4)*(D_Pb_Screen_int^2-(2*rame)^2);
% Total thermal capacitance of the screen
Qs=Cth_ecran_Pb*(pi/4)*(D_Pb_Screen_ext^2-D_Pb_Screen_int^2);
% Total thermal capacitance of the filler
Qf = Cth_filler*((pi/4*(D_Armure_Matelas_ext^2))-...
3*(pi/4*(D_Gaine_PH_ext^2)));
% Total thermal capacitance of the armour
Qa = Cth_Armure*Nfil_Arm*pi*Ep_arm^2/4;
% Total thermal capacitance of the outer serving
Qj = Cth_Gaine_ext*(pi*(Dgaine_ext_ext^2-Dgaine_ext_int^2)/(4));
% Note that IEC does not give precise instructions on how to calculate
% thermal capacitances of the equivalent single-core conductor.
% Nevertheless,G. Anders writes in "Transient analysis of 3-core SL-type
% submarine cables with jacket around each core", Jicable 2015:
% The thermal capacitances of the (equivalent single-core) conductor, the
% insulation and the screen are equal to three times those of the SL cable
% Therefore we redefine:
Qc=3*Qc;
Qi=3*Qi;
Qs=3*Qs;
%% ************************************************************************
%% Loss factors Lambda_1 and Lambda_2
%% ************************************************************************
% 2 Calculation of losses
% 2.1 AC resistance of conductor (p13 IEC 60287-1-1)
% 2.1.1 DC resistance of conductor at 20 degC
% Note 287-1-1:The value of R_DC20 shall be derived directly from IEC 60228.
% Where the conductor size is outside the range covered by IEC 60228,
% the value of R_DC20 may be chosen by agreement between manufacturer and
% purchaser:
% R_DC20=Rho_0_E_Al*Lcable/S_ame_pertes;
% R_DC20 for 1000 mm2 (p9 IEC 60228)?/m!
R_DC20=0.0291e-3;
% Note that IEC 60228 provides R_DC20 in ?/km and IEC 60287-1-1 uses ?/m !
% The d.c. resistance per unit length of the conductor (?/m)at its maximum
% operating temperature ? is given by (p14 IEC 60287-1-1):
R_DC90=R_DC20*(1+Alpha_E_Al*(90-20));
% Type of conductor: Aluminium. Round, stranded.(p31 Table 2 IEC 60287-1-1)
ks=1;
% Although there are no accepted experimental results dealing specifically
% with the coefficient kp for aluminium conductors, it is recommended that,
% for stranded aluminium conductors, the values given for similar copper
% conductors are used. (p31 Table 2 IEC 60287-1-1)
kp=1;
% 2.1.2 Skin effect factor Ys(p14 IEC 60287-1-1)
Xs=sqrt(8*pi*f/R_DC90*ks*1e-7);
% The above formula is accurate providing xs does not exceed 2,8, and
% therefore applies to the majority of practical cases.
if Xs>=2.8
error('Xs>=2.8')
end
% The skin effect factor ys is given by (p14 IEC 60287-1-1):
Ys=Xs^4/(192+0.8*Xs^4);
% The Ys formula is accurate providing Xs does not exceed 2.8, and
% therefore applies to the majority of practical cases.
% 2.1.4 Proximity effect factor yp for three-core cables and for three
% single-core cables (p15 IEC 60287-1-1):
% 2.1.4.1 Circular conductor cables
Xp=sqrt(8*pi*f/R_DC90*kp*1e-7);
% The above formula is accurate providing xp does not exceed 2,8, and
% therefore applies to the majority of practical cases. (p15 IEC 60287-1-1)
if Xp>=2.8
error('Xp>=2.8')
end
% Distance between conductor axes (m) (equal to the external diameter of
% internal sheaths):
s=D_Gaine_PH_ext;
% the distance between the axis of a conductor and the cable centre (m)
c=s/sqrt(3);
% Only if Xp<2.8 else Xp=>2.8?
% The proximity effect factor is given by:
Yp=(Xp^4/(192+0.8*Xp^4))*(D_ame/s)^2*(0.312*((D_ame/s)^2)...
+1.18/(Xp^4/(192+0.8*Xp^4)+0.27));
% Yp=(Xp^4/(192+0.8*Xp^4))*(dc/s)^2*(0.312*((dc/s)^2)...
% +1.18/(Xp^4/(192+0.8*Xp^4)+0.27));
% disp('Yp formula is changed. dc is used instead of D_ame')
% 2.3.10 Cables with each core in a separate lead sheath (SL type)&armoured
% (p24 IEC 60287-1-1)
% For a three-core cable of which each core has a separate lead sheath ???
% is zero and the loss factor for the sheaths is given by:
% d, mean diameter of sheath or screen
d=D_moy_screen;
Xcable=2*w*1e-7*log(2*s/d);
% 2.1 AC resistance of conductor (p13 IEC 60287-1-1)
% The a.c. resistance per unit length of the conductor at its maximum
% operating temperature is given by the following formula
R_AC90_eval=R_DC90*(1+Ys+Yp);
% 2.2 Dielectric losses (applicable to a.c. cables only) p16 IEC 60287-1-1
% the capacitance per unit length (F/m)
% C=eps_r*1e-9/(18*log(D_Matelas_ext/(2*rame)));
% It seems the formula with D_Isol_ext/(D_Isol_int is correct:
C=eps_r*1e-9/(18*log(D_Isol_ext/(D_Isol_int)));
% The dielectric loss per unit length in each phase is given by p16 IEC
% 60287-1-1:
Wd=C*w*(U_0)^2*tan_delta;
% The cross section of screen
S_ecran=(pi/4)*(D_Pb_Screen_ext^2-D_Pb_Screen_int^2);% Check no data in IEC
% Recran=Rho_0_E_Pb*Lcable/S_ecran;
% Rs=Recran; % ou Rgaine? the same for Xcable above
% the resistance of the cable sheath or screen at 20 °C (?/m).
% p17 IEC 60287-1-1 but formula is taken from Cablesizer
% https://www.cableizer.com/documentation/R_sh/ :
Rso=Rho_0_E_Pb/S_ecran;
% the resistance of the armour at 20°C:
% p25 IEC 60287-1-1 but formula is assumed based on the resistance of the
% cable sheath or screen at 20 °C (?/m).
Rao=Rho_E_AG/(Nfil_Arm*pi*(Dfil_Arm/2)^2); % Hypothesis based on SEMREV
%--------------------------------------------------------------------------
% While cycle for calculation of lambda1 and lambda2
iteration=0; % iteration number (used for debugging only)
Error=1; % any arbitrary number >=1e-3
Tmax=90+273.15; % maximal operating temperature of conductor
% First guess of steady-state current (needed for calculation of sheath and
% armour temperature in the while cycle):
I_adm_NEW=I_admissible(Wd,T1,T2,T3,T4,3,R_AC90_eval,0.001,0.001,Tmax,Tamb);
I_adm_OLD=I_adm_NEW;
while (Error>=1e-3)
% maximum operating temperature of sheath or screen is given by
% p17 IEC 60287-1-1:
Temp_sc = 90-(I_adm_NEW^2*R_AC90_eval+0.5*Wd)*T1;
% The resistance of the sheath or screen at its maximum operating
% temperature is given by p17 IEC 60287-1-1:
Rs = Rso*(1+Alpha_E_Pb*(Temp_sc-20));
% Ratio of the total losses in sheath respectively to the total
% conductor losses (or losses in one sheath to the losses in one
% conductor)(p24 IEC 60287-1-1)
Lambda1=(Rs/R_AC90_eval)*(1.5/(1+(Rs/Xcable)^2));
% NB: Here Lambda1 is actually Lambda1'. In general Lambda1 is
% calculated as Lambda1=Lambda1'+Lambda1'' (see p16 IEC 60287-1-1)but
% in accordance with p24 IEC 60287-1-1 Lambda1''==0. Thus, we assume
% that Lambda1=Lambda1';
% The maximum operating temperature of the armour is given by p25 IEC
% 60287-1-1:
Temp_AG = 90-((I_adm_NEW^2*R_AC90_eval+0.5*Wd)*T1+(I_adm_NEW^2*...
R_AC90_eval*(1+Lambda1)+Wd)*3*T2);
% The resistance of the armour at its maximum operating temperature is
% given by p25 IEC 60287-1-1:
Rarm=Rao*(1+Alpha_E_AG*(Temp_AG-20));
% 2.4.2.5 SL type cables (p29 CEI 60287-1-1)
% Where the armour is over a SL type cable, the screening effect of the
% sheath currents reduces the armour loss. The formula for ?2 given in
% 2.4.2.3.1 shall be multiplied by the factor
Lambda1_dash=(Rs/R_AC90_eval)*1/(1+(Rs/Xcable)^2);
Factor=(1-(R_AC90_eval/Rs)*Lambda1_dash);
% where Lambda1_dash is obtained from 2.3.1. % 2.3.1 is for Two single-core cables, and three single-core cables
% Ratio of the total losses in armour respectively to the total
% conductor losses (or losses in one sheath or armour to the losses
% in one conductor) p28 CEI 60287-1-1
Lambda2=1.23*(Rarm/R_AC90_eval)*(2*c/Dmoy_Arm)^2*1/(1+(2.77*Rarm...
*1e6/(w))^2);
Lambda2=Factor*Lambda2;
if stainless_steel==1
Lambda2=0; % if stainless steel is used
end
% Steady-state rating of cable
I_adm_NEW=I_admissible(Wd,T1,T2,T3,T4,3,R_AC90_eval,Lambda1,Lambda2,...
Tmax,Tamb);
Error=abs(I_adm_NEW-I_adm_OLD);
iteration=iteration+1;
I_adm_OLD=I_adm_NEW;
end
if stainless_steel==1
disp('Stainless steel is assumed for armour. Lambda2=0 ie no losses in armour')
end
% Final (eventual) steady-state current
I_adm=I_admissible(Wd,T1,T2,T3,T4,3,R_AC90_eval,Lambda1,Lambda2,Tmax,Tamb);
%% ************************************************************************
%% Equivalent thermal resistances and capacitances
%% ************************************************************************
% Appendix A provides a detailed method for reducing a multiple thermal
% circuit to a two-cell circuit.
% 4.2.2.2 Representation of common types of cable (p27. IEC 60853-2)
% (f) Armoured cables with each core in a separate lead sheath (SL type)
% (p31 IEC 60853-2)
% Note T1 is the thermal resistance of the equivalent single-core cable.
% T1, is one-third of the value for one of the cores of the three-core
% cable as given in IEC Publication 287, i.e TA=T1/3
TA=T1/3;
% qs ratio: losses in conductors and screens/losses in conductors
qs=Lambda1+1;
% Not specified in IEC.
% Source: (https://www.cableizer.com/documentation/q_s/)
% Tf,thermal resistance of filling i.e.thermal resistance of filling
% between cores of SL type cable
Tf=T2; % no reference on how Tf is calculated thus assumed equal to T2
Tf=0.020344899802;
% qa, ratio: losses in (conductors + sheaths + armour)/losses in conductors
qa=Lambda1+Lambda2+1;
% Not specified in IEC standard.
% Source qa=Lambda1+Lambda2+1(https://www.cableizer.com/documentation/q_a/)
% Another source confirming it is ANDERS, George J, George, GEORGALLIS.
% "Transient analysis of 3-core SL-type submarine cables with jacket around
% each core" 2015
TB=qs*Tf+qa*T3; % intiial IEC formula
% TB=3*(qs*Tf+qa*T3); % Anders modification in "Transient analysis of
% % 3-core SL-type submarine cables with jacket around
% % each core" 2015
% disp('Anders modification of TB is used')
QA=Qc+pVW*Qi; % Attention Qc and Qi of equivalent conductor or not
% QA=3*(Qc+pVW*Qi); % Attention Qc and Qi of equivalent conductor or not
% disp('QA is multiplied by 3. Hypothesis')
Ts=0;
QB=(1-pVW)*Qi+((Qs+0.5*Qf)/qs)+((qa*T3/(qs*T2+qa*T3))^2)*((0.5*Qf/qs)+...
((Qa+Qj)/qa)); % intiial IEC formula
% QB=(1-pVW)*Qi+((Qs+Qf/6)/qs)+((3*qa*T3/(qs*T2+3*qa*T3))^2)*(((Qf/6)/qs)+...
% ((Qa+Qj)/(3*qa))); % Anders modification in "Transient analysis of
% % % 3-core SL-type submarine cables with jacket around
% % % each core" 2015
% disp('Anders modification of QB is used')
%% Partials transients calculation
% 4.2.3 Calculation of cable partial transient (p33 IEC 853-2)
% The transient response of a cable circuit to a step function of load
% current, considered in isolation, that is with the right-hand pair of
% terminals in Figure 2 short-circuited, is found as follows:
M0=0.5*(QA*(TA+TB)+QB*TB);
N0=QA*TA*QB*TB;
a=(M0+sqrt(M0^2-N0))/N0;
b=(M0-sqrt(M0^2-N0))/N0;
Ta=(1/(a-b))*((1/QA)-b*(TA+TB));
Tb=TA+TB-Ta;
%--------------------------------------------------------------------------
% Ensuring the column vector
if (size(Iech_load,2) > size(Iech_load,1)) % If number of columns > number of rows
Iech_load = Iech_load';
end
if isempty(optimization)
% extract all values of load
I_ech=Iech_load;
elseif optimization==1 %if optimization is in progress
% extract last 96 values of load for J day
I_ech=Iech_load(days*96-95:end,1);
else % something wrong
error('Check the optimization status');
end
if memoring==1
% extract J day and the rest (zero arrays)
I_ech=Iech_load(days*96-95:end,1);
end
% Duration of one current value in seconds
duration=900;
%-------------------------Initial code-------------------------------
tstart=900; % choosing tstart=3600 instead of 0.001 allows avoiding
% the problem with temperature_correction_on_losses.m later
% timestep=5; % time step, 5 seconds
timestep=900; % time step, 1 hour = 3600 s
tend=duration;
t=tstart:timestep:tend*length(I_ech); % time IEC in seconds
t=t';
Decalage_down=zeros(duration/timestep,1); % time offset for down current
%--------------------------------------------------------------------------
if preloading==1
I_ech=I_ech(1:35040);
end
%--------------------------------------------------------------------------
% Superposition principle is applied in this "for cycle":
for k=1:length(I_ech) % for each current value
if memoring==1 % if memoring status of IEC60853_2 is switched on
if k<=(days*96+1)-(days*96-95) % continue k until 96 values are passed
% Define the time offset for starting point
Decalage=zeros(duration/timestep*k-duration/timestep,1);
% Define I as k-th current value from I_ech
I=I_ech(k);
% Find cable and environment trainsient as well as attainement factor
[Theta_ct,Theta_Et,Alpha_t]=losses_partial_transients(I);
% 4.4 Calculation of the complete temperature transient (p43 IEC853-2)
% After calculating separately the two partial transients and the
% conductor to cable surface attainment factor the total transient rise
% above an ambient temperature is obtained:
% a) for buried configurations by simple addition of the cable and
% modified environment partial transients;
% the complete temperature transient (p43 IEC853-2)
Theta=Theta_ct+Alpha_t.*Theta_Et;
if k==1 % first current value in load profile
% Shift the complete temperature transient for 1 hour ahead and
% make it negative
Theta_down=-[Decalage_down;Theta];
% save only studied interval
Theta_down=Theta_down(1:length(t),1);
% save to matrix for down step current
theta_VECTOR=Theta+Theta_down;
else % second and next current value in load profile
% Shift the complete temperature transient (up step) in time to its
% start point
Theta=[Decalage;Theta];
% save only studied interval
Theta=Theta(1:length(t),1);
% Shift the complete temperature transient for 1 hour ahead and
% make it negative
Theta_down=-[Decalage_down;Theta];
% save only studied interval
Theta_down=Theta_down(1:length(t),1);
% save to matrix for down step current
theta_VECTOR=theta_VECTOR+Theta+Theta_down;
end % end of if k == 1
else % if k>(days*24+1)-(days*24-23)
break % stop for cycle and go to next section of IEC60853_2.m
end % end of if k<=(days*96+1)-(days*96-95)
else % if memoring = []
% Define the time offset for starting point
Decalage=zeros(duration/timestep*k-duration/timestep,1);
% Define I as k-th current value from I_ech
I=I_ech(k);
% Find cable and environment trainsient as well as attainement factor
[Theta_ct,Theta_Et,Alpha_t]=losses_partial_transients(I);
% 4.4 Calculation of the complete temperature transient (p43 IEC853-2)
% After calculating separately the two partial transients and the
% conductor to cable surface attainment factor the total transient rise
% above an ambient temperature is obtained:
% a) for buried configurations by simple addition of the cable and
% modified environment partial transients;
% the complete temperature transient (p43 IEC853-2)
Theta=Theta_ct+Alpha_t.*Theta_Et;
if k==1 % first current value in load profile
% Shift the complete temperature transient for 1 hour ahead and
% make it negative
Theta_down=-[Decalage_down;Theta];
% save only studied interval
Theta_down=Theta_down(1:length(t),1);
% save to matrix for down step current
theta_VECTOR=Theta+Theta_down;
else % second and next current value in load profile
% Shift the complete temperature transient (up step) in time to its
% start point
Theta=[Decalage;Theta];
% save only studied interval
Theta=Theta(1:length(t),1);
% Shift the complete temperature transient for 1 hour ahead and
% make it negative
Theta_down=-[Decalage_down;Theta];
% save only studied interval
Theta_down=Theta_down(1:length(t),1);
% save to matrix for down step current
theta_VECTOR=theta_VECTOR+Theta+Theta_down;
end % end of if k == 1
end % end of if memoring==1
end % go to next k or end of for cycle
%% Thermal memory consideration in the optimization problem
% if optimization==1
% if isempty(thermal_memory) % no thermal memory
% %do nothing
% else % there is a thermal memory
% theta_VECTOR=theta_VECTOR+thermal_memory(days*96-95:days*96,1); % J day
% end
% end
%% Correction of R as a function of T
Error=1e-3;
MAE=100;
if preloading==1
I_ech=Iech_load;
end
% The conductor steady-state temperature rise above ambient:
% Source: 1.4.1 Buried cables where drying out of the soil does not occur
% % or cables in air (1.4.1.1 AC cables) (p10 IEC 60287-1-1)
Theta_inf_dielectric=(0.^2.*R_AC90_eval+0.5*Wd)*T1+(0.^2.*R_AC90_eval*(1+Lambda1)+...
Wd)*3*T2+(0.^2.*R_AC90_eval*(1+Lambda1+Lambda2)+Wd)*3*(T3+T4);
while (MAE>=Error) % Mean absolute error in degrees of C
if preloading==1
I_ech=Iech_load;
end
[theta_VECTOR,MAE]=temperature_correction_on_losses_new(theta_VECTOR,f,ks,kp,D_ame,s,I_ech,duration,...
timestep,Decalage_down,Rso,Rao,Xcable,c,Dmoy_Arm,w,stainless_steel);
end
theta_VECTOR_ISGT=theta_VECTOR;
if isempty(optimization)&& isempty(memoring)&& ~isempty(thermal_memory_preload)
theta_VECTOR=theta_VECTOR+thermal_memory_preload; % thermal_memory_preload should be at studied horizon!
end
if optimization==1
if isempty(thermal_memory) % no thermal memory
%do nothing
else % there is a thermal memory
theta_VECTOR=theta_VECTOR+thermal_memory(days*96-95:days*96,1); % J day
end
end
Theta_at=theta_VECTOR;
if memoring==1 % if memoring is activated
theta_VECTOR_global=theta_VECTOR; % prepare the transient for thermal memory
if isempty(thermal_memory) % if thermal_memory does not exist
thermal_memory=theta_VECTOR_global; % save transient as thermal_memory
else % there is already a thermal memory
thermal_memory(days*96-95:end,1)=theta_VECTOR_global+thermal_memory(days*96-95:end,1); % save transients for period [studied day : end of horizon]
end
Theta_at=thermal_memory(days*96-95:end,1);
% if isempty(TEMPERATURE_D_1)
% TEMPERATURE_D_1{1,1}=thermal_memory;
% else
% TEMPERATURE_D_1{end+1,1}=thermal_memory;
% end
end
%% Memoring
% if memoring==1 % if memoring is activated
% theta_VECTOR_global=theta_VECTOR; % prepare the transient for thermal memory
% if isempty(thermal_memory) % if thermal_memory does not exist
% thermal_memory=theta_VECTOR_global; % save transient as thermal_memory
% else % there is already a thermal memory
% thermal_memory(days*96-95:end,1)=theta_VECTOR_global+thermal_memory(days*96-95:end,1); % save transients for period [studied day : end of horizon]
% end
% end % end of if memoring== 1
% if isempty(optimization)&& isempty(memoring)
% MAE=100;
% while (MAE>=Error) % Mean absolute error in degrees of C
% [theta_VECTOR,MAE]=temperature_correction_on_losses_new(theta_VECTOR,f,ks,kp,D_ame,s,I_ech,duration,timestep,Decalage_down);
% end
% if ~isempty(thermal_memory_preload)
% theta_VECTOR=theta_VECTOR+thermal_memory_preload(377*96+1:377*96+days*96,1);
% end
% Theta_at=theta_VECTOR;
% elseif optimization==1 % if optimization is activated
% theta_VECTOR1=thermal_memory;
% theta_VECTOR1(days*96-95:days*96,1)=theta_VECTOR; % theta_VECTOR considers thermal_memory
% MAE=100;
% while (MAE>=Error) % Mean absolute error in degrees of C
% [theta_VECTOR,MAE]=temperature_correction_on_losses_new(theta_VECTOR1(days*96-95:days*96,1),f,ks,kp,D_ame,s,I_ech,duration,timestep,Decalage_down);
% disp('theta_VECTOR1(days*96-95:days*96,1)???')
% theta_VECTOR1(days*96-95:days*96,1)=theta_VECTOR; % theta_VECTOR after correction R(T)
% Theta_at=theta_VECTOR1(1:days*96,1);
% end
% elseif memoring==1 %if memoring is activated
% theta_VECTOR1=thermal_memory(days*96-95:end,1);
% MAE=100;
% while (MAE>=Error) % Mean absolute error in degrees of C
% [theta_VECTOR1,MAE]=temperature_correction_on_losses_new(theta_VECTOR1,f,ks,kp,D_ame,s,I_ech,duration,timestep,Decalage_down);
% Theta_at=theta_VECTOR1;
% end
% end
%% Preparing the temperature profile
% Initial temperature in Kelvin (K) considering dielectric losses
Tinit=Tamb+Theta_inf_dielectric;
% Where Tinit was defined earlier as Tinit=Tamb (see beginning of script)
if isempty(thermal_memory)&& isempty(optimization)&& isempty(memoring)&& isempty(thermal_memory_preload) %if thermal_memory
% Temperature in degC at the beginning of transient t=0
Temperature=Tinit-273.15;
elseif isempty(optimization)&& isempty(memoring)&& ~isempty(thermal_memory_preload)
Temperature=Tinit+thermal_memory_preload(1)-273.15;
elseif optimization==1 %& ~isempty(thermal_memory)
if days>1
Temperature=Tinit+thermal_memory(days*96-96)-273.15;
else % days==1 or other
Temperature=Tinit-273.15;
end
elseif memoring==1 %& ~isempty(thermal_memory)
if days>1
Temperature=Tinit+thermal_memory(days*96-96)-273.15;
else % days==1 or other
Temperature=Tinit-273.15;
end
elseif isempty(optimization) && isempty(memoring)
Temperature=Tinit-273.15;
else% otherwise
if days>1
Temperature=Tinit+thermal_memory(days*96-96)-273.15;
else % days==1 or other
Temperature=Tinit-273.15;
end
end
% Temperature profile in degC (final result)t>0
% Temperature(end+1:length(Theta_at)+1,1)=Theta_at+Tamb+Theta_inf_dielectric-273.15;
Temperature(end+1:length(Theta_at)+1,1)=Theta_at+Tamb+Theta_inf_dielectric-273.15;
disp('Temperature is changed !!!')
disp(['MAE ',num2str(MAE)])
% Note that if Iech_load had a size [24x1] then Temperature will have
% the size [25x1]. The fisrt value of Temperature corresponds to
% temperature value at t=0 and the second one to t=3600s whereas the first
% value of Iech_load corresponds to t=3600s
% The maximal temperature of submarine cable
T_max=max(Temperature);
% if optimization==1