コンテンツへスキップ

Scilab/Xcosが必要になったので入れて使ってみます。動作の確認として、MATLABでやったのと同じ解析をしてみました。

Scilab/Xcosが必要になったので入れて使ってみます。windows10,Scilab6.0.2です。動作の確認として、MATLABでやったのと同じ解析をしてみました。

目次

  • 導入
    • scilabとは
    • Xcosとは
  • インストール
  • 使ってみる
    • ブロックを並べる
    • 設定
    • シミュレーション
    • 構文
      • polyとは
      • syslinとは
      • bodeとは
      • dampとは
  • MATLAB/SimlinkとScilab/Xocsの比較
    • MATLAB/Simlink
    • Scilab/Xcos

導入

Scilabとは

Scilabはフリーのオープンソフトウェアで、端的に言ってしまえばMATLABのフリー版です。様々な数値解析ができます。

詳しくは公式にて

MATLABに関してはこちらの記事をご覧ください。=>MATLAB使い方個人的まとめ

Xcos

動的な連続/離散時間システムのモデリングやシミュレーションができます。 MATLABにとってのSimlinkです。グラフィカルにブロックを並べることでモデリングでき、コンパイルやシミュレーションもまとめてできる優れものです。

X2C(https://www.lcm.at/virtuelle-entwicklung/x2c/)

インストール

公式ページ: Download Scilab 6.0.2 -scilab

言語設定でjapaneseを選択し、デフォルトのままfullバージョンでインストールします。

インストールページのスクショ

使ってみる

MATLABと同様に、コマンドウィンドウ、Xcosウィンドウ、SimNoteのウィンドウをそれぞれ使って解析を行います。アプリケーション/Xcosで起動するとXcosウィンドウとライブラリが現れます。

Scilabの画面
Xcos

ブロックを並べてブロック線図を作り、sin波を表示させてみます。

ブロックを並べる

まずは信号源、出力/表示などからsin波発生器、CSCOPE、CLOCKを見つけて右クリックかドラッグでXcos上に配置し、矢印をドラッグして配線します。 スコープは赤い矢印があるものを選びましょう。これはイベントシミュレーションを駆動するイベント信号になります。

ブロック図の作成

設定

それぞれのブロックで右クリックすると項目が現れて、その中からブロックパラメータを開きます。 CFSCOPEを設定します。[参考1]

  • Color or mark vector: スコープで表示するグラフ線の色またはマークの連続によるグラフを指定する。入力信号数にかかわらず8セットの数字が必要。数字が正のとき、連続線で色が指定でき、負のとき、マーク形状が指定できる。(plot2dコマンドに準ずる)
  • Output window number: グラフ・ウィンドウ番号の指定。
  • Output window position: グラフ・ウィンドウの表示位置指定。左上の位置をx,yのドットで示す。
  • Output window size: グラフ・ウィンドウの表示サイズをドット数で指定。
  • Ymin: グラフ・ウィンドウのY軸最小値
  • Ymax: グラフ・ウィンドウのY軸最大値
  • Refresh period: グラフ・ウィンドウのX軸(時間)最大値 シミュレーション設定時間がこの値より大きい場合、グラフをクリアし再度続きの時間から描画する。基本的使い方では、シミュレーション時間と一致させる。
  • Buffer size: 入力データを一時的に蓄えておくメモリ領域のサイズを指定。グラフ表示速度に大きく影響するので、変化を観察したい場合以外は、少なくともRefresh period以上に設定する。最小は2。
  • Accept herited events: (inherited?) 継承の無効(0)・有効(1)を指定。1にするとイベント・ポート赤が削除される。アクティベーションが上流で行われる場合、継承が可能。

SCOPEの積算時間を10に、最小値~最大値を-10~10にしました。 クロックとsin発生器のパラメータも設定します。初期化時間を0.0に、ゲイン(振幅)5にします。

シミュレーションの設定

シミュレーション

シミュレーションを開始します。シミュレーション/パラメータ設定から、[参考1]

  • Final integration time: シミュレーション終了時間[単位は任意]
  • Real time scaling: 正の実数nを入力すると、シミュレーション設定単位時間tがn×t秒としてシミュレーションが実行される。n=1とした場合、シミュレーションの設定単位時間が1秒として解析表示される。
  • Integrator absolute tolerance: ソルバーの絶対誤差公差
  • Integrator relative tolerance: ソルバーの相対誤差公差
  • Tolerance on time: ODE(常微分方程式)ソルバーの最小時間公差
  • Max integration time interval: ソルバー呼び出しの最大時間間隔。"too many calls"というエラーメッセージが出た場合、値を小さくする。(解析精度を落とす。)
  • Solver 0(CVODE) - 100(IDA): 数値ソルバーの選択。
  • maximum step size: ソルバーがとる最大時間ステップ

最終積算時間を1.0E1.0、すなわち10sにしました。 AEBならA*10^Bになります。

document/test_0620/image.scgはこのようになりました。指定した通りに動いていることが確認できます。

document/ test_0620/image.scg の出力

構文

ビギナーのための科学技術計算ソフト(Scilab)の使い方講座 3. 付録: Scilab言語の基本構文 -SimCircuit

Scilab コンソールの使い方の基本

「;」を後ろにつけると結果を表示しないみたいですね。

polyとは

[参考2]によれば、ploy関数は以下のように書きます。

p=poly(実数もしくは行列,"変数")
-->a=[-3 -4];
-->p=poly(a,"x")
 p  = 

            2
   12 +7x +x 

コンソールに上のように入力すると、12+7x+x^2を返します。((x+3)(x+4)=0, のときx=-3,-4) pを多項式変数として扱う事ができます。

-->s=poly(0,'s')
s =

    s
polyをコンソールに入力したときの出力

syslinとは

[参考3]によれば、syslin関数は以下のように書きます。

q=syslin(連続時間システムc or離散時間システムd,多項式)
--> s=poly(0,'s');

--> P=syslin('c',1/(s+1))
 P  = 

           
     1     
   ------  
           
   1 + s   

bodeとは

[参考4]によれば、syslin関数は以下のように書きます。

bode(syslinで定義された線形システムの多項式,x軸の最小値,x軸の最大値,x軸の刻み,['コメント'])

コメントにタイトルなどをかきます。 こだわらないなら関数のみ与えればOKです。

dampとは

[参考5]によれば、syslin関数は以下のように書きます。

[wn,z] = damp(連続時間のsys関数)

MATLAB/SimlinkとScilab/Xcosの比較

某授業のレポート1でやった問題の伝達関数を使って確認していきます。

問題: オープンループの伝達関数L(s)が以下の時の、回ループと閉ループでのボーデ線図をそれぞれ表示します。

$$L(s)=\frac{0.1}{s(s+0.5)(s+1)}$$

MATLAB

  • オープンループ

$$L(s)=\frac{0.1}{s^3+1.5s^2+0.5s}$$

MATLABオープンループでのボーデ線図

位相余裕:59° 利得余裕:17.5dB

  • フィードバック $$L(s)=\frac{0.1}{s^3 + 1.5s^2 + 0.5s + 0.1}$$

A. 周波数応答

MATLABフィードバックでのボーデ線図

B. ステップ応答

設定がうまくいかず、ぎこちないですがこんな感じ。最終値は1です。

MATLABフィードバックでのステップ応答

Scilab

  • オープンループ

複数行のコマンドを何度も使用する場合、アプリケーション/SciNoteに書くとスクリプトのようにまとめて実行できます。 見やすさのために出力なしで実行しました。

SciNoteを使用している様子

A. 周波数応答

まずはコマンドで実行してみます。

コード

s=poly(0,'s');
p=syslin('c',0.1/(s*(s+0.5)*(s+1)))
bode(p)
show_margins(p)
[gm,fp]=g_margin(p) disp("ゲイン余裕[dB]") disp(gm) disp("位相交差角周波数[Hz]") disp(fp) 
[pm,fg]=p_margin(p) disp("位相余裕[°]") disp(pm) disp("ゲイン交差角周波数[Hz]") disp(fg) 

出力

ゲイン余裕[dB]

   17.501225

 位相交差角周波数[Hz]

   0.1125395

 位相余裕[°]

   59.289834

 ゲイン交点角周波数[Hz]

   0.0293667

ゲイン余裕:位相が-180°と交差する周波数(位相交差角周波数)においてゲインが0dBより何dB小さいか

位相余裕:ゲインが0dBと交差する周波数(ゲイン交差角周波数)

オープンループボーデ線図は

Scilab オープンループでのボーデ線図

MATLABと一致してますね。

B. ステップ応答

Scilab オープンループでのステップ応答
  • フィードバック

A. 周波数応答

コマンドからやります。フィードバックした場合の伝達関数G(s)は以下のようになるので

$$G(s)=\frac{L(s)}{1+L(s)}$$

コード

s=poly(0,'s');
p=syslin('c',0.1/(s*(s+0.5)*(s+1)))
g=p/(1+p)
disp(g)
bode(g)
show_margins(g)
[gm,fp]=g_margin(g) disp("ゲイン余裕[dB]") disp(gm) disp("位相交差角周波数[Hz]") disp(fp) 
[pm,fg]=p_margin(g) disp("位相余裕[°]") disp(pm) disp("ゲイン交差角周波数[Hz]") disp(fg) 
[wn,z] = damp(g); disp("固有周波数[Hz]") disp(wn) disp("減衰係数") disp(z) 

出力

                           
            0.1            
   ----------------------  
                    2   3  
   0.1 + 0.5s + 1.5s + s   

 ゲイン余裕[dB]

   16.258267

 位相交差角周波数[Hz]

   0.1125395

 位相余裕[°]

   114.71226

 ゲイン交差角周波数[Hz]

   0.0313482

 固有周波数[Hz]

   0.2964606
   0.2964606
   1.1378001

 減衰係数

   0.6108736
   0.6108736
   1.

同じ曲線を描いていることを確認しました。

Scilab フィードバックでのボーデ線図

B. ステップ応答

最終値1で、100秒まで計算しています。

Scilab フィードバックでのステップ応答

コードは[参考6]を参考にさせて頂きました。

雑記

  • Xcosのライブラリで固まったらEscキーを押しましょう!ドラッグより右クリックの方がいいです。
  • 周波数応答を調べるとき、関数に伝達関数を与えています。なのでXcosでブロック線図を描こうと思ったらブロックを並べて伝達関数を計算させて、それを取得して関数に入れる、という手順を踏むと思うので回りくどいです。ここでやった例はすぐ伝達関数を計算できるので、コマンドでやった方が早いと思いました。
  • 何気に変数ブラウザ便利ですね。

参考

1.ビギナーのための科学技術計算ソフト(Scilab)の使い方講座 2.3 Xcosによるシミュレーション・ツールとしての使い方 (GUI型ブロック線図入力モジュール) -SimCircuit

2. Scilab関数の説明:polyとは

3. Scilab関数の説明:syslinとは

4. Scilabでのボード線図の描き方

5. damp -scilab

6. scilabで制御系設計

無料磁場解析ソフトFEMMチュートリアルで使い方を確認しつつその時の記録を残しておきます。

チュートリアルで使い方を確認しつつその時の記録を残しておきます。

目次

  • FEMMとは
  • 動的磁場有限要素法とは
  • インストール
  • 磁場解析をしてみる:空心コイルソレノイドの磁場解析
  • 解析結果
  • まとめ
  • 雑記
  • 参考

FEMMとは

低周波数の磁気学と静電気学について,2次元平面と軸対象問題を解くためのソフトウェアです。Windows上で動作します。

ソルバという、複数の変数を含む数式において目標とする値を得るための最適な変数の値を求めることができる機能がついています。FEMMのソルバは 動磁場有限要素法を使用しているものみたいです。

Kelvin変換,Robin境界,SDI境界,IABC境界など多数のOpen boundary 境界が実装されている、みたいですが、解析方法とかよくわかりません。

チュートリアルにはLuaという言語と統合できるみたいなことが書いてありますが、PythonでもpyFEMMpyfemmを使って動かせるみたいです。 FEMM自体はC言語で動いているみたいです。

なお、この記事では磁場に関しては深く言及しません。

 ソルバー -weblio辞書

動的磁場有限要素法 とは

Finite Element Method (有限要素法)は、微分方程式を近似的に解くための数値解析方法です。複雑な形状。性質を持つ物体を小部分に分割することで近似し、全体の挙動を近似して計算機に計算させようとする方法です。要素の選択とメッシュの作り方を近似していきます。 構造力学や流体力学などの様々な分野で使用されているみたいです。

mesh

要素の円卓とメッシュの作り方を指定してあげます。 FEMMは固体要素を扱うことが多く、メッシュを正三角形で行っています。一般的に、メッシュを細かくすればするほど計算精度が良くなりますが、その分計算時間が掛かります。また、三角形よりも四角形、の方がj計算精度がいいです。

端的に言えば、各要素で編微分方程式を解いてるみたいです。(詳細についてしゃべるのを避けた)

インストール

Finite Element Method Magnetics :Download

からダウンロードします。

download画面

磁場解析をしてみる:空心コイルソレノイドの磁場解析

以下がバージョンが違いますが、チュートリアルです。

FEMM 3.4 Magnetostatic Tutorial1

研究 FEMM チュートリアル.pdf -Evernote

日本語に翻訳してあるものを読みながら、重要な部分だけやります。

今回の解析対象はこれです。

解析対象

2Dの解析ソフトなので、FEMMで一回に解析するのはこれの断面になります。これは軸対称のモデルなので描くのは左側の四角になります。

新規モデルの作成

File/NEWでMagnetic problemを選択します。静磁気問題になります。次に問題の種類を指定していきます。axisymmetric(軸対称)を選択し、Length Unitをmm、Frequencyを0Hzにします。

平面で解析したい場合はPlanarを選択します。

境界線、コイルの描画

有限要素ソルバは対象のオブジェクトを含む空間の有限領域を指定して解析するみたいです(?)。なので、空間を指定していきます。 視界の設定をします。View/Keyboadから、Bottom:-4,Left:0,Top:4,Right:4に指定します。スクリーンは指定された領域を含む 最も小さい長方形に設計しなおされます。

次にOperate on nodesを選択します。小さいブロックボックスで左端にあるこれです。ノードを(0,4),(0,-4),(0,0)にマウスから設置します。Tabを使用すると、キーボードから入力できます。

次に、Operate on segmentsを選択します。左から2番目です。(0,4),(0,-4)をクリックすると線が引かれます。画像は見にくかったので変更しました(上記参照)。 次に、Operate on arc segmentsを選択します。

左から3番目です。(0,4),(0,-4)をクリックすると軸対称な円弧が描かれます。propertiesが現れますが、そのときArc Angleを180、Max segment Degreeを25にします。FEMMでは曲線は直線を繋げたものとして表されるので、解像度のように考えればよいでしょう。こおでは2.5で充分であるとします。

もし失敗した場合はEdit/Undoか、上の切り返す矢印マーククリックして選択した後にメニューのX印を押しましょう。

同じ手順で内側の四角を描きましょう。Operate on nodesを押して、(0.5,-1)(1.5,-1)(1.5,1)(0.5,1)にノードを置く。 その後Operate on segmentsで線で結びます。

材料や特性の関連付け

Operate on Block Labelsを選択します。コイル領域とコイル領域の外側の空気部にブロックラベルを置きます。

このプログラムは、問題を形状、材料、特性によって解くので、それらを関連付けていきます。

メニューからProperties/Materials Library を選択しますこの中から、Air とCopper AVW Magnet wire/18 AGWを右のModel Materialsにドラック&ドロップで加えます。OKを押します。

コイルに回転属性を追加します。

properties/CircuitsからAdd Propertyを円卓し、new circuitに適当な名前をつけます。Seriesを選択して巻き線領域になるように指定します。これは時間調和問題に使用される設定です。 Circuit Currentに1を入力します。

円の中、四角の外のブロックラベルを右クリックで選択し、スペースか上のバーの指を指しているようなマークでpropertyを開きます。

Block typeを先ほどのAirにし、Let Triangle Choose Mesh Sizeのチェックを外して0.1にします。このとき、1辺最大0.1の正三角形で内部を満たそうとします。

四角の中はBlock type 18 AVW、Mesh size 0.1、In circuit Coil、Number of turn 1000を入れます。これで、時計回りに1000巻きの18材のコイルが指定されます。右ねじの方向が正。-1なら逆時計回りに1。

境界条件

境界はboundaryです。Prorerties/Boundaryを選択し、Add Propertyを選択し、nameをABCに、BC typeをMixesにします。ABCは(チュートリアルによると)境界の無いオープンスペースのイ ンピーダンスに近似した,-漸近的な境界条件. を作成している事を意味するらしいです。

c0 coefficient(係数) と c1 coefficientは以下で表されます。

$$ \frac{1}{\mu_r\mu_o}\frac{\partial A}{\partial n} + c_0 A + c_1 = 0$$

$$A:磁気ベクトルポテンシャル、\mu_r:境界に近い領域の比透磁率$$

$$\mu_o:真空の透磁率、n:境界に対するの標準の方向$$

ここでは、以下のようにします。

$$c_0=\frac{1}{\mu_r\mu_o R}$$ $$c_1=0$$

Rは球の外径です。今回の入力は、c0 coefficient に 1/(uo4mm)、c1 coefficientに0を入力します。uoに対する真空の透磁率の数値やmに当たる0.001を代入し、OK、OKをクリックします。 ちなみに1inch=0.0254mです。 空間の有限領域に着いて、無限境界空間にあるコイルによって作成されたフィールドをモデル化出来ました。

解析の開始

境界条件を指定するためにOperate on arc segmentsモードにします。外側の境界を意味する円弧を右クリックで選択し、スペースでpropertiesを開きます。Boundary condにABCを選択します。 バーの黄色のMesh generaterをクリックします。すると、このように塗りつぶされます。

FEA(有限要素法)による解析を実施します。turn the crankを押して解析を行ってください。 今回の解析は一瞬でしたが、プログレスバーが動いていない様子であった場合は境界条件の指定がうまくいっていないかも知れないので、条件をもう一度見直す必要があるかも知れません。

解析結果

眼鏡のアイコンをクリックすると解析結果がみれるタブが現れます。上のバーにポイント、等高線(コンター)、エリアを選択できます。

一点の値

左クリックで任意の点をクリックすることで、そのポイントの様々なフィールドプロパティがFEMM outputウィンドウに現れます。TABでより精密にポイントを選択できます。 消してしまった場合はView/outputで再度表示出来ます。上の方(0,4)真ん中右の方(3,0)を指定したときはこんな感じでした。

(0,4)
(3.0)

コイルの末端特性

コイルアイコンをクリックします。Coilプロパティはこのようになります。問題が線形で、一つの電流しかないので。Flux/Currentの結果はコイルのインピーダンスと解釈できます。Flux/Current = 0.897763 mH、Voltage/Current =R = 0.13155 Ωですね。

等高線(コンター)に沿ったフィールド値のプロット

コンターアイコンを押し、(0,4)(0,0)(0,-4)付近をクリックしてグラフアイコンをクリックします。デフォルトでは磁束密度の振幅がプロットされますが、ほかにも選択できます。 ポイントが選ばれた順番に、進行方向からみて↑右側を選択するようなプログラムになっています。

磁束密度のプロット

デフォルトでは白黒ですが、ツールバーでカラーにできます。

雑記

jMAGという有料のが研究室にあるので、使っていきたいですね。これはノートに入れられますが、jMAGは3D解析が出来るので。今後は下をやりたいですね。

  • pythonで様々なパターンを自動的に解析できないか探す
  • 他の解析ソフトを使用して今回の結果を検証

参考