電気通信大学大学院電気通信学研究科博士(工学)の学位申請論文

(1)

多目的大空間におけるファジィモデルを用いた空調制御に関する研究

植田俊克

電気通信大学大学院電気通信学研究科博士(工学)の学位申請論文

2009 年 3 月

(2)

多目的大空間におけるファジィモデルを用いた空調制御に関する研究

博士論文審査委員会

主査　本多　中二　　教授

委員　萩野　剛二郎　教授

委員　新　　誠一　　教授

委員　板倉　直明　　准教授

委員　酒井　孝司　　准教授

(3)

著作権所有者植田　俊克

2009

(4)

A Study on HVAC Control Based on a Fuzzy Model for Multipurpose Large Indoor Spaces

Toshikatsu Ueda

Abstract

In large indoor spaces, air volume is large, and there is great distance between the occu- pied zones and air outlets as well. As a result, it is easy for time delays and instabilities to arise with HVAC (Heating, Ventilating and Air-conditioning) systems for large spaces. In addition, large spaces such as dome type baseball stadiums and event halls are used for a variety of purposes, and there are limitations on the locations where the sensors for an HVAC control system can be installed. As a result, it is difficult to obtain representative temperature and air velocity distributions of occupied zones. This paper proposes an HVAC control system using a fuzzy model as a method to resolve these problems. A fuzzy model is constructed by using data obtained from estimating various parameters such as outlet air temperature, outlet air volume, and human body load, which affect the temperature distribution within a space, by using CFD (Computational Fluid Dynamics) simulations. Fuzzy models adopt a simplified fuzzy reasoning method, and can estimate nonlinear temperature characteristics with high accuracy. Rules for fuzzy models are automatically generated using data from CFD simulations. The HVAC control presented here estimates the averaged temperature of an occupied zone based on the tempera- ture of the sensor location using a fuzzy model, and conducts PID control so that the estimated temperature becomes the set temperature. From the results of control simulation using CFD, this HVAC control based on a fuzzy model is able to achieve control with little offset as compared with conventional PID control where the temperature at the sensor location is the representative temperature of the occupied zone.

This paper is comprised of six chapters, and a short summary of each chapter is given below.

Chapter 1 consists of the background behind this study and an abstract of the method proposed in this study.

Chapter 2 describes fuzzy modeling and CFD, which are elemental technologies in this study, and also explains the control method (FM-PID) based on the fuzzy model.

(5)

Chapter 3 presents a case study using CFD simulation based on 9 cases—consisting of 3 different usage purpose and 3 sensor positions—as parameters. In this chapter, a comparison is made between FM-PID control in a CAV (constant air volume) system and PID control, where the sensor temperature serves as the representative temperature. Based on the results of this comparison, it was demonstrated that with FM-PID control, it is possible for the temperature in an occupied zone to serve as the set temperature without depending on the usage purpose or sensor location.

Chapter 4 presents a case study similar to the one in Chapter 3, and also reviews control performance in a VAV (variable air volume) system. A VAV system is superior in energy conservation as compared to the CAV system, and is used widely in offices, etc. However, since temperature distribution changes greatly when the VAV system is used in large, indoor spaces, there are very few examples of the VAV system being used in such spaces. Based on the results of control simulation using CFD, it was demonstrated that FM-PID control can be suitable for VAV systems in large, indoor spaces.

Chapter 5 introduces a control simulation with the VAV system that takes the introduction of outside air into consideration, which demonstrates the high control performance and energy conservation effects of FM-PID control.

Chapter 6 consists of the results of this study as well as future issues for research.

(6)

多目的大空間におけるファジィモデルを用いた空調制御に関する研究

植田　俊克

論文概要

　体育館，ドーム球場，コンサートホール，イベントホール，アトリウムなどに代表される大空間は気積が大きいという共通の特徴を有し，さらにドーム球場，

イベントホールは様々な用途に利用される．

　大空間空調ではその気積の大きさから吹出し口と居住域が離れている場合が多く，時間遅れや制御の不安定が生じやすい．さらに多目的に用いられる大空間では使用用途によって外乱である人員密度分布が異なることから，適切なセンサー位置が未知であることが多く，温度分布の不均一や居住域の温度と設定値にずれが生じやすい．そのため大空間での空調制御を円滑に運用するためには熟練のオペレータによる長い経験と勘に因るところが多い．また，京都議定書が発効され，膨大なエネルギーを消費する大空調には省エネルギーへの取り組みが必須となっている．

　そこで，本研究では多目的に使用されるドーム球場において，センサー位置や使用用途の影響を受けずに快適かつ省エネルギー性の優れた空調制御法を提案し，数値流体力学( 以下，C F D とする) 解析によりその有効性を検証する．提案する制御法は，外乱に応じた空間の温度分布をあらかじめ求めておき，温度分布を予測しながら制御を行うもので，温度分布の予測にはファジィモデルを，制御に関しては P I D 制御を用いる．ファジィモデルは非線形性の強いデータに対して高精度に表現でき，得られたメンバーシップ関数からデータ構造を容易に理解できるという特徴を有している．また，P I D 制御に関しては既設の大空間で多く採用されているため，今までの膨大な知見をそのまま用いることができ，さらに本制御法の既設建物への導入が容易になるという利点がある．本論文は全６章から構成されており，各章の概要を以下に示す．

(7)

　第１章『序論』では，本論文の研究背景および既往の研究と本研究の目的，

本制御法の概要を述べる．

　第２章『ファジィモデルを用いた P I D 制御』では，本研究の要素技術であるファジィモデリングと C F D 解析に関して述べた後，本制御法を提案する．ファジィモデリングは簡略化ファジィ推論法を用いて，データから自動的にルールを作成することができる方法としている．また，C F D はデータを作成する際と制御シミュレータとして用いており，実現象を詳細に予測する．提案する制御法は操作量( 吹出し温度，吹出し風量) と検出量( 温度センサー値) ，検出器( 温度センサー) 位置座標を入力値としたファジィモデルを用いて居住域の平均温度を予測し，居住域の平均温度が目標値となるように P I D 制御を行うもので，

任意の検出器の位置で任意の居住域の空調制御が可能となっている．

　第３章『定風量方式』では，大空間空調で最も一般的に用いられている風量を一定として吹出し温度を制御する定風量方式において，使用用途３パターン，センサー位置３パターンの９ケースに対して C F D 解析による制御シミュレーションを行い，センサー温度を代表温度とする P I D 制御と本制御法を比較することで，本制御法が制御性を損なうことなく，いかなる使用用途，センサー位置においても，居住域温度をほぼ目標値とすることができることを示す．さらに，さまざまな外乱に対して柔軟に対応可能であることを示す．　第４章『変風量方式』では，居住域温度と吹出し温度との差を可能な限り大きくしつつ吹出し風量を制御する変風量方式において第３章と同様の解析を行う。変風量方式は事務室で多く採用され，定風量方式よりエネルギー消費が少ないという特徴を有するが，空間の温度分布が不均一になりやすく大空間ではほとんど採用されていない。そこで，任意の位置での空間温度を予測できる本制御法が，変風量方式の大空間空調に有効であることを示す．

　第５章『外気導入を考慮した変風量方式』では，本制御法を用いた変風量方式の実建物への導入を想定し，人体が発生源となる二酸化炭素濃度が許容値以下となる必要外気量を確保した上での制御性を検証し，安定した制御が可能であることを示した後，省エネルギー効果を試算する．

　第６章『結論』では，本研究の成果と今後の課題および展望をまとめる．

(8)

第１章　序論

1 . 1 　背景

　ドーム球場，イベントホール，体育館，音楽ホール，エントランス，アトリウムなどの大空間は，気積が大きいという共通の特徴を有する．さらに，ドーム球場，イベントホール，体育館では，競技および競技観戦，イベント，コンサート，

集会など多目的に利用される．

　大空間は近年まで屋外で行われていたイベントを全天候型の室内へ移動した例が多く，音楽ホールを除いた初期の大空間は，屋外よりも快適であれば良いといった低いグレードの空調で十分であった．しかし，次々と新しい大空間が誕生し，より質の高い空調に対する要求が年々高まっている[ 1 ] ．

　大空間の空調を行う際，気積が大きく，吹出し口と居住域の距離が離れていることから，居住域に温度ムラが生じやすく制御時には時間遅れが問題となる．さらに多目的に利用される大空間の空調では，用途により人員密度分布が異なるため，用途ごとの気流・温度分布を考慮した空調を行う必要がある．しかし，実際には主用途のみで空調計画される場合も多く，その他の用途に関しては運用開始後の室内温度データを用いて設定温度，制御パラメータなどを調整しているのが現状である．調整期間は，四季による変化を考慮するため，最低でも 1 年，場合によっては各用途の頻度や試行錯誤により数年間が必要となる．このような状況の中で少子高齢化により熟練した設備管理者が減少しており，経験的に設定値を変更することが困難となりつつある．また，大空間のほとんどが無柱空間であるため温度センサーの設置位置が限定されてしまうことが空調制御を一層困難なものとしている．

　一方，地球環境問題に果たす空調の役割は大きく，建物の冷暖房が民生業務部門の二酸化炭素排出量の 1 / 3 を占める[ 2 ] ことからも省エネルギー対策は必須事項である．その中でも大空間は，全体の割合は少ないものの，気積の大きな空間の空調に膨大なエネルギーを必要とする．

　以上から，これからの大空間の空調には，快適性と省エネルギー性の両立が求められる．また建物の構造的な耐久年数に対し，空調設備の寿命は短く，今まさ

(12)

にバブル期に建築された事務所ビル，大空間の設備更新が始まったところである．

よって新築建物だけでなく既設建物の空調設備の改修が容易であることもこれからの空調設備に必要な要素となる．

1 . 2 　既往の研究と本研究で提案する制御法

　1 0 0 m を超える規模の大空間建造物の空調システムに関する実規模実験を行うことはその規模の大きさからほぼ不可能であり，古くから模型実験が行われてきた．

1 9 6 4 年に代々木国立屋内競技場で勝田ら[ 3 ] がはじめて大空間の空調設計に模型実験を用いて以来，勝田・土屋による一連の研究[ 4 - 6 ] を経て強制対流場での近似相似則が整備され，縮尺相似模型と A r ( A r c h i m e d e s ) 数の一致を前提とした実験が最も実績のある検証方法となっている．また，近年ではコンピュータの高速化と数値解法の高速化・高精度化により N a v i e r - S t o k e s の運動方程式を基礎方程式とした流れの数値解析( 以下，C F D 解析とする）が，模型実験に代わる新たな予測手法として注目されており，村上ら[ 7 ] は模型実験と C F D 解析との比較を，さらに尾関ら[ 8 ] は実大実験と C F D 解析の比較を行うことで精度の確認を行っている．実規模の実験が困難で模型実験に頼らざるを得ない大空間空調の実現象予測に，実スケールでの流れの挙動を再現できる C F D 解析は有効な手段となっており，現段階において C F D 解析は模型実験と同程度の精度を有するまでになっている．

　一方，既往の大空間空調制御の研究において，岩本ら[ 9 , 1 0 ] は模型実験により適切な温度センサー配置に関する検討を行い，さらに橋本[ 1 1 ] が行った P I D 制御パラメータのチューニングに用いたニューラルネットワーク手法を予測制御へ適用し，時間遅れの低減を図っている．さらに，樋渡・庄司・大野ら[ 1 2 - 1 4 ] は適切な温度センサー位置決定のため，換気効率指標を用いて竣工後の実測と C F D 解析を行い，C F D 解析の精度を確認した上で加藤らの温熱環境形成寄与率[ 1 5 ] ( 以下，C R I ) を C F D 解析で求め，温度センサー配置の適正化を図っている．しかし，

多目的使用時にセンサー温度と居住域温度が異なる場合の研究は行われていない．

　そこで本研究では，あらかじめ熱負荷に応じて空間温度分布が推定可能なモデルを用いて，任意の居住域を制御対象とすることのできる制御法を提案する．任

(13)

意の位置での温度を推定する方法としては，C F D 解析を用いる方法があるが，C F D 解析は温度推定に膨大な計算時間を要し，リアルタイム性が求められる空調制御に導入することは不可能である．そこで，あらかじめ実験・C F D 解析等により得られた温度分布データを用いて予測モデルを作成することとする．モデルの作成に関しては，C R I を用いることも可能であるが，C R I は気流分布を固定しているため，冬期の大空間空調で生じる人体負荷の増加に伴う暖房から冷房へのダイナミックな気流の変化への対応には適さない．したがって，本研究では構造が容易で，メンバーシップ関数により温度分布性状が大局的に把握できるファジィモデルを用いることとする．ファジィモデルは冷暖房空調時の室内温度分布データのように非線形性の強いデータに対して高精度に表現できるという特徴を有している[16]．

　ファジィモデルの導出（以下，ファジィモデリング）に関しては簡略化ファジィ推論法を用いてルールを自動的に調整する荒木らの方法[ 1 7 ] を基本としてデータの特性に合わせて改良を加える．すなわち最急降下法を用いて推論誤差がもっとも小さくなるようにメンバーシップ関数を決定する．その後，誤差の最も大きいデータ部に新たにルールを作成して推論誤差をさらに小さくする．この操作を繰り返すことにより，モデルを自動的に高精度で構築することが可能となる．

　制御に関しては非線形性の強い制御対象に対し，ファジィ推論に基づくファジィ制御が有効であることが示されている[ 1 8 - 2 3 ] ．ファジィ制御は制御則をファジィ推論形式で記述するが，ファジィ推論のルール作成に試行錯誤を必要とすることと既設建物への導入のしやすさを考慮し，本研究では実績の豊富な P I D 制御を用いることとする．

　以上より本研究で提案する制御法は，測定されるデータからファジィモデルを用いて用途毎に異なる居住域の平均温度を推定し，推定値と目標値の偏差から P I D 制御を行う方法となる．

(14)

第２章　ファジィモデルを用いた PID 制御

2 . １　ファジィモデリング

2 .1 .1 　簡略化ファジィ推論法

　人間が日常行っている様なあいまいな推論をファジィ集合を用いて模擬することをファジィ推論と呼び，一般に以下のように定式化できる．

ルールi： IF x₁ is A_1i and ･･･ and x_n is A_ni

THEN y_i=B_i　(i=1,2,･･･,r) ･･････(1)

ここで，IF 以降 TH EN までの部分を前件部，TH EN 以降の部分を後件部といい，x は前件部変数( 入力変数) ，y は後件部変数( 出力変数) ，A_nはx_nでのファジィ集合，B は y での結論部，iはルール番号，rはルール数，nは前件部変数の数である．ファジィ推論には結論部B_i をファジィ集合で与える方法[ 2 4 ] ，入力変数の関数とする方法 [ 2 5 ] ，実数c_iで与える方法[ 2 6 ] があり，実数で与える方法を簡略化ファジィ推論法という．簡略化ファジィ推論法は他のファジィ推論法に比べ , モデリングの自由度が制限されるものの実用上十分な精度を有しており，近年利用が盛んである [ 2 0 ] ．さらに結論部が実数であるため，ルールの自動生成に適している[ 2 6 , 2 7 ] ．　簡略化ファジィ推論プロセスは，以下のように示される[ 2 8 ] ．

① 適合度h_iを求める．

h

_i

= Π

ⁿ_k=1

μ A

_ni

( x

_n⁰

)

･･････(2)

　ここで，μ A_{n i}(x_n⁰) はA_{n i}のx_n⁰におけるメンバーシップ値である．

② 加重平均を求める．

Σ

^r_i=1

h

_i

c

_i

y

₀

= 　　　･･････(3) Σ

^r_i=1

h

_i

　ここで，y₀ は出力値である．

　例としてA のメンバーシップ関数を三角型，r= 4 ，n= 2 ，入力値をx_{1 0}，x_{2 0}，出力値をy₀ とした場合の推論プロセスを図 2 - 1 に示す．ファジィモデルは入力空間をファジィ分割し，各空間毎に入出力関係を構築するため，複雑であったり数式で近似することが困難な入出力関係を持つシステムに有効な方法である．さらに，

(15)

後件部実数値とメンバーシップ関数の関係 1 c₁

c₂

c₃ c₄ y₀

x₁₀ x₂₀

A₁₁ A₁₂

a₁₁

a₁₂

0 A21A22 a21

a22 01

y₀= (h₁c₁+h₂c₂+h₃c₃+h₄c₄)/(h₁+h₂+h₃+h₄)

① 適合度(h_i)

② 出力値(y₀)

A₁₁

0 1

x₁₀ a₁₁

A21

0 1

x₂₀ a21

c₁

A₁₁

0 1

x₁₀

a₁₁ A₂₂

0 1

x₂₀ a₂₂

c₂

A12

0 1

x₁₀ a₁₂

A₂₁

0 1

x₂₀ a₂₁

c₃

A₁₂

0 1

x10

a₁₂ A₂₂

0 1

x20

a₂₂

c₄

h₁=a₁₁a₂₁

h₂=a₁₁a₂₂

h₃=a₁₂a₂₁

h₄=a₁₂a₂₂

図 2 - 1 　簡略化ファジィ推論プロセス( ２入力１出力系)

(16)

多変数への拡張が容易で，メンバーシップ関数の分割割合で出力変数の性質が直感的に類推できるという特徴を有している．

2 .1. 2　ルール生成法

　ファジィモデルは，ルールの記述に用いるメンバーシップ関数，分割方法によって精度が大きく異なるため，ルールの作成には逐次的にルールを生成する荒木らの方法[ 1 7 ] を基本とし，データの特性に合わせて改良を加えることで精度向上を図る．図 2 - 2 に１入力１出力系におけるルールの生成手順を示す．図中R はメンバーシップ関数で分割された領域である．荒木らの方法は簡略化ファジィ推論法を用いて，そのパラメータであるルール後件部の実数値を調整しながら逐次的に新たなルールを自動的に獲得する方法で，パラメータの調整は最急降下法，ルール

(b) 本法

図 2 - 2 　後件部実数値の生成例( １入力１出力系) (a ) 荒木らの方法

Y

X

R₁

Y

X

R₁ R₂

Y

X

R₁ R₂ R₃ A₁ A₂ A₁ A₂A₃ A₄ A₁ A₂A₃ A₄A₅ A₆

Y

X

R₁

Y

X

R₁ R₂ Y

X

R₁ R₂ R₃

A₁ A₂ A₁ A₂A₃A₄ A₁ A₂A₃～ A₆

(17)

の生成は図 2 - 2 ( a ) に示すように，前件部領域の中で最も推論誤差の大きい領域を等分割するように１つずつ生成していく．前件部変数が格子状に配置されていない場合に最急降下法を用いて誤差が最も小さくなる後件部実数値を決定するため，精度の高いモデリングが可能となる．一方本提案法は，図 2 - 2 ( b ) に示すように，ルール生成の際，前件部領域の中でもっとも推論誤差の大きい値の位置でメンバーシップ関数を分割する点が特徴となっている．これにより，非線形性の強いデータに対して効率よくルールを作成できる．反面，ノイズの多いデータに対しては推論誤差が大きくなるという欠点を有するが，C F D 解析結果や関数近似のようにノイズの発生しにくいデータに関しては有効である．

　図 2 - 3 にモデリングのフローチャートを示す．メンバーシップ関数は三角型であり、まず各入力データの最大最小値を端点とした初期メンバーシップ関数を作成する．次にモデルの誤差が最大となる部分でメンバーシップ関数を分割し，収束判定条件を満たすまでこれを繰り返すことで，より精度の高いモデルを生成していく．以下に手順を示す．

①入出力データの読み込み：多入力１出力データを読み込む．その際，トレーニングデータとチェッキングデータには分割せず，すべてのデータをモデリングに使用する．

②初期ファジィモデルの作成：各入力データの最大・最小値を端点とした初期メンバーシップ関数を作成する．後件部実数値が出力データに無い場合は最急降下法により求める．

③ファジィモデルの最大誤差算出：ファジィモデルから得られるすべての入力データ位置での推定値と出力データを比較し，最大誤差を求める．

④収束判定：収束判定条件を満足すればモデリングを終了する．収束判定条件は，最大誤差による場合，最大誤差が前回よりも大きくなった場合，

ルール数による場合のいずれかを選択することとする．

⑤最大誤差位置でメンバーシップ関数分割：収束判定条件を満足しない場合は，最大誤差位置でメンバーシップ関数を分割し，出力データを後件部実数値に代入する．その際，得られたメンバーシップ関数の中で後件部実数値が出力データに無い場合は，前回のファジィモデルの推定値を代入する．

(18)

⑥最急降下法にて後件部実数値決定：出力データを代入していないすべての後件部実数値に最急降下法を適用する．なお，分割範囲のデータが１個以下と少ない場合は，過度なフィッティングを避けるため最急降下法を使わない．

⑦繰り返し：③へ戻る．

図 2 - 3 　ファジィモデリングフロー

注^† 初期ファジィモデルは各入力データの最小値・

最大値を端点とし、出力データが存在しない場合は最急降下法にて後件部実数値を決定

最大誤差位置でメンバーシップ関数分割入出力データ読込

初期ファジィモデル^†の作成ファジィモデルの最大誤差算出

収束 no

yes END

最急降下法にて後件部実数値決定 START

(19)

2 . 2 　室内気流・温度場の予測

2.2.1　CFD 　

　CFD は，Conputational Fluid Dynamics の頭文字をとっており，日本語では数値流体力学，あるいは計算流体力学と訳されているが，現在では”C F D ”とそのまま使う場合が多くなっている．C F D を駆使してコンピュータで流れ場を計算することを C F D 解析と呼んでいる．

　C F D では，流れの支配方程式を数値的に解いて，流れの現象を解明することになるが，膨大な計算量を必要とし，2 0 年ほど前まではスーパーコンピュータを用いて解析が行われていた．しかし近年のコンピュータ性能の飛躍的な向上により，

現在ではパーソナルコンピュータを用いて室内の流れ場を解くことが可能になっている．

　C F D の研究は，理学・工学・医学など幅広い分野で行われており，理論的な研究の他，様々な流れ現象を C F D で解き，実験と比較することで精度向上が図られている．応用研究としては，気象学，航空宇宙学，機械工学，土木・建築・環境工学，医療工学などが挙げられる．また，実現象の予測ツールとして，設計開発，

性能・環境改善などに実用化され，今や実験に代わる手法として数値実験と呼ばれるまでになっている．さらに，航空宇宙分野や原子力関連など，実験に膨大な時間と費用を必要としたり，実験自体が不可能な場合の現象解明手段として不可欠な存在となっている[ 2 9 ] ．

　本研究で扱っている建築設備分野においては，ファン，ポンプ，バルブ，吹出し・吸込み口などの流体を扱う製品開発のみならず，空調空間の流れ・温度場予測，通風・換気などの非空調空間の換気性能把握，ボイラ排気の近隣への影響，

エアコン室外機周辺の温度予測などに利用され，扱う物理量も水流・気流・温度だけではなく，湿度，粉塵濃度，揮発性化学物質濃度など多岐に渡っている[ 3 0 - 4 0 ] ．なお，本研究では C F D 解析に商用コード『S T R E A M 』（( 株) ソフトウエアクレイドル）を用いる．

(20)

2.2.2　CFD の数値解法

　流れを数値的に解く際，まず扱う流体が圧縮性流体であるか非圧縮性流体であるかによって扱う支配方程式が異なる．圧縮性流体( 気体) の場合は気体の状態方程式と N a v i e r - S t o k e s の運動方程式( 以下，N - S 方程式とする) が，非圧縮性流体の場合は連続の式と N - S 方程式が支配方程式となる．さらに温度を解析する場合はエネルギー方程式が，濃度に関しては濃度方程式が必要となる．空調空間を C F D 解析する場合，扱う流体は空気( 気体) であり，本来は圧縮性を示すが，空調空間の流れは非圧縮性を示すことが知られている[ 4 1 ] ．建築設備分野で圧縮性流体として C F D 解析を行うのは，温度差の著しい流れ，特殊ガスが床を這うような２流体間の気体の重さを考慮する流れ，著しい圧力変動を生じる流れのような特殊用途に限られる．非圧縮性流体の温度を含めた支配方程式を表 2 - 1 に示す．

　C F D 解析では，流体を微少な格子に分割し，その格子（あるいは交点）１つ１つに支配方程式を適用する．格子系は構造格子系と非構造格子系に，構造格子系はさらにデカルト( 直交) 座標系，円筒座標系，一般座標系に分けられ，非構造格子は格子の形が，4 面体，6 面体，ピラミッド，三角柱などに分けられる[ 4 2 ] が，

空調空間では構造格子デカルト座標系が，扱いやすさ，計算速度，計算精度の点から多く用いられており，本 C F D 解析でも構造格子デカルト座標系を用いる．

　離散化は，代表的には有限差分法[ 4 3 ] ，有限体積法[ 4 4 ] ，有限要素法[ 4 5 ] が挙げられるが，こちらも扱いやすさ，計算速度，計算精度，計算メモリの点から有限体積法とする．また，有限体積法に用いる運動方程式の移流成分項( 表 2 - 1 中，

( 4 ) 式左辺第２項) の差分スキームは，計算安定化のために，流れの上流側に重み付けを行う必要があるが，その近似精度が流れ場全体の精度に大きく影響するため，本解析では 3 次精度を有する Q U I C K スキーム[ 4 6 ] とする．

　次に，連続の式と N - S 方程式を連立して解く必要が生じるが，その方法を計算アルゴリズムと呼び，M A C [ 4 7 ] 系と S I M P L E [ 4 8 ] 系に大別される．本 C F D では計算の安定性，計算速度の観点から S I M P L E 系の S I M P L E C 法[ 4 9 ] を用いる．

(21)

質量保存の式 x 0

u

j

j =

¶

運動量保存の式

(

₀

)

_i₃

i j j

i j

j j

j

i i g T T

x u x u x

x p 1 x

u u t

u n b d

r - -

ïþ ïý ü ïî

ïí ì

÷÷ ø ö çç

è æ

¶ + ¶

¶

¶ + ¶

¶ - ¶

¶ = + ¶

¶

エネルギー保存の式

÷÷ ø ö çç

è æ

¶

= ¶

¶ + ¶

¶

j r j j

j

x T P x x

Tu t

T n

拡散物質保存の式

÷÷ ø ö çç

è æ

¶

= ¶

¶ + ¶

¶

j m j j

j

x D C x x

Cu t

C

　　　　　　　　　　　　　　　　　　　　・・・(4) 　

・・・(5) 　　　　　　　　　　　　

・・・(6)

・・・(7) 表 2 - 1 　支配方程式

ここに、u_i^：x_i方向の速度[m/s], x_i：位置座標[m],

t：時間[s],

p：流体の圧力[Pa], r：流体の密度[kg/m³],

n

：動粘性係数[m²/s],

g：重力加速度[m/s²], b：体膨張率[1/℃],

d

ij：Kroneckerのデルタ,

T：流体の温度[℃], T₀：流体の基準温度[℃], Pr：Prandtl数[-],

C：拡散物質の濃度[m³], D_m：拡散係数[m²/s], 下添え字はEinsteinの総和規約に従う

(22)

2.2 .3　乱流

　空調分野で扱う流れは，ごく一部を除き乱流である．乱流は大小様々な渦で構成され，互いに干渉して分裂し，最終的には熱となる．C F D を用いて N - S 方程式を厳密に解くためには，最小の渦を再現できるスケールで格子分割する必要があり，空調空間に適用するためには現在のスーパーコンピュータを用いても不可能である[ 5 0 ] ．そのため，コンピュータを用いて実用的な時間で解を得るために，

乱流のモデル化を行う．現在，主流となっている方法は k - ε型 2 方程式モデル [51](以下，k- εモデルとする)と Large Eddy Simulation[52](以下，LES とする) である．k - εモデルは，N - S 方程式の変数の瞬時値を，平均量と変動成分に分割し，N - S 方程式に代入後，時間( アンサンブル) 平均することで得られるレイノルズ方程式中で未知の項であるレノルズ応力項を渦動粘性係数ν_t を用いて拡散項に統合し，さらに渦動粘性係数が乱流エネルギー k とエネルギー逸散率εで表現されるとしたモデルである．k - εモデルの基礎方程式を表 2 - 2 に示す．なお，

k- εモデルにおける N-S 方程式の拡散項は表 2-2 式(9)中の右辺第 2 項である．LES は，比較的流れ場に左右されない小さいスケールの乱れだけを統計的に処理し，

大きい渦は N - S 方程式を直接解くもので，k - εモデルに比べモデル化の度合が少ないため精度が高い．しかし，L E S は時間平均操作を施さないため非定常解析が必須であり，さらに大きいスケールの渦を表現するためには k - εモデルに比べて膨大な格子数を必要とすることから多大な計算時間を必要とする．よって本研究では，比較的高精度で，実用的な計算時間で解が得られる k - εモデルを用いることとする．

(23)

質量保存の式

x 0 U

j

j =

¶

運動量保存の式

(

₀

)

_i₃

i j j

i e j i

j j

i i g T T

x U x

U k x

3 2 P x x

U U t

U n b d

r - -

ïþ ïý ü ïî

ïí ì

÷÷ ø ö çç

è æ

¶ + ¶

¶

¶ + ¶

÷÷ø çç ö

è æ +

¶ - ¶

¶ = + ¶

¶

乱流エネルギーｋの保存の式

e

n ÷÷ + + -

ø ö çç

è æ

¶

= ¶

¶ + ¶

¶

k k j k j j

j P G

x k x

x kU t

k

粘性逸散率εの保存の式

(

^e

)

e n e

e e

e l k 3 k 2

j j

j

j C P C G C

k x x

x U

t ÷÷+ + -

ø ö çç

è æ

¶

= ¶

¶ + ¶

¶

エネルギー保存の式

÷÷ ø ö çç

è æ

¶

= ¶

¶ + ¶

¶

j T j j

j

x T x

x TU t

T n

拡散物質保存の式

÷÷ ø ö çç

è æ

¶

= ¶

¶ + ¶

¶

j c j j

j

x C x

x CU t

C n

表 2 - 2 　k - ε型２方程式モデル

・・・(8)

・・・(9)

・・・(10)

・・・(11)

・・・(12)

・・・(13)

ここに、

U ^：x_i方向の平均速度[m/s], P：平均圧力[Pa], k：乱流エネルギー[m²/s²],

e

：粘性逸散率[m²/s],

nt：渦動粘性係数[m²/s],

Pk：kの生産項, G_k：浮力生産項, 下添え字はEinsteinの総和規約に従う

n_t = C_D ^ke² ,

j i i

j j

i t

k x

U x U x

P U

¶

÷÷ ø ö çç

è æ

¶ + ¶

¶

=n ¶ , _i₃

i t

k x

g T G

T

s d b n

¶

= ¶ ,

t

e n n

n = + ,

T t

T s

n n

n = + ,

k t

k s

n n

n = + ,

e

e s

n n

n = + ^t ,

C t

C s

n n

n = + ,

9

T =0.

s , s_k =1.0, s_e = 1.3, s_C = 0.9,

09 0

C_D = . , C_l = 1.44, C₂ = 1.92, C₃ = 0

(24)

2 . 3 　制御法

2 . 3 . 1 　ファジィモデルの入出力変数

　空調空間の温度分布を決定する要因として，熱負荷およびその分布，吹出し風量，吹出し温度が挙げられる．本提案法では，まず熱負荷，熱負荷分布，吹出し温度，吹出し風量をパラメータとした定常 C F D 解析を行い，メッシュ分割位置での個々の空気温度を求める．なお，パラメータとしての熱負荷については，大空間空調の熱負荷の大部分を占める人体負荷のみを扱うこととする．これにより得られる変数は人体負荷，人体負荷分布，吹出し温度，吹出し風量，任意の３次元空間位置 x ，y ，z での温度とその空間位置 x ，y ，z の８変数となり，これらのデータを用いて多入力１出力のファジィモデルを構築する．

　次に外乱を熱負荷とした空調制御を考え，測定値や入力値から居住域平均温度を出力するファジィモデルを作成する．測定可能なデータを人体負荷分布，吹出し温度，吹出し風量，センサー温度とその位置とし，予測するデータを居住域平均温度とする．ここで，人体負荷分布( 客席位置) は使用用途毎にあらかじめ決められているため，人体負荷分布は別扱いとする．以上より，ファジィモデルは入力変数を吹出し温度，センサー温度，吹出し風量，温度センサー位置( ３次元成分) ，出力変数を居住域平均温度とする６入力１出力となる．さらに非等温流れ場の浮力の影響は空調に用いる温度域においては温度差にのみ依存するため，吹出し温度とセンサー温度との差は熱負荷と線形な関係となる．よって，入力変数を吹出し温度とセンサー温度との差，吹出し風量，温度センサー位置の５入力とすることができる．図 2 - 4 に C F D 解析データとファジィモデルの入出力変数を示す．

(25)

Travr

a ) 定常 C F D 解析で得られるデータ( 網掛けは C F D 解析の出力値）

T₀ Q x y z

Tr(x,y,z)

FM^* Travr

b ) ファジィモデルの入出力データ

dT₀(x,y,z) Q

x y z

FM^*

c ) 浮力による温度特性を考慮したファジィモデルの入出力データ

ここに、 q ：人体負荷

モード：人体負荷分布 T₀：吹出し温度

Q ：吹出し風量 x ：空間の位置座標 y ：空間の位置座標 z ：空間の位置座標

Tr(x,y,z) ： (x,y,z)での温度 Travr ：居住域平均温度 dT₀(x,y,z) ： =T₀-Tr(x,y,z) q モード T₀ Q x y z Tr(x,y,z)

( 網掛けデータは前処理で作成）

図 2 - 4 　C F D 解析データとファジィモデル入出力変数

* F M は各モードごとに作成( F M : ファジィモデル)

(26)

2 . 3 . 2 　制御アルゴリズム

　提案する制御法( 以下，F M - P I D 制御とする) は，ファジィモデルを用いてあらかじめ居住域の平均温度を推定し，居住域平均温度が目標値となるように制御を行うもので，用いるファジィモデルは前項の入力変数を計測または入力することにより居住域平均温度を推定する．居住域平均温度推定後の制御部はファジィ制御をはじめ種々の方法が考えられるが，既設建物への導入の容易さから P I D 制御とする．図 2 - 5 に F M - P I D 制御の制御ブロック図を，表 2 - 3 に P I D 制御の基本式を示す．

PID 制御部制御対象

検出部ファジィ推論部

目標値

外乱( 人体負荷)

Trave T₀

Tsns (X_Tsns,Y_Tsns,Z_Tsns)

T’rave +

-

ここに、 T₀：吹出し温度 Q ：吹出し風量　　 Tsns ：室温センサー温度

(X_Tsns,Y_Tsns,Z_Tsns) ：室温センサー位置 Trave : 居住域平均温度　　

T’rave ：居住域平均温度推定値 Q

図 2-5　FM-PID 制御ブロック図

表 2-3　PID 制御

・・・(14)

( ) ( ) ( ) ( )

÷÷ø çç ö

è

æ + +

= ^K ^e ^t _T

ò

^e ^d ^T ^de_dt^t

t

m ^t _D

I

p 0

1

t t

ここで、

t：時間

( )

t

m ：制御出力, e

( )

t ：偏差 (＝目標値－現在値), Kp：比例ゲイン, T_I：積分時間, T_D：微分時間

(27)

第３章　定風量方式

3 . 1 　緒言

　定風量(以下，CAV とする． CAV：Constant Air Volume)方式は最も古典的な空調方式であり，大空間やクリーンルームなどの比較的特殊な空間に現在も多く採用されている．定風量方式は吹出し風量を一定とし，吹出し温度を制御する方式で，以下のような利点がある[ 5 3 ] ．

　・換気量を定常的に十分確保できる．

　・風量が十分かつ安定しているため，温湿度・気流分布・空気質等に関し安　　定した質の高い空調ができる．

　・設備費が廉価で保守管理が容易である．

また，欠点としては以下が挙げられる．

　・室温制御用温度センサーの設置位置が決めにくい．

　・年間を通して最大負荷に対応可能な送風量で空調するためファンの搬送動　　力消費量が大きい．

　このような理由から，ほとんどの大空間では C A V 方式が採用されている[ 5 4 - 5 7 ] ．また，吹出し温度の制御は温度センサーを用いた P I D 制御が一般的であるが，温冷感指標である S E T^*制御や P M V 制御が一部で試みられている[ 5 8 ] ．　そこで，本章では大空間で最も一般的に用いられている C A V 方式において，制御対象をドーム球場とした場合の，F M - P I D 制御の制御性を P I D 制御と比較することとする．

(28)

3 . 2 　C A V 方式の制御

　図 3 - 1 に C A V 方式の空調機周り制御システム( 温度のみ) を示す．室内あるいは吸込み口に取り付けられた温度センサーの計測値を用いて調節器( 図中，D D C ) にて冷水コイル，温水コイルの電動二方弁( 図中，M V ) の開度を調整し，吹出し空気温度を制御する．なお，図中には示していないが，換気量確保のために一定量の外気導入と，排気を行う必要があり，外気導入及び排気用のダクトにはダンパが必要である．また，コイル能力の選定の際，吹出し温度の上下限値を設定する必要があり，本制御シミュレーションでは通常の空調設計で一般的に用いられる値である設計温度± 1 0 ℃とする．

図 3 - 1 　C A V 方式制御システム

（温度のみ、発停・湿度その他の制御を含まず）

DDC

MV MV TE

C C

H C

冷水温水

FAN 給気

還気外気

排気

ここに、 TE ：室内温度検出器 DDC ：デジタル調節器 MV ：電動二方弁 CC ：冷水コイル HC ：温水コイル FAN ：送風ファン

(29)

図 3- 2 　制御対象

吹出し口

吸込み口

吸込み口温度センサーc 温度センサーw

温度センサーa スタンド部

アリーナ部２０ｍ

３０ｍ

８ｍ８ｍ８ｍ

１０ｍ x

z y

3 . 3 　制御対象及びシミュレーション条件

3 .3.1 　制御対象

　図 3 - 2 に示すドーム球場の１スパンを制御対象とする．寸法は 5 0 m ( L ) × 1 0 m ( W )

× 2 4 m ( H ) で，スタンド部とアリーナ部を設ける．空調は吹出し口をスタンド部後方に 1 箇所，吸込み口を天井とスタンド前方に各 1 箇所設置し，吹出し方向は水平とした．各制気口の風量・寸法を表 3 - 1 に示す．温度センサーは無線センサーを想定してアリーナ部中央に設けた場合( 図 3 - 2 中，温度センサー a ），スタンド部に設けた場合( 図 3 - 2 中，温度センサー c ），スタンド部後方の壁面に設けた場合( 図 3 - 2 中，温度センサー w ）の３通りとする．さらに使用用途として図 3 - 3 に示すモード A ( 競技観戦時) ，モード B ( コンサート時) ，モード C ( 小集会) の３通りの人体負荷分布を想定する．また空調期間を冬期，設定温度を 2 2 ℃とし，人体負荷の増加とともに暖房から冷房へと変化する空調空間の制御性を検討することとする．なお，モード A の居住域はスタンド部床上高さ 2 . 0 m まで，モード B の居住域はスタンド部とアリーナ部の床上高さ 2 . 0 m まで，モード C の居住域はアリーナ部高さ 2 . 0 m までとする．

(30)

図 3 - 3 　人体負荷発生位置　 (a) モード A

(b) モード B

(c) モード C 人体負荷人体負荷

人体負荷

寸法 [m] 風量 [m³/h]

吹出し口 2.0×1.0 18,750 吸込み口(上部) 2.0×1.0 3,750 吸込み口(下部) 2.0×1.0 15,000

表 3 - 1 　吹出し・吸込み条件

(31)

期間冬期

用途モードA, モードB, モードC 人員密度[人/m²]^† 0.0,0.5,1.0,1.5,2.0,2.5,3.0

熱通過負荷[W/m²] 50 (天井部)

空間分解能 18(x)×5(y)×22(z) = 1,980 セル　注 ^† 人体負荷発熱量は70 W/人(顕熱のみ)とした。

3.3.2　CFD 解析条件

　ファジィモデル構築用データの作成と制御シミュレーションはすべて C F D 解析により行うこととする．C F D 解析は，2 . 2 . 2 ，2 . 2 . 3 項で示したとおり乱流モデルをｋ - ε型 2 方程式モデル，解法を構造格子系の有限体積法による S I M P L E C 法とし，差分スキームは運動方程式の移流項のみ Q U I C K スキームとする．計算メッシュ数は 1 ，9 8 0 セル，壁面境界条件は図 3 - 1 に示す X Z 面の両側とアリーナ側の Y Z 面にフリースリップ，他は対数則とし，人体負荷は床上第 1 セルの流体中に体積発熱させる．制御シミュレーションにおける非定常解析時の時間間隔は，クーラン条件[ 5 9 ] をほぼ満足する 0 . 2 5 s とする．

3 . 3. 3 　ファジィモデル用データ

　C F D を用いた定常解析により，ファジィモデル構築用のデータを作成する．

モード A ，モード B ，モード C それぞれについて，人体負荷を 0 . 0 ～ 3 . 0 人 / m²までとした 7 通りの C F D 解析を行う．人体負荷は顕熱のみを考慮し，7 0 W / 人[ 6 0 ] を，

冬期の想定として暖房熱通過負荷 5 0 W / m²を天井面に一定値として与える．C F D 解析条件を表 3 - 2 に示す．

表 3 - 2 　ファジィモデル作成用 C F D 解析条件　

(32)

3 .3 .4　ファジィモデル

　前項で得られた C F D 解析結果からファジィモデルを構築する．まず，制御対象の特性に合わせて不要な入力変数を削除する．必要な温度分布は居住域のみであることから，3 次元成分の温度分布データを床面第一セルのみの 2 次元データとして抽出し，さらに各温度センサー位置の y 座標が等しいことから，温度センサー位置は x 方向のみとなる．また，C A V 方式では吹出し風量は一定であるため，

最終的に 2 . 3 . 1 項に示した 5 入力 1 出力ファジィモデルは，センサー温度と吹出し温度の差，温度センサー位置 x の 2 つを入力変数，居住域平均温度を出力変数とした 2 入力 1 出力のファジィモデルとなる．モデリングの終了条件は試行錯誤の結果から，後件部変数の最大誤差が 2 % 以下となった場合とする．なお，誤差は後件部変数の最大値から最小値を引いた値で標準化した値である．得られたメンバーシップ関数を図 3 - 4 に示す．図に示されるメンバーシップ関数が密に分割されている部分が非線形性の強い出力データが存在する部分であり，居住域平均温度と吹出し温度との差に対してはすべてのモードで冷房側に多い．X 方向位置に関しても全モードでスタンド部側に多いことから，冷房時のスタンド部で非線形性の強い出力データが存在していることが分かる．このようにメンバーシップ関数の分割配置からデータの大局的な特性が把握できることがファジィモデルの特徴の一つとなっている．図 3 - 5 ～ 7 に各モードにおける定常 C F D 解析とファジィモデルの入出力関係を示す．モード A での定常 C F D 解析結果( 図 3 - 5 ( a ) ) は，それぞれの位置( 以下，x とする）において吹出し温度と x での温度との差( 以下，d T 0 とする）が 4 ℃の時，吹出し温度と居住域平均温度との差( 以下，d T r とする）が 4 ℃で同じとなっており，x における d T 0 と d T r の関係はほぼ線形となっている．

直線の傾きは x が 8 . 7 5 m でほぼ 1 となっており，x が小さいほど小さい値を示す．

これらの入出力関係に対し，ファジィモデル出力値( 図 3 - 5 ( b ) ) もほぼ同様の傾向を示し高精度に入出力関係をモデリングしているが，d T r が - 8 ℃付近でのすべての x の d T 0 が - 10 ℃となっていること，x が 1 6 .2 5 ～ 2 8 .7 5 m で d T 0 が - 5 ℃付近で線形にならないことが相違点として挙げられる．前者は初期メンバーシップ関数を固定値としてモデリングを行ったためであり，後者は x が 1 6 . 2 5 m 以上，d T 0 が約 - 5 ℃以上，d T r が - 5 ℃以下の領域にデータが無いためにメンバーシップ関数を

(33)

0 1

-6.35 4.13

0 1

1.25 48.33

0 1

-15.42 4.13

0 1

1.25 48.33

0 1

-12.79 4.08

0 1

1.25 48.33

dT

x

dT

x

dT

x

( a ) 吹出し温度と位置 x での温度との差 ( モード A )

(b) 位置 x (モード A)

( c ) 吹出し温度と位置 x での温度との差 ( モード B )

(d) 位置 x (モード A)

( e ) 吹出し温度と位置 x での温度との差 ( モード C )

(f) 位置 x (モード A)

図 3 - 4 　各前件部変数のメンバーシップ関数

(34)

-8 -6 -4 -2 0 2 4 6

-12 -10 -8 -6 -4 -2 0 2 4 6

吹出し温度と位置xでの温度との差[℃]

吹出し温度と居住域平均温度との差 (dTr)[℃]

1.25 3.75 6.25 8.75 11.25 13.75 16.25 18.75 21.25 23.75 26.25 28.75 31.67 35 38.33 41.67 45 48.33 位置x[m]

(a) 定常 CFD 解析

( b ) ファジィモデル

図 3 - 5 　モード A における入出力関係

-8 -6 -4 -2 0 2 4 6

-12 -10 -8 -6 -4 -2 0 2 4 6

1.25 3.75 6.25 11.25 16.25 18.75 21.25 23.75 26.25 28.75 38.33 48.33 位置x[m]

分割していく過程で x = 1 6 . 2 5 m の値を参照しながら後件部実数値を線形補間して求めたことが原因と考えられる．この領域には入出力データが存在しないため計算上誤差として検出することは困難である．本研究では x が 1 6 . 2 5 ～ 2 8 . 7 5 m での温度センサー位置は設定していないため制御シミュレーションには影響を与えないが，精度を求めるために収束判定条件を厳しくしてメンバーシップ関数を密に分割すると，逆にこのように入出力データが存在しない場合が生じ，検出困難な誤差を生じることに留意する必要がある．モード B での定常 C F D 解析結果( 図 3 - 6

(35)

-20 -16 -12 -8 -4 0 4 8

-24 -20 -16 -12 -8 -4 0 4 8

1.25 3.75 6.25 8.75 11.25 13.75 16.25 18.75 21.25 23.75 26.25 28.75 31.67 35 38.33 41.67 45 48.33 位置x[m]

(a) 定常 CFD 解析

( b ) ファジィモデル

図 3 - 6 　モード B における入出力関係

-20 -16 -12 -8 -4 0 4 8

-24 -20 -16 -12 -8 -4 0 4 8

1.25 3.75 6.25 8.75 16.25 18.75 21.25 23.75 38.33 48.33 位置x[m]

( a ) ) は，d T 0 と d T r が 0 ℃の位置で交差しており，傾きは x が小さいほど大きくなっている．また，交差点において線形性が崩れており，その傾向は x が小さいほど大きい．ファジィモデル出力値（図 3 - 6 ( b ) ) も上記の傾向を再現しており，

本モデリング手法が入出力データの関係を高精度に表現できることが分かる．ただし，x が 1 . 2 5 m で d T 0 が - 1 3 ℃付近の線形性が著しく定常 C F D 解析値と異なっており，これはモード A と同じ理由によると考えられる．なお，x が 1 . 2 5 m の位置は温度センサー w であるが，実用的な空調吹出し温度差は 1 0 ℃より小さいため，

電気通信大学大学院電気通信学研究科 博士(工学)の学位申請論文

多目的大空間におけるファジィモデルを用いた 空調制御に関する研究

植田 俊克