-
Notifications
You must be signed in to change notification settings - Fork 0
/
blpcm.go
573 lines (512 loc) · 15.5 KB
/
blpcm.go
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
//This file is part of EESLISM.
//
//Foobar is free software : you can redistribute itand /or modify
//it under the terms of the GNU General Public License as published by
//the Free Software Foundation, either version 3 of the License, or
//(at your option) any later version.
//
//Foobar is distributed in the hope that it will be useful,
//but WITHOUT ANY WARRANTY; without even the implied warranty of
//MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.See the
//GNU General Public License for more details.
//
//You should have received a copy of the GNU General Public License
//along with Foobar.If not, see < https://www.gnu.org/licenses/>.
/* binit.c */
package eeslism
import (
"fmt"
"math"
"os"
"strconv"
"strings"
)
/* 壁体デ-タの入力 */
func PCMdata(fi *EeTokens, dsn string, pcm *[]*PCM, pcmiterate *rune) {
N := PCMcount(fi)
s := "PCMdata --"
if N > 0 {
*pcm = make([]*PCM, N)
for j := 0; j < N; j++ {
var PCMa = new(PCM)
PCMa.Name = ""
PCMa.Condl = -999.0
PCMa.Conds = -999.0
PCMa.Crol = -999.0
PCMa.Cros = -999.0
PCMa.Ql = -999.0
PCMa.Tl = -999.0
PCMa.Ts = -999.0
PCMa.Tp = -999.0
PCMa.Iterate = false // 収束なしがデフォルト
//PCMa.iterateTemp = false // デフォルトは見かけの比熱だけで収束
PCMa.IterateTemp = true
PCMa.NWeight = 0.5 // 収束計算時の現在ステップ温度の重み係数
PCMa.AveTemp = 'y'
PCMa.DivTemp = 1
PCMa.Ctype = 2 // 二等辺三角形がデフォルト
// パラメータの初期化
PCMa.PCMp.a = -999.0
PCMa.PCMp.B = -999.0
PCMa.PCMp.b = -999.0
PCMa.PCMp.bl = -999.0
PCMa.PCMp.bs = -999.0
PCMa.PCMp.c = -999.0
PCMa.PCMp.d = -999.0
PCMa.PCMp.e = -999.0
PCMa.PCMp.f = -999.0
PCMa.PCMp.omega = -999.0
PCMa.PCMp.skew = -999.0
PCMa.PCMp.T = -999.0
PCMa.IterateJudge = 0.05 // 収束判定条件は前ステップ見かけの比熱の5%
PCMa.Spctype = 'm' // モデル形式をデフォルトとする
PCMa.Condtype = 'm'
ct := &PCMa.Chartable
ct[0].PCMchar = 'E' // 引数0はエンタルピー
ct[1].PCMchar = 'C' // 引数1は熱伝導率
for i := 0; i < 2; i++ {
ct[i].itablerow = 0
ct[i].T = nil
ct[i].Chara = nil
ct[i].filename = ""
ct[i].tabletype = 'e'
ct[i].fp = nil
ct[i].lowA = 0.0
ct[i].lowB = 0.0
ct[i].upA = 0.0
ct[i].upB = 0.0
ct[i].minTempChng = 0.5 // 最低温度変動幅
}
(*pcm)[j] = PCMa
}
}
for i := 0; i < N; i++ {
PCMa := (*pcm)[i]
s = fi.GetToken()
if s[0] == '*' {
break
}
PCMa.Name = s // PCM名称
for {
s = fi.GetToken()
ce := strings.IndexByte(s, ';')
if ce >= 0 {
s = s[:ce]
}
if st := strings.IndexByte(s, '='); st >= 0 {
dt, _ := strconv.ParseFloat(s[st+1:], 64)
if strings.HasPrefix(s, "Ql") { // 潜熱量[J/m3]
PCMa.Ql = dt
} else if strings.HasPrefix(s, "spcheattable") {
PCMa.Chartable[0].filename = s[st+1:]
PCMa.Spctype = 't'
} else if strings.HasPrefix(s, "table") {
if s[st+1] == 'e' {
PCMa.Chartable[0].tabletype = 'e'
} else if s[st+1] == 'h' {
PCMa.Chartable[0].tabletype = 'h'
}
} else if strings.HasPrefix(s, "conducttable") {
PCMa.Chartable[1].filename = s[st+1:]
PCMa.Condtype = 't'
} else if strings.HasPrefix(s, "minTempChange") {
PCMa.Chartable[0].minTempChng = dt
PCMa.Chartable[1].minTemp = dt
} else if strings.HasPrefix(s, "Condl") { // 液相熱伝導率[W/mK]
PCMa.Condl = dt
} else if strings.HasPrefix(s, "Conds") { // 固相熱伝導率[W/mK]
PCMa.Conds = dt
} else if strings.HasPrefix(s, "Crol") { // 液相容積比熱[J/m3K]
PCMa.Crol = dt
} else if strings.HasPrefix(s, "Cros") { // 固相容積比熱[J/m3K]
PCMa.Cros = dt
} else if strings.HasPrefix(s, "Tl") { // 液体から凝固が始まる温度[℃]
PCMa.Tl = dt
} else if strings.HasPrefix(s, "Ts") { // 固体から融解が始まる温度[℃]
PCMa.Ts = dt
} else if strings.HasPrefix(s, "Tp") { // 見かけの比熱のピーク温度[℃]
PCMa.Tp = dt
} else if strings.HasPrefix(s, "DivTemp") { // 比熱数値積分時の温度分割数
PCMa.DivTemp = int(dt)
} else if strings.HasPrefix(s, "Ctype") { // 見かけの比熱の特性曲線番号
PCMa.Ctype = int(dt)
} else if strings.HasPrefix(s, "a") { // これ以降は見かけの比熱計算のためのパラメータ
PCMa.PCMp.a = dt
} else if strings.HasPrefix(s, "b=") {
PCMa.PCMp.b = dt
} else if strings.HasPrefix(s, "c") {
PCMa.PCMp.c = dt
} else if strings.HasPrefix(s, "d") {
PCMa.PCMp.d = dt
} else if strings.HasPrefix(s, "e") {
PCMa.PCMp.e = dt
} else if strings.HasPrefix(s, "f") {
PCMa.PCMp.f = dt
} else if strings.HasPrefix(s, "B") {
PCMa.PCMp.B = dt
} else if strings.HasPrefix(s, "T") {
PCMa.PCMp.T = dt
} else if strings.HasPrefix(s, "bs") {
PCMa.PCMp.bs = dt
} else if strings.HasPrefix(s, "bl") {
PCMa.PCMp.bl = dt
} else if strings.HasPrefix(s, "skew") {
PCMa.PCMp.skew = dt
} else if strings.HasPrefix(s, "omega") {
PCMa.PCMp.omega = dt
} else if strings.HasPrefix(s, "nWieght") {
PCMa.NWeight = dt
} else if strings.HasPrefix(s, "IterateJudge") {
PCMa.IterateJudge = dt
} else {
Eprint("<PCMdata>", s)
}
} else {
if s == "-iterate" {
PCMa.Iterate = true
*pcmiterate = 'y'
} else if s == "-pcmnode" {
PCMa.AveTemp = 'n'
} else {
Eprint("<PCMdata>", s)
}
}
if ce != -1 {
break
}
}
// テーブルの読み込み(見かけの比熱)
if PCMa.Spctype == 't' {
TableRead(&PCMa.Chartable[0])
}
// テーブルの読み込み(熱伝導率)
if PCMa.Condtype == 't' {
TableRead(&PCMa.Chartable[1])
}
// 入力情報のチェック
if PCMa.Spctype == 'm' {
var Tin, Bin, bsin, blin, skewin, omegain, ain, bin, cin, din, ein, fin int
var Qlin, Condlin, Condsin, Crolin, Crosin, Tlin, Tsin, Tpin int
Tin = dparaminit(PCMa.PCMp.T)
Bin = dparaminit(PCMa.PCMp.B)
bsin = dparaminit(PCMa.PCMp.bs)
blin = dparaminit(PCMa.PCMp.bl)
skewin = dparaminit(PCMa.PCMp.skew)
omegain = dparaminit(PCMa.PCMp.omega)
ain = dparaminit(PCMa.PCMp.a)
bin = dparaminit(PCMa.PCMp.b)
cin = dparaminit(PCMa.PCMp.c)
din = dparaminit(PCMa.PCMp.d)
ein = dparaminit(PCMa.PCMp.e)
fin = dparaminit(PCMa.PCMp.f)
Qlin = dparaminit(PCMa.Ql)
Condlin = dparaminit(PCMa.Condl)
Condsin = dparaminit(PCMa.Conds)
Crolin = dparaminit(PCMa.Crol)
Crosin = dparaminit(PCMa.Cros)
Tlin = dparaminit(PCMa.Tl)
Tsin = dparaminit(PCMa.Ts)
Tpin = dparaminit(PCMa.Tp)
// 必須入力項目のチェック
if Condlin+Condsin+Crolin+Crosin+Tlin+Tsin != 0 {
fmt.Printf("<PCMdata> name=%s Condl=%f Conds=%f Crol=%f Cros=%f Tl=%f Ts=%f\n",
PCMa.Name, PCMa.Condl, PCMa.Conds, PCMa.Crol, PCMa.Cros, PCMa.Tl, PCMa.Ts)
}
// モデルごとに入力値をチェックする
if PCMa.Ctype == 1 || PCMa.Ctype == 2 {
if Qlin+Tsin+Tlin != 0 {
fmt.Printf("<PCMdata> name=%s Ql=%f Ts=%f Tl=%f\n",
PCMa.Name, PCMa.Ql, PCMa.Ts, PCMa.Tl)
}
} else if PCMa.Ctype == 3 {
if Qlin+Tpin+Tin+Bin != 0 {
fmt.Printf("<PCMdata> name=%s Ql=%f Tp=%f T=%f B=%f\n",
PCMa.Name, PCMa.Ql, PCMa.Tp, PCMa.PCMp.T, PCMa.PCMp.B)
}
} else if PCMa.Ctype == 4 {
if Tpin+ain+bin != 0 {
fmt.Printf("<PCMdata> name=%s Tp=%f a=%f b=%f\n",
PCMa.Name, PCMa.Tp, PCMa.PCMp.a, PCMa.PCMp.b)
}
} else if PCMa.Ctype == 5 {
if Tpin+bsin+blin+ain != 0 {
fmt.Printf("<PCMdata> name=%s Tp=%f bs=%f bl=%f a=%f\n",
PCMa.Name, PCMa.Tp, PCMa.PCMp.bs, PCMa.PCMp.bl, PCMa.PCMp.a)
}
} else if PCMa.Ctype == 6 {
if Qlin+Tpin+skewin+omegain != 0 {
fmt.Printf("<PCMdata> name=%s Ql=%f Tp=%f skew=%f omega=%f\n",
PCMa.Name, PCMa.Ql, PCMa.Tp, PCMa.PCMp.skew, PCMa.PCMp.omega)
}
} else if PCMa.Ctype == 7 {
if ain+bin+cin+din+ein+fin != 0 {
fmt.Printf("<PCMdata> name=%s a=%f b=%f c=%f d=%f e=%f f=%f\n",
PCMa.Name, PCMa.PCMp.a, PCMa.PCMp.b, PCMa.PCMp.c, PCMa.PCMp.d, PCMa.PCMp.e, PCMa.PCMp.f)
}
}
}
}
}
// PCMの物性値テーブルの読み込み
func TableRead(ct *CHARTABLE) {
if ct.filename == "" {
return
}
fp, err := os.Open(ct.filename)
if err != nil {
fmt.Printf("<PCMdata> xxxx file not found %s xxxx\n", ct.filename)
return
}
ct.fp = fp
// 設定されている行数を数える
var c byte
var row int
row = 0
// Count the number of rows
for {
_, err := fmt.Fscanf(ct.fp, "%c", &c)
if err != nil {
break
}
if c == '\n' {
row++
}
}
// Close the file temporarily
err = ct.fp.Close()
if err != nil {
fmt.Println("<PCMdata> ファイルのクローズに失敗")
return
}
// Allocate memory
ct.T = make([]float64, row)
ct.Chara = make([]float64, row)
if ct.T == nil || ct.Chara == nil {
fmt.Println("<PCMdata> メモリの確保に失敗")
}
// Reopen the file
fp, err = os.Open(ct.filename)
if err != nil {
fmt.Println("<PCMdata> ファイルのオープンに失敗")
return
}
ct.fp = fp
var st, tt int
var spheat, prevheat, prevTemp float64
prevheat = 0.0
prevTemp = -999.0
spheat = 0.0
for i := 0; i < row; i++ {
T := &ct.T[i]
Char := &ct.Chara[i]
// 温度の読み込み
_, err := fmt.Fscanf(ct.fp, " %f ", T)
if err != nil {
fmt.Println("<PCMdata> 温度の読み込みに失敗")
break
}
// テーブルの下限温度
if st == 0 {
ct.minTemp = *T
prevTemp = *T
st = 1
}
// 昇順に並んでいないときのエラー表示
if *T <= prevTemp && tt == 1 {
fmt.Printf("i=%d 温度データが昇順に並んでいません T=%f preT=%f\n", i, *T, prevTemp)
}
var dblTemp float64
// 特性値の読み込み
_, err = fmt.Fscanf(ct.fp, " %f ", &dblTemp)
if err != nil {
fmt.Println("<PCMdata> 特性値の読み込みに失敗")
break
}
*Char = dblTemp
// 見かけの比熱の場合
if ct.tabletype == 'h' {
*Char = prevheat + spheat*(*T-prevTemp)
}
prevheat = *Char
prevTemp = *T
spheat = dblTemp
// テーブルの上限温度
ct.maxTemp = *T
ct.itablerow++
tt = 1
}
// Close the file
err = ct.fp.Close()
if err != nil {
fmt.Println("<PCMdata> ファイルのクローズに失敗")
}
// 上下限温度範囲以外の線形回帰式を作成
// 下限以下
ct.lowA = (ct.Chara[0] - ct.Chara[1]) / (ct.T[0] - ct.T[1])
ct.lowB = ct.Chara[0] - ct.lowA*ct.T[0]
// 上限以上
ct.upA = (ct.Chara[ct.itablerow-1] - ct.Chara[ct.itablerow-2]) / (ct.T[ct.itablerow-1] - ct.T[ct.itablerow-2])
ct.upB = ct.Chara[ct.itablerow-1] - ct.upA*ct.T[ct.itablerow-1]
}
// 初期化されているかをチェックする
func dparaminit(A float64) int {
if math.Abs(A-(-999.0)) < 1e-5 {
return 1
} else {
return 0
}
}
func PCMcount(fi *EeTokens) int {
N := 0
add := fi.GetPos()
for fi.IsEnd() == false {
s := fi.GetToken()
if strings.HasPrefix(s, "*") {
break
}
if strings.HasPrefix(s, ";") {
N++
}
}
fi.RestorePos(add)
return N
}
// 固相、液相物性値と潜熱量からPCM温度の物性値を計算する(比熱、熱伝導率共通)
// 熱伝導率の計算時はQl=0とする
func FNPCMState(Ctype int, Ss, Sl, Ql, Ts, Tl, Tp, T float64, PCMp *PCMPARAM) float64 {
var Qse, Qla, Tls float64
Tls = Tl - Ts
// 顕熱分の補間
if T > Ts && T < Tl {
Qse = Ss + (Sl-Ss)/Tls*(T-Ts)
} else if T <= Ts {
Qse = Ss
} else {
Qse = Sl
}
// 潜熱分の補間
Qla = 0.0
Ql = 0.0
// 熱伝導率計算の場合は潜熱ゼロ
if Ctype == 0 {
Qla = 0.0
} else if Ctype == 1 { // 潜熱変化域潜熱比熱一定
if T > Ts && T < Tl {
Qla = Ql / (Tl - Ts)
}
} else if Ctype == 2 { // 二等辺三角形
if T > Ts && T < Tl {
if T < (Tl+Ts)/2.0 {
Qla = 4.0 * Ql / (Tls * Tls) * (T - Ts)
} else {
Qla = 4.0 * Ql / (Tls * Tls) * (Tl - T)
}
}
} else if Ctype == 3 { // 双曲線関数
Temp := math.Cosh((2.0 * PCMp.B / PCMp.T) * (T - Tp))
Qla = 0.5 * Ql * (2.0 * PCMp.B / PCMp.T) / (Temp * Temp)
} else if Ctype == 4 { // ガウス関数(対象)
Temp := (T - Tp) / PCMp.B
Qla = PCMp.a * math.Exp(-0.5*Temp*Temp)
} else if Ctype == 5 { // ガウス関数(非対称)
Temp := (T - Tp) / PCMp.B
if T <= Tp {
Temp = PCMp.bs
} else {
Temp = PCMp.bl
}
Qla = PCMp.a * math.Exp(-((T-Tp)/Temp)*((T-Tp)/Temp))
} else if Ctype == 6 { // 誤差関数歪度
//Temp := math.Exp(-(T - Tp) * (T - Tp) / ((2.0 * PCMp.omega) * PCMp.omega))
Qla = Ql / math.Sqrt(2.0*math.Pi) * math.Exp(-(T-Tp)*(T-Tp)/((2.0*PCMp.omega)*PCMp.omega)) * (1.0 + math.Erf(PCMp.skew*(T-Tp)/(math.Sqrt(2.0)*PCMp.omega)))
} else if Ctype == 7 { // 有理関数
if T < 0 {
Qla = 0.0
} else {
Qla = math.Pow(T, PCMp.f) * (PCMp.a*T*T + PCMp.B*T + PCMp.c) / (PCMp.d*T*T + PCMp.e*T + 1.0)
}
}
return (Qse + Qla)
}
/* ------------------------------------------------------ */
// PCMの状態値計算(家具用)
func FNPCMStatefun(Ctype int, Ss, Sl, Ql, Ts, Tl, Tp, oldT, T float64, DivTemp int, PCMp *PCMPARAM) float64 {
var dblTemp, dTemp, TPCM, Qld, Qllst float64
dTemp = (T - oldT) / float64(DivTemp)
Qllst = 0.0
// 前時刻からの温度変化が小さい場合は時間積分はしない
if math.Abs(T-oldT) < 1e-4 {
Qllst = FNPCMState(Ctype, Ss, Sl, Ql, Ts, Tl, Tp, (T+oldT)*0.5, PCMp)
} else {
// 見かけの比熱の時間積分
for i := 0; i < DivTemp+1; i++ {
TPCM = oldT + dTemp*float64(i)
Qld = FNPCMState(Ctype, Ss, Sl, Ql, Ts, Tl, Tp, TPCM, PCMp)
dblTemp += Qld
}
Qllst = dblTemp / float64(DivTemp+1)
}
return Qllst
}
// PCMの温度から見かけの比熱を求める(テーブル形式)
// Told:前時刻のPCM温度、T:暫定現在時刻PCM温度
func FNPCMstate_table(ct *CHARTABLE, Told, T float64, Ndiv int) float64 {
var oldEn, En, Chara float64
// 前時刻と暫定現在時刻のPCM温度が同温の場合
if math.Abs(T-Told) < ct.minTempChng {
Tave := 0.5 * (T + Told)
// 見かけの比熱の計算
if ct.PCMchar == 'E' {
oldEn = FNPCMenthalpy_table_lib(ct, Tave-0.5*ct.minTempChng)
En = FNPCMenthalpy_table_lib(ct, Tave+0.5*ct.minTempChng)
// 見かけの比熱の計算
Chara = (En - oldEn) / ct.minTempChng
} else {
// 熱伝導率の計算
Chara = FNPCMenthalpy_table_lib(ct, Tave)
}
} else {
if ct.PCMchar == 'E' {
oldEn = FNPCMenthalpy_table_lib(ct, Told)
En = FNPCMenthalpy_table_lib(ct, T)
// 見かけの比熱の計算
Chara = (En - oldEn) / (T - Told)
} else {
// 熱伝導率の計算
var dblTemp, Tpcm, dTemp float64
dblTemp = 0.0
dTemp = (T - Told) / float64(Ndiv)
// ToldからTまでを積分する
for i := 0; i < Ndiv; i++ {
Tpcm = Told + dTemp*float64(i)
dblTemp += FNPCMenthalpy_table_lib(ct, Tpcm)
}
Chara = dblTemp / float64(Ndiv+1)
}
}
return Chara
}
func FNPCMenthalpy_table_lib(ct *CHARTABLE, T float64) float64 {
var prevTpcm, Tpcm, enthalpy, preventhalpy *float64
var retVal float64
if T < ct.minTemp {
// テーブルの最低温度よりPCM温度が低い場合は端部の特性値を線形で外挿
retVal = ct.lowA*T + ct.lowB
} else if T > ct.maxTemp {
//テーブルの最高温度よりPCM温度が高い場合は端部の特性値を線形で外挿
retVal = ct.upA*T + ct.upB
} else {
for i := 0; i < ct.itablerow; i++ {
Tpcm := &ct.T[i]
enthalpy := &ct.Chara[i]
if *Tpcm > T {
break
}
prevTpcm = Tpcm
preventhalpy = enthalpy
}
// 線形補間
retVal = *preventhalpy + (*enthalpy-*preventhalpy)*(T-*prevTpcm)/(*Tpcm-*prevTpcm)
}
return retVal
}