JDKインストールして、
環境変数指定して、
コード書いて、 javacでコンパイルして（クラスファイル作って）、 javaで実行→エラーでできなかった。

原因と思われるエラー

おれのエラー

Exception in thread “main” java.lang.UnsupportedClassVersionError: hello : Unsupported major.minor version 60.0

こういうエラーは、javacとjavaのバージョンが違ったときに出るらしい。
実際下記のようにバージョンを確かめるコマンドを打つと、javacとjavaのバージョンが違った。

>javac -version
javac 16.0.1

>java -version
java version “1.6.0_26”
Java(TM) SE Runtime Environment (build 1.6.0_26-b03)
Java HotSpot(TM) 64-Bit Server VM (build 20.1-b02, mixed mode)

16と1.6だそうです。違いますね。

さっき手動でインストールした（解凍したファイルを移動しただけ）～\jdk-16.0.1\binには、16.0.1.0が入っている。
だから、16が正だ。

1.6 はどこから出てきたの？

原因の発見

コマンドラインで、where javaしたら原因が見つかった。
whereについてはこのサイト。
下記のように、コマンドラインでjavaっていったら、2つあるらしい。

>where java
C:\Windows\System32\java.exe
C:\pg\jdk-16.0.1\bin\java.exe

多分上の方を優先して使ってるんだと思う。

javacは一つみたい。

>where javac
C:\pg\jdk-16.0.1\bin\javac.exe

入れた覚えないんだけどな？
なんか変な依存関係があって、ある日突然正常に動かなくなったら嫌だな。

　対処

とりあえず、環境変数の
%JAVA_HOME%\binを
%SystemRoot%\system32よりも上に持ってきた。
上に書いてあるほど、前に書いてあって、優先されるようです。

下記のように、コマンドラインでjavaと入力したら、16が使えるようになった！

>where java
C:\pg\jdk-16.0.1\bin\java.exe
C:\Windows\System32\java.exe

>java -version
openjdk version “16.0.1” 2021-04-20
OpenJDK Runtime Environment (build 16.0.1+9-24)
OpenJDK 64-Bit Server VM (build 16.0.1+9-24, mixed mode, sharing)

…Androidアプリを作りたくて、Android Studioの勉強の前に、JDKでお勉強してるんだが…
JDKを導入して、環境変数変えたけど…
Android Studioのjava.exeとか、javac.exeとかに影響しないのかな…
java.exeも、javac.exeもどうやら、C:\Program Files\Android\Android Studio\jre\binに入っているらしいので、androis studioはそれを使っているんだろうか？
android studioは、コマンドラインでjavaと打って起動するjava（system32に入っているやつとか）を使ってないのかな？。そもそもsystem32にはいつ何で入ってたんだろう？

Hello worldできた。

>java hello
Hello.

電気回路シミュレーターのLTSpiceを使って、時間が経つと電気抵抗が上がったり下がったりする素子を作りたい。
操作すると抵抗値が変化する可変抵抗とかではない。
時間が変数となって、値が上がり下がりする抵抗値を入れたかったのだ。

検索したら見つけたページ

時間が経つにつれて上がる抵抗は、
抵抗Rを、変数（時間time）を使って式で表す。

valueに

R=10000-100*time

と入力すれば、10kohmから、1秒ごとに100ohm減っていく抵抗を表現できる。

R=を忘れたら動かなかった気がする。

R=10000+100*time

と入力すれば、10kohmから、1秒ごとに100ohm増えていく抵抗を表現。

静電容量Cの経時変化

さっきのコンデンサ版。
こちらもやはり、LTSpiceを使って、時間が経つと静電容量が上がったり下がったりする素子を作りたい。
時間が変数となって、値が上がり下がりする静電容量を入れたかったのだ。

こちらも上のページと同じサイトで、解決策を見つけた。

なぜかCを直接入力できず、電荷量で指定することになる。

時間が経つにつれて上がるコンデンサ（キャパシタ）は、 電荷量Qを、変数（電圧xと時間time）で表すことになる。
電圧xを入れないと、静電容量っぽい振舞ではなくなる。
Q=CV

普通の1uFのコンデンサを作りたければ

Q=1u*x

xの所には電圧が入るので、入れ忘れてはいけない。

Q=静電容量*xと書く感じ。

Q=1u*x+1n*time*x

一秒あたり1 nFずつ容量が増えていく。
Q=CV
図では2 Vを印加しているので、一秒ごとに2nクーロン電荷が多くたまるようになるから、実質2 nAの電流が流れる感じになっている。

三角波を式で書く

三角波についてはPWLを使うのがお約束らしいが、だめだ。
GUIを使って値を指定するのが気に入らない。（大変。）

式を使って書く。（こっちの方が大変だった…）

結論を先に書くと、

componetから電圧源を選択するとき、bvを使用して（voltageではだめ）、

V=(abs(int(time)-time+1/2)-0.25)*4

と書くと、三角波になる。

入力した式の意味

式を作る作戦

1. int(time)として階段関数を作る。

2. そこに-timeして鋸歯状波を作る（-側にオフセットがかかっている。）

3. +1/4して（オフセットを引いて）0を中心とした鋸歯状はを作る。
4. 絶対値をとると、三角波ができる（＋側にオフセットがかかっている。波高は0.25）
5. *4して波高を調節。

手順1~2を図にすると、

手順3~4を図にすると、

最後の2枚の図は、pythonで書いてみた。

一応コードを載せておく。

import numpy as np
import matplotlib.pyplot as plt


x=np.linspace(0,4,1001)

y1= x.astype(np.int32)
y2=x
y3=y1-y2


f,(ax1,ax2,ax3) = plt.subplots(3,1)

ax1.axhline(y=0, xmin=-10, xmax=10, c="k")
ax2.axhline(y=0, xmin=-10, xmax=10, c="k")
ax3.axhline(y=0, xmin=-10, xmax=10, c="k")

ax1.plot(x, y1)
ax2.plot(x, y2)
ax3.plot(x, y3)

ax1.set_ylim(-4,4)
ax2.set_ylim(-4,4)
ax3.set_ylim(-4,4)

plt.show()



f,((ax1,ax2),(ax3,ax4)) = plt.subplots(2,2)

y4=y3
y5=y3+0.5
y6=np.abs(y3)
y7=np.abs(y5)

ax1.axhline(y=0, xmin=-10, xmax=10, c="k")
ax2.axhline(y=0, xmin=-10, xmax=10, c="k")
ax3.axhline(y=0, xmin=-10, xmax=10, c="k")
ax4.axhline(y=0, xmin=-10, xmax=10, c="k")

ax1.plot(x, y4)
ax2.plot(x, y5)
ax3.plot(x, y6)
ax4.plot(x, y7)

ax1.set_ylim(-4,4)
ax2.set_ylim(-4,4)
ax3.set_ylim(-4,4)
ax4.set_ylim(-4,4)

plt.show()

LTSpiceを使って矩形波などを書く方法については、以前書いていた。（2016年！4年も前からマイナー記事を書き続けている。）

今年は新型コロナ一色の年だったが、私生活的には意外といい年だったと思う。

• Intel RealSense Viewerを使用

• pyrealsenseのインストール

の二つをやって手間取った。
2個目はいろいろやって、結局うまくいった。
原因はPythonのバージョン？

実行した順に書いていった。

Intel RealSense Viewerが起動しない。

そもそもIntel RealSense Viewerが起動しない。
下記エラーメッセージが出て（しかも全文表示されない）、はじまらない。

Backend in rs2_create_device(info_list:0000026F257C1900, index:0):
MFCreateDeviceSource(_device_attrs, &_source) returned: HResult 0x80070005: “アクセスが拒否されました。”

→win10のカメラのプライバシー設定が必要。
windows updateか何かのはずみで、Win10のカメラのプライバシー設定が入るらしい。
設定＞プライバシー＞カメラあたりで設定を変更してみる。

「Intel RealSense Viewer起動しない問題」はこれで解決した。

pip install pyrealsenseが通らない。

pip install pyrealsense

が通らない。結局通らなかったので諦めた。原因は、python3.8だったせい？

これを読んで下記条件を満たしたのに通らなかった。

• Cythonのインストール
• numpyのインストール
• pycparser のインストール
• wheelのインストール
• 環境変数PYRS_INCLUDES を追加。rs.hのあるフォルダに通す。
• 環境変数 PYRS_LIBS を設定。

結論

pip install pyrealsense

は通らない。

pythonのバージョンが悪いのかも…

上記設定はそのままで、
pipで自動でお任せしないで、.whlを自分でダウンロードしてpipした。
ついでにpythonのバージョンを3.7にした。

https://pypi.org/project/pyrealsense2/#files から、whlをダウンロード。現段階で、このサイトのリストにPython3.7までのファイルしかなかった。
自分はPython3.8を使っていた。
だからPython3.7にインストールしなおした

pip install （先ほどダウンロードした.whlのパス）

Windowsなら、ファイルブラウザ上でShift＋右クリックして「パスをコピー（A）」を押せば、クリップボードに絶対参照でパスがコピーされる。

結論

pyrealsenseがインストールできて解決。

まとめ

Python のバージョン（3.8）が悪かったのか、
pip install pyrealsenseじゃダメだったのか。
不明。

もしかすると、python3.7にするだけでpip install pyrealsenseは通ったのかも。

実験上手な奴は、条件を一つだけ変える。
実験下手な奴は、条件をたくさん変える。

僕は断然後者です。

このページの目的

Blenderを勉強する

３DCGを作りたくなった。

PS4のDetroit become humanをプレイして、リアルできれいなCGが作りたくなった。だからBlenderを勉強することにした。

ちょっと前はレンダリング方程式を勉強したくなったが、今は休止中。

欲求・目的・目標

欲求:全部できるようになりたい

これがあいまいだと、投げ出しそう。

• モデリング、レンダリング、ライティング、UVいろいろ分野があるが、全部勉強していきたい。

だから、

• 公式チュートリアル的なのをこなしたい。

  
  

目的:公式チュートリアルのトレース

YouTubeのBlenderの公式チュートリアルが英語なので、トレースし、日本語ブログ形式で解説する。

目標:記事10個

• ノック形式で10個こなす。

高速フーリエ変換より遅い。
でも、離散フーリエ変換より早い。
だから、「中速フーリエ変換」と名付けとく。

かかる時間
高速フーリエ変換<<<中速フーリエ変換<離散フーリエ変換
こんな感じ。

原理

離散フーリエ変換では、各周波数の三角関数と、入力信号の積の総和計算している。

三角関数の周期性を考慮して、入力信号f(x)を加法定理でまとめてから計算する。
→かけ算の数が減り早くなる。はず。

最後のほうに、exp(-ikn/N)をNを何個か振った状態で載せておくので、にらめっこして原理を察してほしい。
考え方は、途中まで高速フーリエ変換と一緒なはず。

低速フーリエ変換（普通の離散フーリエ変換）

import time
import numpy as np
arr = np.random.rand(1024)

start = time.time()

import math

N = len(arr)
W = [[math.e**(-1j*2*math.pi*k*i/float(N)) for i in range(N) ] for k in range(N) ]
out = [sum([arr[i]*W[k][i] for i in range(N)]) for k in range(N)]

elapsed_time = time.time() - start
print ("elapsed_time:{0}".format(elapsed_time) + "[sec]")

中速フーリエ変換

import time
import numpy as np
arr = np.random.rand(1024)

start = time.time()

import math

def order(Ns):
    if len(Ns) == 1:
        return Ns
    N1=Ns[1::2]
    N2=Ns[0::2]
    N2 = order(N2)
    return N1+N2

N = len(arr)
out = [0]*len(arr)#np.zeros(10)

N0 = range(N)
order1 = order(N0)

N1 = arr[:len(arr)/2]
N2 = arr[len(arr)/2:]
N3 = [N1[i]-N2[i] for i in range(len(N1))]
N4 = [N1[i]+N2[i] for i in range(len(N1))]
N5=[]
i1 = 0
i2 = 1

while i1 < len(arr):
    W = [math.e**(-1j*2*math.pi*order1[i1]*i/float(len(arr))) for i in range(len(arr)/2*i2)]
    N5.append(sum([N3[i]*W[i] for i in range(len(N3))]))
    if i1==len(arr)-2:
        N5.append(sum(N4))
        break
    if order1[i1] > order1[i1+1]:
        N1 = N4[:len(N4)/2]
        N2 = N4[len(N4)/2:]
        N3 = [N1[i]-N2[i] for i in range(len(N1))]
        N4 = [N1[i]+N2[i] for i in range(len(N1))]
        i2+=1
    
    i1+=1

for i in order1:
    out[order1[i]] = N5[i]

elapsed_time = time.time() - start
print ("elapsed_time:{0}".format(elapsed_time) + "[sec]")

高速フーリエ変換

import time
import numpy as np
arr = np.random.rand(1024)
start = time.time()

import math
def FFT(f):
    N = len(f)
    if N == 1:
        return f[0]
    f_even = f[0:N:2]
    f_odd = f[1:N:2]
    F_even = FFT(f_even)
    F_odd = FFT(f_odd)
    W_N = [math.e**(-1j * (2 * math.pi * i) / float(N)) for i in range(0, N // 2)]
    F = [0]*N
    try:
        F[0:N//2] = [F_even[j] + [W_N[i] * F_odd[i] for i in range(len(W_N))][j] for j in range(len(W_N))]
    except IndexError:
        F = [F_even + W_N[0] * F_odd]
    try:
        F[N//2:N] = [F_even[j] - [W_N[i] * F_odd[i] for i in range(len(W_N))][j] for j in range(len(W_N))]
    except  IndexError:
        F.append(F_even - W_N[0] * F_odd)
    return F

a1 = FFT(arr)

elapsed_time = time.time() - start
print ("elapsed_time:{0}".format(elapsed_time) + "[sec]")

N=2の場合のexp(-ikn/N)

左が実数で右が虚数。
青丸がサンプル点。

N=4の場合のexp(-ikn/N)

N=8の場合のexp(-ikn/N)

N=16の場合のexp(-ikn/N)

おれはいったい、何をやっているのだ。

Calculixとは

フリーの解析ソフト。

• ccx　ソルバー
• cgx　プリポスト

環境

• WSL

• Xming

• Calculix 2.11

• gmsh

  gmshの機能を使ってメッシングすることがある。

• imagimagick

  hcpyで画像作成するとき、convertするので使う。

• pip2

• matplotlib

• gnuplot

ノック11本目: exampleの実行 Contact/Hertz2D

データの入手

公式サイトからのリンクで行けるGithubのexample。

下記サイトのリンクから、自分のCalculixのバージョンにあったファイルを探す。

https://github.com/mkraska/CalculiX-Examples/releases

Contact

Hertz2D

平面ひずみ要素で、四角と1/4円をぶつけて接触させる。その時の応力分布をみる。
強制変位と荷重の2パターンで1/4円を移動させている。

ccx2.11では正しい結果がでない？

ccx2.11では収束しなかった。（頑張って条件をいじって結果をだしても、READMEのコンターの分布と同じにならなかった。）
ccx2.16のwindows版を使ったら、一発で通った。結果はREADMEの結果と同等のものが作れた。

READMEにあるmonitor.pyの結果と同じように収束しなかった。
このEXAMPLE全体に言えるが、初期設定では収束しないやつとかがある。
これもccxのバージョンのせいだと思っている。

内容

Hertzの接触を、FEMで求める。
球体とか円柱が、他の面等に接触した時の応力分布等を表したHertzの式というのがあるらしい。
周波数のHertz（Hz）と同じ人が考えたらしい。

まずはディレクトリの移動。

cd ~/CalculiX-Examples-2.11/CalculiX-Examples-2.11/Contact/

「~」は、サンプルを回答した場所のアドレス。

サンプルの実行

cgx -b pre.fbd

をすると、メッシュ、境界条件などのファイルが出力される。
これらのファイルは、.inpで読み込まれて使用する。

ccx Hertz

とすれば、同じ結果が出てくる…はずだった。

収束しない

STEP1は強制変位で動かしている。
STEP2は荷重変位で動かしている。
強制変位で動かしたほうが、収束はしやすいそうな。

STEP1は収束した。
ただ、STEP１はREADMEに書いてあるようなイテレーション回数で、収束しない。
STEP2の荷重で動かしているほうは収束しないで終わる。

なんで？

ccx2.11では結果が違う。

終息したSTEP1の結果だけを見てみる。

READMEにあるような、対称な応力分布になってくれない。片側の応力が出てこない。
接触面で応力が出ているけれども、サンプル画像では、内側で応力が発生していた。

ccx2.16(windows版)を使ってみる。

READMEと同じ結果が出てくる。
きっとこれが正しい結果なのだろう。同じ結果なので載せない。

ccxのバージョンは新しいのに限る。
以下の「余計なこと」は、すべてccx2.16でやった。

いろいろやってみる。

master slaveを入れ替えてみた

下記のようにマスタースレーブをひっくり返した。

*CONTACT PAIR, INTERACTION=contact, TYPE=SURFACE TO SURFACE
Ssblk,Sscyl

結果は変わらなかった。

平面応力要素でやった（デフォルトは平面ひずみ要素）

平面ひずみだから意味のある解析なのですが、あえて平面応力でやってみました。

ファイルall.mshの

*ELEMENT, TYPE=CPE8, ELSET=Eall

を

*ELEMENT, TYPE=CPS8, ELSET=Eall

に書き換えた。

cgx Windows版をメインで使わない理由

windowsで最新版(2.16)が使えるのに、わざわざlinux上の2.11を使う理由は、

• cgxのGUIのメニューが一部使用できない箇所があった
  （常にではなく、frdのロード時（？）など限られたタイミングで起こる。）。
  下記画像、（左）WSL（Xming）のcgx2.11の右クリックメニュー（右）windows版のcgx2.16の右クリックメニュー
  windows版ではdatasetsが出ない。載せていないが、矢印を展開した先でも似たような現象が散見された。

だけど、

• ccxは、windows版でも可。
  ccxはむしろ新しい方がいい。

Calculixとは

フリーの解析ソフト。

• ccx　ソルバー

• cgx　プリポスト

自分の環境

• WSL

• Xming

• Calculix 2.11

• imagimagick

データの入手

公式サイトのリンクから行けるGithubのページCalculixのexample。

下記サイトのリンクから、使用しているCalculixのバージョンに合ったファイルを探す。

https://github.com/mkraska/CalculiX-Examples/releases

ノック10本目: exampleの実行 Kasten

Kasten

Kastenは、ドイツ語で"箱"。

大変形、塑性、不整のある座屈（小さい内圧）、シェル要素とreadmeの書き出しにあった。Buckling with imperfectionは不整のある座屈と翻訳するらしい。

箱が座屈するモデル。 S8 (8-節点 四角形シェル要素)

 cd ~/CalculiX-Examples-2.11/CalculiX-Examples-2.11/Kasten

「~」は、サンプルを解凍したアドレス。

exampleには存在しないけれど

1/4モデルを作成して、

2. 材料非線形解析なしの場合

3. 非線形解析なしの場合。

を用意した。

1/4モデルの作成

1/4モデル作成のため、pre.fbdを下のように書き直した。

 pnt p1 0 0 0
 pnt p2 50 0 0
 pnt p3 0 50 0
 line l1 p1 p2 12
 line l2 p1 p3 12
 ​
 seta lower all
 swep all upper tra 0 0 150 20
 flip A002
 ​
 seta miryz l L004
 seta mirxz l L005
 ​
 plot la all
 plus pa all
 plus sa all
 rot y
 rot d 10
 rot r 10
 frame
 hcpy png
 #sys mv hcpy_1.png geo.png
 ​
 div all auto 3
 elty all qu8
 mesh all
 send all abq
 ​
 send all abq pres -0.01
 send lower abq nam
 send upper abq nam
 send mirxz abq spc 246
 send miryz abq spc 156

1/4モデルの為に、Kasten.inpの*boundary部分を下記のように書き直した。

 *boundary
 Nlower,1,6
 Nupper,1,6
 *include,input=mirxz_246.bou
 *include,input=miryz_156.bou

1/4モデルの為に、post.fbdをミラーコピーするように下記のように書き直した。

 copy all base1 mir x
 move base1 tra 100 0 0
 copy all base2 mir y
 move base2 tra 0 100 0

1/4モデルの解析結果

非線形解析なしの場合

非線形がない場合を見るために、Kasten.inpのStepのnlgeomを削った。

材料非線形がない場合

材料非線形がない場合を見るために、Kasten.inpの*plasticと、その下3行を削った。

あと、上記だけだと収束しなくなったので、*controlを追加して、収束条件を甘くしてしまった。

 *CONTROLS,PARAMETERS=FIELD
 5.0e-1,1e-1,1e-2,,2e-8,1.e-5,1.e-3,1.e-8

フリーの解析ソフト。

• ccx　ソルバー
• cgx　プリポスト

環境

• WSL

• Xming

• Calculix 2.11

• gmsh

• imagimagick

• pip2

• matplotlib

データの入手

公式サイトからのリンクで行けるGithubにあるCalculixのexample。

下記サイトのリンクから、自分のCalculixのバージョンにあったファイルを探す。

https://github.com/mkraska/CalculiX-Examples/releases

ノック9本目: exampleの実行 Linear/

Mesh1

アルミ製皿の固有値解析。
ポストのshapes.fblを使うと、複数の固有値のアニメーションを一気に作れる。
要素はIncompatible mode eight-node brick element (C3D8I)らしい。

ディレクトリの移動

cd ~/CalculiX-Examples-2.11/CalculiX-Examples-2.11/Linear/Mesh1

「~」は、サンプルを解凍した場所のアドレス。

実行

cgx -b pre.fbl
ccx modal
cgx -b shapes.fbl

pre.fblの内容

手打ちで作ってみたら、途中、GUIからframeを使わないと、全体が見えなくなった。

下記で円弧が出現してた。作成途中でできたD002、D003を通る、中心点D001の円弧の命令？（ 14）で書いたline。

line ! D002 D003 D001 14

1/4を作って、コピーで1/1にする。コピーのコマンドは下記。
最後のxは、yz面に対してミラーコピーの意味。yだったらx-z面

copy all new mir x
copy all new mir y

すべてコピーなので、不要な重複が出る。マージで消す。
マージは下記。

merg p all 0.001
merg l all 0.001
merg s all 0.001

all.mshの内容

先ほどのpre.fblの結果、all.mshが出力される。

all.mshには、節点座標とエレメントが使用する節点の情報が書いてある。
エレメントはIncompatible mode eight-node brick element (C3D8I)を使っていた。

modal.inpの内容

メッシュを読んでエレメントセットを用意して、

マテリアルを指定して、

上記マテリアルとエレメントセットをプロパティに設定して、

下記で、固有値計算になる。

*STEP*
frequency
～省略～
*END STEP

このmodal.inpが、ccxの解析で使う入力ファイルになる。
ccxで解析された結果は、modal.frdとmodal.datに出力される。

プロパティ設定の下の1は、3次元要素なので要らない気がする。あっても動くが。

shapes.fbl

出力した12個の固有値のアニメーションを自動で作るcgxスクリプト。

はじめに、wihleで使用する変数をvaluで用意している。

次に、readでmodal.frdを読み込んでいる。

次に、rotで角度を整えている。

最後に、whileで、すべてのデータセットのアニメーションを作成していた。
Whileについては以前やったので割愛。
whileはendwhileで終わるが、endwhileがテキストエディタの最終行に来たらエラーで止まるのは要注意。

cgxで最短でアニメを描く方法

多分、gifアニメは3行で出力できる。

アニメーションは、moviでフレームを設定したら、dsのaオプションで出力される。
多分imagimagickか何かが必要になると思う。出力されるファイルは、リピートされないgifアニメ。

read modal.frd new
movi frames 30
ds 1 a

Calculixとは

フリーの解析ソフト。

• ccx　ソルバー
• cgx　プリポスト

環境

• WSL

• Xming

• Calculix 2.11

• gmsh

• imagimagick

• pip2

• matplotlib

データの入手

公式サイトからのリンクで行けるGithubにあるCalculixのexample。

下記サイトのリンクから、自分のCalculixのバージョンにあったファイルを探す。

https://github.com/mkraska/CalculiX-Examples/releases

ディレクトリの移動。

cd ~/CalculiX-Examples-2.11/CalculiX-Examples-2.11/Contact/Shell1

「~」は、サンプルを解凍した場所のアドレス。

内容

下のような8節点シェルで作ったモデルの接触解析。
ハニカム（ハチの巣、ハニーコム）構造のサンドイッチの最小構成みたいな形の上に半球。
半球に対して、z方向に-2の強制変位をかけている。

初期設定でうまく流れなかったが、contact pairのmaster,slaveを入れ替えたらうまくいった。

Too many cutbacks

初期.inpのまま解析実行したら、下記エラーで止まった。

too slow convergence; the increment size is decreased to 3.125000e-02
the increment is reattempted

*ERROR: too many cutbacks
best solution and residuals are in the frd file

cutbacksの回数を増やしたら、止まらないんじゃないかと思い、cutbacksを増やす方法を検索。
controlをいじるといいらしい。 controlには引数が多すぎる。cutbackの回数以外は、初期設定にしておいた。
もしかすると、省略したら（何も書かずにコンマで区切ったら）よかったのかも。

下記*controlを*step(stepの外でも内でもいい？)に追加したら、このエラーは避けられたが…

** cutback 15
*CONTROLS, PARAMETERS=TIME INCREMENTATION
4,8,9,16,10,4,,15,,
0.25,0.5,0.75,0.85,,,1.5,

結局別のエラーが出た。
下記は結局出てきたエラー。

too slow convergence; the increment size is decreased to 7.629395e-06
the increment is reattempted

*ERROR: increment size smaller than minimum
best solution and residuals are in the frd file

これを解消する気にはならなかった。
*controlは消して、別の方法を探った。

contact pairがおかしかった

下記の2行目で、slave,masterを決めていた。

*contact pair, interaction=tool, type=surface to surface
Stopneg,Sind

dependent (slave), independent(master)の順番。

これを、下記のように逆さにしたら収束した。

*contact pair, interaction=tool, type=surface to surface
Sind, Stopneg

ちなみにsindは下記。

stopnegは下記。

結果

表示しているのはds 2 e 23で出せる、最悪主応力(worst Principal Stress)
READMEにある結果と最大値最小値一緒だし、同じ結果でしょ。多分。

フリーの解析ソフト。

• ccx　ソルバー
• cgx　プリポスト

環境

• WSL

• Xming

• Calculix 2.11

• gmsh

  gmshの機能を使ってメッシングすることがある。

• imagimagick

  hcpyで画像作成するとき、convertするので使う。

• pip2

• matplotlib

データの入手

公式サイトからのリンクで行けるGithubのexample。

下記サイトのリンクから、自分のCalculixのバージョンにあったファイルを探す。

https://github.com/mkraska/CalculiX-Examples/releases

ノック7本目: exampleの実行 Contact/Shellassembly

Contact

shell_assembly

4種類の.inpが入っていた。

• ペナルティ法でnode to surfaceの接触
• ペナルティ法でsurface to surfaceの接触
• *tieを使ったMPC接触
• *equationを使ったMPC接触

固有値計算してる。

他のサイトでは、「接触を考慮した固有値解析はできない」と書いてある。 )
不安になるが、まぁ、できるんならいいんじゃない…？
いや、良くないのかなぁ。「接触 固有値」でググると、皆、いろいろ工夫してるみたいだし。

どうやら、.inpファイルのstep, perturbationに鍵があるらしい。

READMEざっくり訳

*tieとペナルティ法を使った、MPC接触は、シェルのエッジとシェルの表面の間の、硬い接続（？）になります。 *tieを使った、Node to surfaceペナルティ法やMPC接触は似た結果になります。 ペナルティ接触の固有値の周波数は、ふつう、MPC接触を使った場合よりもきわめて低く、もっともらしいです。ペナルティ接触はcomplianceを足しているからです。

Surface to surface接触は、MPC接触とは全く異なった固有周波数になります。0付近に5つの周波数しかありません（6つの剛体モードがあるべき）。

（surface to surfaceもnode to surfaceのどちらも）ペナルティ接触は、モード解析が必要です。モード解析は、静的なstepに先立って摂動ステップとして行われます。 摂動解析の最初の結果のインクリメントはモードシェイプではないことに気を付けてください。

*equationを使ったMPC接触は、シェルとフェイス接触を、ヒンジ接続にします。（7剛体モード）

う～ん分からない。訳が合ってるのかも不安。

まずはディレクトリの移動。

cd ~/CalculiX-Examples-2.11/CalculiX-Examples-2.11/Contact/

「~」は、サンプルを解凍した場所のアドレス。

最新版だとShell0に入っているサンプル。

操作

cgx -b pre.fbd
ccx tie
cgx tie.frd

こんな感じでtieの部分をequ、pc-ss、pc-nsに替えて実行する。

pre.fbd

seto～setcを使ってセットを作っていた。

上側のプレート1をswepで作って、

下側のプレート2をswepで作って、

qu4（inpでいうs4（4節点シェルエレメント））を指定してメッシュ作成して(all.mshにS4と記入される)、

プレート1のサーフェイスを出力して(S1.sur)、

プレート2のエッジ（L007 ）を出力して(k1.sur)、

すべてのメッシュを出力してして(all.msh)、

プレート1のサーフェイスとプレート2のエッジ(L007)を123456のMPCとして出力していた(k1.equ)。

PgDn, PgUpでplotするsetの切り替えができた！

PgDn、PgUpでセット間を行き来できるのは知らなかった。

plot fa all

と打った後に、PgDn、PgUpすると、作成した各セットを行き来できるので、いちいちplotコマンドを打たなくていい。作ってる最中に便利かも。

.inpを見る

tie.inp

*tie,name=t1,position tolerance=0.1
Sk1,Ss1

とあった

pc-ss.inp

stepの部分

*step,nlgeom
*static
*end step
*STEP,perturbation
*frequency
~~~
~~~
*END STEP

と書いてあった。
「PERTUBATIONってなんだ？」となったので、取説を検索すると、

The parameter PERTURBATION is allowed for FREQUENCY and BUCKLE steps only. If it is specified, the last *STATIC step is taken as reference state and used to calculate the stiffness matrix.

パラメーターPERTURBATION は*FREQUENCY(固有値解析) *BUCKLE(座屈解析)にだけ使用できます。もし、これが指定されると、最後の *STATIC stepの剛性マトリックスが参照され、計算に使用されます。

たしかに、直前に*staticしてるstepがある。その結果を読み込んで、固有値計算をしているらしい。

接触の部分

*surface interaction, name=Klebung
*surface behavior, pressure-overclosure=tied
1.E7
*friction
1,1.E7
*contact pair, interaction=Klebung, type=surface to surface, adjust=1
Sk1,Ss1

pc-ss.inpとは、pressure-overclosure=tiedと、type=surface to surfaceが違う。

*frictionの下の数値が、摩擦係数と粘性slopeらしい。取説には下の様に書いてある。

For face-to-face penalty contact with PRESSURE-OVERCLOSURE=TIED the value of the friction coefficient is irrelevant.

PRESSURE-OVERCLOSURE=TIEDを使ったface to faceペナルティ接触にとっては、摩擦係数は無関係です。

とあるので、1っていう数字は要らないんじゃないかと思ったが、なきゃ無いでエラーが起こるのでそのままにした。

pc-ss.inp

接触の部分

*surface interaction, name=Klebung
*surface behavior, pressure-overclosure=linear
1.E7,10000,10000
*friction
1,1.E7
*contact pair, interaction=Klebung, type=node to surface, adjust=1
Sk1,Ss1

pc-ns.inpとは、pressure-overclosure=linearと、type=node to surfaceが違う。

equ.inp

これは、下記で読み込んでいる、k1.equにの条件が書かれているところに特徴があった。

*include, input=k1.equ

*equation を使って、MPCの式が書かれていた。

*EQUATION4
4
6,1,-1.000000000000
2,1,0.800000000000
11,1,0.200000000000
12,1,0.000000000000

4というのが、式の項の数で、
その下に、ノード番号、拘束方向、変数の係数が書かれている。変数は、ノードの拘束方向の変位。

k1.equを見ると、近い4つのノード同士で式を作っていた。ノードがどうやって選ばれたのかは知らない。

結果

固有値の結果は、.datファイルに出力される。

cgxのwindowの「Time:」って書いてある横の数値は、そのまま周波数として読めそう（CYCLES/TIME）。

*el fileS

PgDn, PgUpでplotするdatasetの切り替えができた！

plotコマンドのところでも書いたが、dsコマンドでも、PgDn, PgUpデータセットの切り替えができる。
下記コマンド等、適当に入力後、PgDn, PgUpで、データセットを切り替えられる。

ds 1 e 1

問題点

cgxのオプション

cgxを使って結果の読取りをする際、cgxの起動時に -bオプションををつけると結果が表示されない。何もつけない（自動）か、-readで出てくると思う。

cgxのdataset

cgxのGUIからdatasetを選択すると、アプリケーションが勝手に終了して下記の表示が出た。

freeglut (cgx): Menu manipulation not allowed while menus in use.

表示させるためには、コマンドから下記の様に入力。

ds 1 e 4

データセット1の4番目の結果を表示。何番に何が入っているかは、出力要求によって違う。

windows版ならGUIは問題なく動くと思う。