MC 日記
TeX メモ書き DD
読む人いないかもしれんけど研究日誌。 My Works Log 2003 2004 2005 2006
Current interest: dislocations, random field and random anisotropy, spin glass, polymer, chiral and frustrated spin models, morphology of fracture and cracks, Monte Carlo methods
★研究キュー: Fracture model [diffusion + cohesive zone + FEM] (active, Phase 2), Dislocation wall formation in 2D (Phase 3), Dislocation Dynamics (active, Phase 1), FSS of SAW based on phi^4 theory(susp), PERM for Self Avoiding Membrane(susp), Random Field XY model(susp),
Phase 1: 文献調査、モデル構築、予備計算
Phase 2: モデル決定、本計算、データ収集
Phase 3: データ整理、論文執筆
Acronyms: RAM=Random Anisotropy Magnet, RF=Random Field, MC=Monte Carlo, RG=Renormalization Group, FP=Fixed Point, SG=Spin Glass, RSB=Replica Symmetry Breaking, FSS=Finite Size Scaling, LRO=Long Range Order, QLRO=Quasi Long Range Order, SAW=Self-Avoiding Walk, AF=Antiferromagnetic, WF=Wilson-Fisher, IASCC=Irradiation Assisted Stress Corrosion Cracking, FEM=Finite Element Method
RR=Referee report, RRC=Reply to Referees' Comments

TODO list
CZ論文

2009/11/09
ESでのVASP鉄螺旋転位芯計算、9x10x2原子でなんとかいけるメドがたった。 プロダクトラン開始。 外側は固定して無理矢理周期境界。境界はぐちゃぐちゃになってるけどまあ無視。 Hのbinding energyとmigration enrgyを出す。 しかしやはり余分なXY方向の力が多少出てしまう。0.4 eV/Angstrom程度。


亀裂計算はいいかげんMD切り上げてFEMに移行。メッシュ切り直し中。
2009/11/02
Symmetric tilt GB でtensileが[100]に近い場合に破面が ギザギザとなりfracture toughnessが上がり、 水素がある場合はそれが抑えられることを確認。ようやくメドがついた。 しかし厚みがある場合は局所的なギザギザは抑えられる。しかし長い距離を 進んでゆっくり緩和すればたまたま集団で蛇行することもあろう。 高温でやるべきか。厚さと温度とスピードをきちんとスケーリングに乗せないと いけないだろうな。
2009/10/26
Kの推定値の統計誤差がひどい件、結局<y/x>なんて推定量がセンスなさすぎた。 線形回帰なら<xy>/<xx>が正しい。
2009/10/23
MDで fracture toughness = stress intensiti facotr の臨界値を出す 方法はだいたい固まった。周囲を固定し原子全体を亀裂周囲の弾性解で 少しずつ変位させていく。K値の増分を2e4 (Pa m^1/2) 程度にすれば 亀裂周辺でも変位は最大で0.01オングストローム程度なんで変な欠陥が 導入されることはない。以下はそうやってひっぱった後の配位

色は初期配置に基づいて分けている。弾性解もいっしょに線で示してある。 亀裂から離れたところでは一致している。 亀裂が進展したことによって周辺では進展しない場合の弾性解からずれる。 したがって全体にかけたK値は使えない。 そこで亀裂周辺の原子の変位を、亀裂先端位置とK値をパラメータとした 弾性解でフィッティングしてローカルなK値を計算した。応力方向の変位 をuy, K=1の時の弾性解の変位をuy0 とし、log(uy/uy0)のばらつきが最小となる 亀裂先端位置を求める。このとき先端位置周辺の原子だけ用いる。 近すぎる原子は弾性解からずれるので除外。具体的には1nm〜2nmの 範囲とした。次に uy/uy0 の平均値からK値が決まる。

それで計算できたが、統計誤差が大きすぎる。数百原子あるはずなんで average out するはずだが、と思ったが、考えてみると亀裂近くの 変位はひじょうに小さく、熱振動にうもれてしまう。 なのでLAMMPSで位置の時間平均をとるようにした。
#10 step おきに観測しそれを300個平均することを4000ステップおきに行う
# すなわちはじめの1000ステップを捨てる
fix posx all ave/atom 10 300 4000 x 
fix posy all ave/atom 10 300 4000 y 
fix posz all ave/atom 10 300 4000 z 
#ファイル pave.step数 に出力
dump posdump all custom 4000 pave.* id type f_posx f_posy f_posz

2009/09/17
Pasadena方面に現状報告。
Currently I'm compiling Fe+H grain boundary cohesive element data with pre-crack from MD simulations and observing critical K1c. Unfortunately K1c does not drop significantly with segregation of H on the grain boundary, contrary to our early expectation, and I'm experimenting various configuration in which H reduces K1c.
動画テスト
あーこまったのう。 別シナリオとしてHがFeの応力誘起BCC->FCC転移を抑制するかと思ってしらべたが、あまり効果無し。

新ESでの転位VASP計算も150原子でLREAL=TRUEでも時間内に電子密度が収束しない。 どうしたものか。
2009/09/03
やはり普通に亀裂周囲の等方弾性体の解を使って変位させ、 限界の応力拡大係数Kcを求めるのがスマートかつ破壊力学屋にも 受け入れやすそうだな。すべてK値で議論できる。 ただし亀裂周囲では激しく変位するのでコントロールは系の外周のみで 行い、これを変位させ固定して緩和するのがよさげだ。ただしゆっくり やらないと固定と内部領域がはがれてしまう。

始め考えていた、界面の開きと水素濃度の熱力学は実際MDやってみると 机上の空論だと分かってきた。亀裂が進展するかどうかは 先端部分での atomic, subatomic なスケールの振舞いで決まる。 この影響が拡大して亀裂がマクロに成長するのは量子力学の観測のような 感じだ。亀裂は一本しかないので多くの可能性のある亀裂について 平均した量というのはあまり意味がない。
マクロな量である、界面の開きと水素濃度といったものに 対応するミクロな状態は、その間でエルゴード性があるわけじゃなくて 固まっているので、マクロな量で統一的に議論するということができない。 ただ、亀裂先端は三次元だと横にひろがっている。 どこか一番弱いところで亀裂が進展し、残りが応力集中でひきずられるという いみで多少のアベレージングはある。 しかし割れ遅れた部分への応力集中が十分でないと亀裂先端は直線から どんどんずれる。

そういった困難な状況で工学屋が唯一拠り所にできるマクロな量がK値 ということだ。ただしこれも万能ではない。
2009/09/01
LAMMPS でトータル応力の計算メモ
compute mystr all stress/atom
compute sxx all reduce sum c_mystr[1]
compute syy all reduce sum c_mystr[2]
compute szz all reduce sum c_mystr[3]
#1000 step おきに出力、10ステップおきに計測し100回の平均を出力
fix fixid1 all ave/time 10 100 1000 c_sxx c_syy c_szz file str-log
実際の応力は出力を系の体積で割らねばならない。units metalだと 1e5倍すれば Pa になる。

亀裂先端が系のとなりの端まで来ている場合の限界応力をMDで 計算することを考える。系はL*Lを2つ張り合わせて上下にひっぱる。 亀裂周辺の弾性解は上下の辺上において横方向の変位はほぼ一定、 縦方向は亀裂真上とその反対の角ではちょうど二倍の変位になる。 実際有限要素計算でとなりの要素が割れていてバネが切れているなら、 亀裂真上のノードと真下のノードの間にエフェクティブに働く バネ定数は半分になっているので、変位としては反対側の二倍に なっていると期待できる。そのような条件で小さい初期亀裂を 入れてMDをやって限界応力orひずみを もとめれば亀裂進展条件がわかる。

有限要素計算で変位が二倍からずれていたら、、どうしたもんか。
また横方向は open boundary になる。
2009/08/18
memo:
Langevin Eq.
ma= -Ux - mv/d +sqrt(kT m / dt d) N
d: dump factor time
N: Noise
U: Potential


vnew = -dt/m Ux +(1-dt/d) vold +sqrt(dt kT/md) N
dt=dとすればMC的 Langevin
2009/08/11
LAMMPS intel cc with math kernel lib Makefile
--- Makefile.intel_mkl.orig	2009-08-11 16:39:53.000000000 +0900
+++ Makefile.intel_mkl	2009-08-11 18:37:27.000000000 +0900
@@ -14,13 +14,13 @@
 
 include		Makefile.package
 
-CC =	        mpiicc	
-CCFLAGS =	$(PKGINC) -O3 -fno-alias -ip -unroll0 -g -DMPICH_IGNORE_CXX_SEEK -DLAMMPS_GZIP -DFFT_FFTW -I/opt/intel/mkl/10.0.011/include/fftw 
+CC =	        /opt/intel/impi/3.2.1.009/bin64/mpiicc	
+CCFLAGS =	$(PKGINC) -O3 -fno-alias -ip -unroll0 -g -DMPICH_IGNORE_CXX_SEEK -DLAMMPS_GZIP -DFFT_FFTW -I/opt/intel/mkl/10.1.3.027/include64/fftw 
 DEPFLAGS =	-M
-LINK =		mpiicc
+LINK =		/opt/intel/impi/3.2.1.009/bin64/mpiicc
 LINKFORT =	
-LINKFLAGS =	$(PKGPATH) $(LINKFORT) -mt_mpi -L/opt/intel/mkl/10.0.011/lib/em64t 
-USRLIB =	$(PKGLIB) /opt/intel/mkl/10.0.011/lib/em64t/libfftw2xc_intel.a
+LINKFLAGS =	$(PKGPATH) $(LINKFORT) -mt_mpi -L/opt/intel/mkl/10.1.3.027/lib/em64t 
+USRLIB =	$(PKGLIB) /opt/intel/mkl/10.1.3.027/lib/em64t/libfftw2xc_intel.a
 FORTLIB =	
 SYSLIB =        $(FORTLIB) -lstdc++ -lpthread -lmkl_em64t -lguide 
 ARCHIVE =	ar

2009/08/05
Σ11 symmetric tilt GB [110]x[113]x[332]をひっぱってSSカーブが伸びるのはFCCに変態してるせいだった。 うわちゃ、また面倒な。BCCを100方向にひっぱるとあるところでFCCになる。 他の界面でもtensileとshearが混ざって[100]方向になればその可能性が出てくる。 まあ破壊起点の推定とかを考えるとそうでないものに集中すればよいのかもしれん。
水素をいれてもあまり変わらんな。結論:Σ11は割れない
2009/07/22
μがd と共に増える、μ_d >0 というのは界面が開けば水素が いごこちよくなることを意味する。 -F_c,d>0なんでF_d,c つまり応力の濃度微分<0となる。 つまり水素が入れば応力が下がってもっと開く。開けばもっと水素が入れる。
と思ったが、応力一定の条件でcが増えた場合のμの増分を計算すると、 Fdc*Fcd という項として表れるのでこの符号は無関係という面白いことになった。 実際にC-Fd面でFcのコンターを書いてみないと分からない。
しかし、粒界面が途中まで割れて途中からひっついてるなんて状況は ありえない気がしてきた。
とりあえずΣ3、11、19それぞれ引っ張ってSSカーブ書いてみたが、 かなり様子が違って面白い。何度が不連続に変わる。ヒステリシス をFEMモデルに採り入れるべきか。そしてギザギザでしばらく伸びる。 応力で一意に状態を指定できない。
2009/07/21
鉄の粒界面+水素の系の応力に対する非線型応答をまとめてFEMで使う。 非線型といっても途中でバネが切れる形になる。両端のずれと水素濃度を変数とした 自由エネルギーF(d,c)をVASP計算から出せば全ての計算ができる。 応力は Fのd微分、水素の化学ポテンシャルはFのc微分。 応力のc微分と学ポテンシャルのd微分は等しいはずだが、MD計算でやると かならずしもそうはいかない。引っ張ってから水素を粒界面におくのと、 水素をおいてから引っ張るのの違い。
またShear応力への応答がめんどくさい。スライドしたり粒界面が移動したり。スライドしたあとも摩擦があるし。 Σ値によっても方向によっても全く違うし。
2009/01/28
一部でメッシュ細かくし形状関数を重ねる手法
J.Fish, Computers and Structures Vol.43, 539 (1992)

これは便利。 これを粒界亀裂先端で使えるようにせねば。

基本的にメッシュ内部に追加自由度をおいて、メッシュ境界では0になるshape function を使ってenrichしてやるのか。progressiveに細かいのをどんどん足していけばいい。 shape funcの空間微分の積の積分を計算する必要が有る。 四角限定なら(1-|x|)(1-|y|)という形か。場合分けが面倒かも。 三次元で三角柱の場合どうすんだ。難しい。
2009/01/05
いろいろdisturbancesによりメインの課題が進まず。
発注したコード用の3Dメッシュで0体積四面体が多数あると報告をうける。 今からメッシュ生成コード作りなおしてたらものすごい手間だ。でもどうもデータの 1-origin, 0-origin のミスのようだと判明した。そうであってくれ!たのむぜ。 将来的にはちゃんとしたメッシュ生成コードは必要なのだが。pet project として spare time に完成させるつもり。
とりあえずメッシュは手を離れた。次のデバグだ。
問題:あるポテンシャルを使って緩和された個体状態の 金属原子の任意の配置があたえられている。金属ー水素ポテンシャルもある。 水素のエネルギー極小位置をすべて列挙せよ

全般的にDelaunay分割した四面体の中心にあるが転位や表面近くなどは違う。
表面などで幽霊原子を導入しDelaunay分割でbulkの四面体構造が保たれるようにする。
ローカルにMDしてエネルギーを計算する。これが幽霊がある場合がなんかおかしい。
2008/11/10
CZ要素論文構成
環境助長亀裂進展では粒界面に脆化元素が侵入し表面エネルギーを低下させ 亀裂進展応力を低下させる。 この効果はDFTなどの量子計算で評価することが出来、 それを有限要素計算に取り入れることで実際の臨界応力を見積もることが出来る。 環境助長亀裂進展の典型的な例として水素脆化と応力腐食割れがある。 どちらも水環境に接する鉄あるいは鋼で発生し、脆化元素は水素と酸素である。 本件級では鉄の粒界面に水素および酸素を配置したDFT計算および熱力学的考察により これら元素の偏積エネルギーおよび脆化効果を統一的に扱える粒界面のモデルを 定時する。 さらにマクロスケールでの計算を可能にするためいわゆる繰り込まれたCZ力を 構築する。これは亀裂先端でのみ細分化されたメッシュを用いその自由度を 最適化することで、与えられた粗自由度に対する亀裂先端での詳細な情報を 繰り込んだエネルギーお呼び力を計算するものである。 簡単な二次元の系についてテスト計算を行い、この面間要素が正しい結果を出すことを示す。
考察:SCCは水素か酸素か?水素環境ではKSCCはあまり下がらない。でも粒界破壊は起こる。これは水素脆化による ものと思われる。酸素は1/10まで表面エネルギーを低下させるのでこちらの方が効果が大きい。
2008/10/23
動画完成
No.1 No.2

2008/10/22
水素の拡散に亀裂が追随する計算もできた。 ある程度水素の力を借りてすすんだ後は亀裂が長くなって 応力集中も√長さで増加するので、その後は瞬間的に割れる。 ここまでの時間が遅れ破壊の破断時間であろう。 もとの亀裂の2倍になるまでの時間を計測、応力に対して プロットすれば実験と合うはずだ。
結局亀裂先端でÅ単位のメッシュを切った。原子より細かいメッシュ意味あんのか、て 感じだが粗いとFEMが収束しないのでしょうがない。 ひずみエネルギーと面間引力の合計を変分で最小化すると、開き具合をδ(x)、引力をF(δ)とおけば δ(x)は d^2δ/dx^2 = K F(δ) という方程式の解になる。

結局全体サイズは100nm×100nmがいまのところ限界だった。もっとメッシュを工夫すれば1μmくらいはいくけど時間がない。
二次元でこのサイズならMDで楽勝じゃん、というのは言わない約束でしょおとっつぁん。すまないなあゲホゲホ
2008/10/20
水素なしの場合&水素ありで拡散なし、計算終了。
実験では水素がフルに入ると粒界面の引力が1/3まで低下する。 VASP計算では1/4まで下がっているが、これをそのまま使用。 McLean の式によりフルに入ることはなく80%くらいで止まる。こととき1/3程度になって実験と合う。
粒界面<>バルク間の水素のやりとりも、原子スケールに一旦落として考えて 時間発展式を構築。うまくいっている模様。
原子一個の単位時間あたりジャンプ確率をJとする。 拡散定数とは D=J a0^2/6 の関係。a0は一回のジャンプの距離。
時間 dt の間に面積A0^2の粒界面(原子一個分)にバルクからくる原子の個数はCb(1-Cgb)(dt J /6) ここでCbはバルクサイトの占有確率、CgbはGB面での占有確率。すでに人がいた場合は入れないのでファクター(1-Cgb)がある。
逆にでていく原子の個数は Cgb exp(-Es/kT) (dt J/6)、Esは原子一個あたりGB偏積エネルギー。これが等しいとおけば McLean の式が出る。 flux が大きすぎてどっちかが空になったり飽和したりする場合は fluxを小さくする。 粒界面に接するバルクの体積要素はメッシュを細かくする必要がある。現在10nmでやってる。
割れた面に水素が吸着していくので、新たに割れようとする面にはGB拡散では水素はいかない。周囲のバルクからかき集めてくるしかない。
2008/10/17
Back of the envelope calculation でバルクの拡散はいらないことが判明。
まず標準的な鉄バルク中の水素濃度。溶存がendothermicで +0.3eV 上がる。 300K でのBoltzman factor は 1e-5 程度。 格子定数は a0=2.85Å、1 unit に水素は12個入る。 一方、表面では 111面だと sqrt(12) a0^2 の面積に最大24個入る。 結晶粒の大きさは100μm程度。この中の濃度1e-5の全バルク溶存水素が粒界面に析出 したとしよう。そうすると界面での occupation ratioは ちょうど1.0程度。
あ、ちょうどいいか。でもマルテンサイト鋼では粒はもっと小さい。 界面が開いて飽和濃度が上昇した時にもっと水素が入るわけだけど、 バルクからの供給だけでは足りない。どちらかというと 面濃度は変化せず周囲のバルク濃度が下がることによって McLean Eq. が満たされることになる。
界面内部での拡散を考えれば面濃度自体が変化できる。それにより 面間力が低下する。

さて粒界面の拡散でやると決めたら離散化スキームを決めねば。
		   /
                  V 
V=======V=======V
                  V
                   \
図のVと書いたFEM頂点に定義するか、その間の===のエレメント中央に定義するか。 通常の線形形状関数を使った定義では保存則を守るのが面倒なので、FVMの定式化により Vか=で定義しその周辺の平均濃度を表す。▽C 計算は濃度差/距離で。三十点をはさむ場合は Vで定義しないと面倒。三十点では常にμが等しくなるよう瞬時に調節する。 fluxを計算した後はそれを担当領域の体積で割って増分とする。
2008/10/16
Bulk の拡散は完成。しかしまだメッシュが粗いっぽい。 1%ほどのstrain で水素一個の溶存エネルギーは0.2eVほど下がる。 300Kだとボルツマン因子はすごく大きい。したがって そのくらいstrain のあるあたりは濃度が数桁大きくなる。

粒界面エネルギーもそれらしきものを作った。 A(C) + B(C) LJ(x) という形でもできた。B(C)は二次式。 C微分が-になる、つまりchemical potential μが得な領域を x-C 平面で 表示し、x=0ではC=0.5まで、 x>>1 ではC=0.8程度まで入るようにした。 さらにB(C)は C=1のときC=0の場合の1/4になるようにした。 μ(x=0,C=0)=Segregation energy, μ(x>>1,C=0)=Absorption energy ともにreasonableな値となる。
時間スケールは、応力緩和、McLeanの式に従うバルク<>界面の水素移動、 バルク拡散、粒界拡散、表面拡散などがある。MDによると 粒界拡散はバルクと同じくらいなので今回は無視。表面拡散は速いが 今回は水素供給源がバルクなので関係ない。 応力緩和は無限に速いとする。McLeanの緩和はバルク拡散よりかなり速いはず。 これだけやればバルク拡散はいらないかも、という気がしてきた。
2008/10/15
全体1μm、亀裂0.1μmの系ですべてやることに決定。 FEMの精度はÅが必要。スケール差10^4はまあ許容範囲だ。 結晶粒は実際は50μmくらいあるんだが、まあよかろう。
4-5GPa で亀裂が進展した。実験では2GPa くらい。これは初期亀裂の 長さの-1/2乗に比例する。現実には1μmほどの析出物が破断し初期亀裂 となる、と考えればオーダーはあっている。

あとは粒界面エネルギー Egb(x, C)を決めねば。開き x と水素濃度C の関数。 これを x で微分すると面間引力、Cで微分すると chemical potential となる。 うまいこと決めてVASP計算と合わせる。x 依存性は A(C) + B(C) LJ(x-x0(C)) てな感じでLJポテンシャルをスケールし移動する。係数をうまく決めて C微分がVASPの結果に合うようにする。

とりあえず今日は拡散を解くコードを追加せにゃ帰れん。
拡散は応力mesh と Dual な meshで計算する。


2008/10/14

くっそー時間足りない。線形バルク弾性+LJ粒界面エネルギーのBFGSによる最小化および SX-6でのベクトル化、初期亀裂導入と亀裂先端でのメッシュ細分化、応力の可視化、まで終了。 まず水素入れずに全体応力5GPa固定、全体のサイズをいろいろ変えて割れるところを探す。

2008/09/30
Segregation Energy が濃度に依存する場合のMcLean の式。 http://en.wikipedia.org/wiki/Segregation_in_materials Bulk サイトN, Bulk原子P, GBサイトn, GB原子p の時エントロピーは kT Log ( NCP nCp)
Stirlingの式 Log(n!)~ 1/2 Log(2πn) + n(Log n -1) ~ n(Log n -1)
Log(xCy) ~ xLogx-yLogy -(x-y)Log(x-y)
Segregation energy E(C) はGB占有率Cの場合のトータルenergy。 E'(C)が占有率Cの時もう一個入る時のエネルギー変化。
最終的に C(GB)/(1-C(GB)) = C(B)/(1+C(B)) exp( E'(C(GB))/kT)
C(B)<<1であれば、C(GB)= C(B) exp(..)/(1+C(B) exp(..))
バルクの拡散を解いてC(B)を計算、粒界面においてC(B)をfixし C(GB)をニュートン法で決める。
ESのVASP計算は最後にCHGCARを出して終わったけど qsubスクリプトにそれをとってくる指定をしないといけなかったらしく 消えてしまった。痛い。結局転位は動かなかった。
プリンストンの同僚がついにFe-Hポテンシャル完成といってきた。 こちらで様々な欠陥とのbinding energyを計算せねば。
三次元 mesh はいちおう切れた。薄い四面体とかあるが、これはしょうがない。
通常の場合なら、必要に応じて表面や辺にもノードを追加していく。しかし 今の場合、ある結晶粒の多面体でそれをやると既にメッシュを生成した、それに 接する多面体でもデローニ分割をやりなおしする必要があり面倒。何回かiteration して収束させればいいんだろうけど。デローニ分割の途中経過を保存しておき、 面にノードを追加した場合は相手の結晶粒でもそのノードを追加する。計算量はO(1)。 幸い各デローニ分割処理ごとに入力や途中経過をひとまとめにしたクラスを作ってある。 あー、しかしstatic なメンバもあったな。このままじゃ使えない。
別の方法としては結晶粒を張り合わせる平面要素で、両側のノードが対応してない ような場合でも補間して計算するという手もある。 まあしかし前処理で全て解決できるのならそうしたい。
フロリダのMMMでは二次元の応力+拡散連成をやる。メッシュはすでにいいのが出来た。 MD計算によれば鉄中水素は粒界面拡散はバルクにくらべてそんなに速くない。 なのでB拡散+GB拡散+McLeanの式によるB-GB接続、でやる。 応力も計算。B拡散は静水圧成分p に影響され flux= D▽C(B) +C(B) VH kT ▽p 。
C(B)を応力計算のメッシュのノードに置くのか三角形の中心に置くのか、決めねばならない。 pは各三角形で計算されるので三角形に濃度を置くと楽。しかし中心というのを 何にするか。重心より外接円中心がいいが、かならずしも内部にあるとは限らない。 重心でも別にいいか。 メッシュ生成の時にそういう条件にしてやればいいが。▽C(B)は隣接する三角形の濃度差/三角形の中心の距離。 外接中心なら中心ー中心の線が境界に垂直なのでfluxの和が計算しやすい。 C(B)▽p は、▽(C(B)p)とするか。あるいはC(B)は二点の平均にするか。
応力によりGBが開くとMcLeanの式のポテンシャルが変化する。C(B)も変化。 バルクの弾性はC(B)には影響されないが、GBの接着力はC(GB)により低下する。
ああああ。できるかな1月で。B拡散+(GB拡散:しない)+McLean+B応力+GB接着力の連成だ。
応力でGB開く>水素がGBに流入>もっとGB開く>もっと水素流入>>>>破断。
水素拡散が律速となるが、鋼の場合の速さはよく分からん。純鉄単結晶ならデータはあるが。
オーダーが1日であれば実験とあう。
2008/09/24
ES 鉄転位計算、前年度のデータをよく見直す。 転位が動くぎりぎりのところではかなり原子の緩和ステップを100とか取ってる。 それでもstress は実は減り続けてたのを途中で打ち切ってた。がびーん。 人の計算なのでよく把握してなかった。今年のデータを見ると、 動きそうなstrain でけっこうstressが減ってる途中でやめていた。 なので続行決定。あと2回くらいしか時間がないが、収束した後は早いんで なんとかなるだろう。
ノートPCは修理を依頼しに電話したとき立ち上げたらWINDOWS起動までいったんで おや、と思ったが、今日修理にきてもらい、はやりHDがイカれてた。出張に持っていく つもりなんでプレゼン直前で落ちたりしたらかなわんのでHD入れ換え。
メッシュ生成は三次元までいちおうできた。しかし平らな四面体ができる 問題がやはり出た。メッシュ生成のレビュー論文で三次元特有の問題としてあげて あったので把握はしていたが。
Jonathan Richard Shewchuk, www.cs.berkeley.edu/~jrs/
2008/09/19
二分木クラスでコア吐いたんでこいつを整理してて終わった。 だいぶすっきりした。
#ifndef BTREEH

#define BTREEH

#include 
#include 
#include 

// Manage list of pointers which is associated with several integers
class Btree
{
   public:
       // Find a pointer associated with a key of integers, using binary tree search.
       // Returns NULL if not found
       void *Find(int *key);
       // insert a pointer associated with a key of integers into the list
       void Insert(int *key, void *iptr);
       // delete a pointer associated with a key of integers into the list
       void Delete(int *key);
       void DeleteAll();
#ifdef ALLOW_STR
       // use string as a key instead of integer
       void *Find(char *keystr);
       void Insert(char *keystr, void *iptr);
       void Delete(char *keystr);
#endif
       // total number of pointers
       int  NData();
       // get a complete list of pointers
       void GetAll(void **ptr);
       // init list class. 
       Btree(int nk, bool do_sort0);// nk=number of key integer (0 if key is string) 
                                    // order of int is irrelevant if do_sort0=true
       ~Btree();

       void Loop1();
       bool Loop2();
       void Loop3();
       void *LoopVar();
#define BTFOR(x) for((x)->Loop1();(x)->Loop2();(x)->Loop3())
       // example:
       // Btree bt(2, false);
       // BTFOR(&bt)
       // {
       //     someclass *p = (someclass *)(bt.LoopVar());
       // }

   private:
       int nkey;
       // root of bi-tree, head of simple pointer cascade, current loop var
       class BtNode *Root, *Head, *loop, *loopnext;
       //returns  the parent's pointer which points to the hit
       class BtNode **FindAux(int *key);
#ifdef ALLOW_STR
       class BtNode **FindAux(char *key);
#endif
       void DeleteNode(BtNode **p);
       int nnodes;
       bool do_sort;
};

#endif

2008/09/18
MACへのデータ移行でいろいろ時間を取られる。 StuffItがアホでディレクトリ構造を解凍しないとか。 もらった論文はメールの添付のままSylpheedのデータで保存して Sylpheedで検索とかしていて便利だったが、OSXで入れようとして 挫折。標準メーラーにimportする方法もググって見つけたがうまくいかず。
メッシュ切るためのクラスいろいろを、書き直してたりして 機能的には進歩ないがコードはかなり短くなりバグも入りにくくしている。 この時間の投資が後で回収されるかどうか。

2008/09/17
雑用WINマシンは見捨ててMACに乗り換えた。なんだパワポもMS-*も入ってるじゃん。 FuguというSFTPツールを入れて環境が95%復活。快適じゃ。
昨日できなかったFeHポテンシャルのメールを出す。
クラスの実装、テスト 。二分木探索のクラスを昔つくったが、とても重宝する。 多面体から頂点へのポインタを作るときとか:

void Polyhedron::InitpVtx()
{
    int f_idx, v_idx;
    Btree bt(1);// 整数1個をキーとする二分木を初期化

    for(f_idx = 0;f_idx < fn;f_idx++)
    {
	Face *f = pFace[f_idx];
	for(v_idx = 0; v_idx < f->vn; v_idx++)
	{
	    Vtx *v = f->pVtx[v_idx];
	    int key = v->idx;
	    bt.Insert(&key, v);
	}
    }

    pVtx = (Vtx **)malloc(bt.NData() * sizeof(Vtx *));
    vn = bt.NData();

    bt.GetAll((void **)pVtx);
}
という感じで簡単にdouble countingを防げる。
2008/09/16
簡単にアップロードできるようになったんでこちらの日記を再開。
本日やること:
ESの結果帰ってくる。残りノード時間と相談し次どうするか決める。 作りかけ鉄水素ポテンシャルを使い、screw dislocation にstrainを かけたとき通常より小さいstrainで転位が動くような水素の位置を みつけ、それを使ってESのVASPで徐々にstrainを増やして緩和させてきた。 水素なしでは16%-17%で動いた。水素無しMDではほぼおなじstrainで動いた。
http://www.jamstec.go.jp/esc/projects/fy2006/46_kaburaki.pdf
http://www.jamstec.go.jp/esc/projects/fy2007/33_kaburaki.pdf
しかし水素いれて15%まで来たがまだ動かないのでかなり冷汗。 まあ水素によって芯の構造がかなり変わってるので、動かなけりゃまあ その構造を成果とする。
鉄水素ポテンシャルの共同研究者に論文のポイントとなる点をまとめて 案として提案する。 メインは、いまあるポテンシャルは古い鉄EAM用なんで新しいのが求められてる、そしてVASPで いろいろ結果が出ているのでそれを採り入れる、という点。VASPの結果が正しい保証は ないが、少なくともVASPやる時の初期配置作成に使えばすげー速く収束する。
多結晶のメッシュを切るためのC++の構造体の整備。
1:ボロノイ分割の頂点にノードを置く。シリアル番号を付ける。 ファイルに番号と座標を吐く。
2:同エッジを等分割しノードを置く。シリアル番号を付ける。各エッジに対応したバッファを作り、 ノード数、シリアル番号*数、位置*数を記憶。 ファイルに番号と座標を吐く。
3:同面をメッシュで切り、生成した新しいノードにシリアル番号を付ける。 まず面の各頂点をバッファに詰め、各エッジのノードのデータも詰める。実際には バッファへのポインタの配列となる。これを二次元Delaunay分割のルーチンに渡して 新しいノードを生成。シリアル番号を降る。 各面に対応したノードデータのバッファを作成。 ファイルに番号と座標を吐く。
4:体積内部について同じことをする。 頂点、エッジ、面を詰めてDelaunay。 四面体の個数と、四面体-頂点のシリアル番号の対応表を作る。 データは各結晶粒ごとにファイルに吐いてから捨てる。
最終データは、ノード数、ノード座標のリスト、四面体数、四面体からノード番号*4へのリスト
緊急事態発生。WIN雑用専用ノートPC、コールサインgollumがBSoDにより再起不能。 ダメコンを急げ。データ損傷、先週バックアップしたので軽微。 WINマシン修復までOpenOfficeで凌げ。Linuxからだとプリンタのパスワード使えないので プリンタが使えん。パワポのライセンスをさがしとけ。修理代が職場の予算から 出るか確認した後修理業者に連絡せよ。見積り金額5万。新しいの買った方が安かったりする。 1 1 1 1