Twin Screw Simulator (Ver.9.0.0)
改良成果資料

2022/09/28
株式会社HASL
Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

○改良成果一覧① ユーザ定義多成分移流方程式の解析機能の実装( Ai + ui∇ ) fi = Bi

Ai , Bi : ユーザ定義任意関数（i=1~n)f i : ユーザ定義未知関数（i=1~n)n : ユーザ定義方程式数
u : 流速ベクトル∇ : ナブラ演算子

用途：

濃度計算
履歴情報計算，滞留時間分布(RTD)計算
化学反応モデル計算(反応押出解析)
繊維破断，液滴分裂，溶融可塑化，etc.

② ユーザ対応改良成果③ TSSカスタマイズ機能

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved1

① ユーザ定義移流方程式の解析機能の実装

Df
=g
Dt

：ユーザが自由に定義可能な移流方程式（連立偏微分方程式）
D / Dt : 物質微分(Material derivative)

Screw extruder

Chemical reactor

Df ∂f=+ ui∇f = g
Dt ∂t

Df ∂f
=
+ ui∇f = g
Dt ∂t

Uniformity

Steady state∂f= 0  f ( z) ?
∂t

∇f = 0  f (t ) ?
Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved2z移流方程式の離散化法∂f=g
Chemical reactor∂t
Backward differencet = n∆tt = 0 t = ∆t

f

n −1

f

z = n∆zz = 0 z = ∆z

f

t

∂f
u
=g
∂z

Screw extruderun

f n − f n −1
= g n  f n = f n −1 + ∆tg n
∆t

n −1

f

n

z

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved⇕ 等価 ∆z∵ ∆t =

Upwind differencef n − f n −1∆z nu
= g n  f n = f n −1 +
g
∆z
u

3

u

定常移流方程式の離散化式
ui∇f = g

fT

qC →T

fL

q L →C

fC ?

定常状態における濃度収支式（2D 風上差分表現）：
fC ( qC →T + qC → R ) − f L qL →C − f B qB →C = vol×g
流量収支式（非圧縮条件）：

qL →C + qB →C = qC →T + qC → R

fR

qC → Rq B →CfB

f L qL →C + f B qB →C + vol×g
fC =q L →C + q B →C
メリット：全体系で構成された連立方程式を解くこと
で，濃度分布を容易に定量化可能．
デメリット：要素分割刻み幅よりスケールの小さい濃度不均一性を表現できない．（解析結果が要素分割密度に依存）

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved4非定常移流方程式の離散化式∂f+ ui∇f = g∂t
後退差分

f n − f n −1
+ ui∇f n = g∆t
1 n −1
 1
 nf
 + ui∇  f = g +
∆t
 ∆t

1
1 n −1
A= ,B = g +f∆t∆t( A + ui∇ ) f = B
Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved5ユーザ定義ルーチン
表1 ユーザが自由にカスタマイズ可能なサブルーチン

ユーザ定義ルーチン名

機能1) initialsetforchem

化学反応種の初期設定及び解析条件の設定2) viscal

粘度計算3) tempcal1

1D温度計算4) chemfscal

化学反応式の左辺及び右辺荷重ベクト
ルの設定

5) chemvariable

化学反応成分の代数的関係式の計算

6) chemwrite

化学反応解析結果のファイル出力

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

6コーリングシーケンス
1) initialsetforchem:
移流方程式（化学反応モデル）の初期条件設定熱流動場との連成解析の場合（非連成解析の場合は，chemcalを非ニュートン
反復計算終了後のループ外でコール）

gamcal:
流速，ひずみ速度の計算
fluxcal:
節点流量の計算4) chemfscal:
移流方程式左辺係数行列と右辺生成項の設定

unfillcal:
部分充満状態の計算
3) tempcal1:
温度の計算

matrixsetchem
有限要素離散化方程式の構成

chemcal:
移流方程式(化学反応モデル）の計算

matrixsolve
有限要素離散化方程式の求解n収束?

5) chemvariable:
解析結果の代数的関係式の計算

y6) chemwrite:
移流方程式解析結果の出力
Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

移流方程式収束反復計算ループ2) viscal:
粘度の計算

化学種数反復計算ループ

非ニュートン熱流動反復計算ループ

heleshaw:
熱流動場の計算

returny収束?

n
図1 コーリングシーケンス

7

公開情報
表2 解析モデル要素情報
変数名

内容

nelem(1)要素数

iinletnodenum

周方向要素分割数

axialdiv(1)軸方向要素分割数

inqc(1)

要素構成節点数(inqc(1)=4)

ndiv

肉厚方向の要素分割数

entab(i,j,1)

要素jを構成するi番目の節点番号(i=1~inqc(1),j=1~nelem(1))

height(i,1)

要素iの肉厚(i=1~nelem(1),単位cm)

vol(i,1)

要素iの体積(i=1~nelem(1),単位cm3)

iinletnodenum×(axialdiv(1)-1)+1
jiinletnodenum×axialdiv(1)=nelem(1)展開i1●(i,j)

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

iinletnodenum

位置(i,j)の要素番号:iinletnodenum×j+i，
i=1~iinletnodenum, j=0~axialdiv(1)-18

表3 解析モデル節点情報
変数名

内容

nnode(1)

節点数

xnode(i,1),ynode(i,1),znode(i,1)

節点iのx,y,z座標(i=1~nnode(1),単位cm)

Iinletnodenum×axialdiv+1

iinletnodenum×(axialdiv+1)=nnode(1)

j

右端の
節点番号
は左端と
等しい

展開
i
1

●
(i,j)

1
iinletnodenum

位置(i,j)の節点番号:iinletnodenum×j+i，
i=1~iinletnodenum, j=0~axialdiv(1)Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved9表4 境界条件,成形条件情報
変数名

内容

tempinlet

流入温度（℃）

pprescribedinlet

流入口設定圧力（MPa）

pprescribed

流出口圧力（MPa)qprescribed

押出量（kg/h)

flowrateinlet

流入口設定流量(cc/s)

rpm

スクリュ回転数(rpm)href(1,ie,1)

要素ieのスクリュ面熱伝達係数(W/cm2/K)

tbound(1,ie,1)

要素ieのスクリュ面境界温度(℃)

href(2,ie,1)

要素ieのバレル面熱伝達係数(W/cm2/K)

tbound(2,ie,1)

要素ieのバレル面境界温度(℃)

iboundt(in,1)

節点inの温度境界条件（0:拘束，非0:自由）ibound(in,1)

節点inの圧力境界条件（0:拘束，非0:自由）Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved10表5 物性情報
変数名

内容

rhoh(i,ie,1)

要素ieの層(i=1~ndiv+1)の密度(g/cm3)

cph(i,ie,1)

要素ieの層(i=1~ndiv+1)の比熱(J/g/K)

vish(i,ie,1)

要素ieの層(i=1~ndiv+1)の粘度(Pa・s)
i=1 : スクリュ表面,i=ndiv+1：バレル表面

表6 解析結果要素情報
変数名

内容

temh(i,ie,1)

要素ieの層(i=1~ndiv+1)の温度(℃)

gam(i,ie,1)

要素ieの層(i=1~ndiv+1)のひずみ速度(J/g/K)

uh(i,ie,1),vh(i,ie,1),wh(i,ie,1)

要素ieの層(i=1~ndiv+1)の流速ベクトル成分(cm/s)

tempe(ie,1)

要素ieの温度

fille(ie,1)

要素ieの充満率

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved11表7 解析結果節点情報
変数名

内容

pres(in,1)

節点inの圧力(Pa)temp(in,1)

節点inの温度(℃)

visn(in,1)

節点inの粘度(Pa・s)gamn(in,1)

節点inのひずみ速度(s-1)

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved12表8 ユーザ定義変数
変数名

内容

chemcnumber

解析対象とする移流方程式の本数（化学種数）

chemvnumber

解析で考慮する配列変数の数

commonvnumber

ルーチン間で共用するスカラー変数の数

chemcname(i)

化学種の名称（i=1~chemcnumber)

chemvname(i)配列変数の名称（i=1~chemvnumber)

chemc(i,ie)

化学種の要素濃度(i=1~chemcnumber,ie=1~nelem(1))

chemcn(i,in)

化学種の節点濃度(i=1~chemcnumber,ie=1~nelem(1))

chempar(i,ie)

化学種依存の要素変数(i=1~chemvnumber,ie=1~nelem(1))

chemparn(i,in)

化学種依存の節点変数(i=1~chemvnumber,ie=1~nelem(1))

commonvpar(i)ルーチン間で共用するスカラー変数(i=1~commonvnumber)

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

13

Test sample ：化学反応種単一成分の濃度計算（理論精度検証）
Explicit
Iteration :20

Df= −κ f , κ = 0.05
Dt

Implicit
Iteration :1

Df= u i∇f l = −κ f l −1
DtDf
= u i∇f l = −κ f l  (κ + u i∇) f l = 0
Dt

Qex :10kg/h

Phead:3MPa

f (0) :1
表9 ケーススタディ条件

Df= −κ f 2 
Dt

Case

Screw rotation speed [rpm]

1
2
3

100
200
400
Explicit
Implicit

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights ReservedDf= u i∇f l = −κ ( f l −1 ) 2
Dt2
Df= u i∇f l = −κ f l  κ f l −1 + u i∇ f l = 0
Dt

( )

14

(

)

理論解:

df
κ
u
= −κ f  f ( z ) = f (0) exp(− z )dzuu, κ > 0
u

f (0)

．．．

f n −1 f n f n +1

i=0

∆z

z

i=n
何れの離散化解もn → ∞, ∆z → 0

前進差分：陽解法uf n +1 − f nn= −κ f n  f n = (1 − κ∆z / u ) f (0)∆z

後退差分（風上差分）：陰解法
nf − f n −1
1


u n
= −κ f n  f n =  f (0)∆z 1 + κ∆z / u 

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

15の極限操作によって理論解と一致x

 1
∵ lim 1 +  = e
x →±∞
 x

Average residence time [s]Flow direction

100 rpm

200 rpm

Screw length [mm]400 rpm

図2 平均滞留時間分布予測結果の比較

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

図3 部分充満状態予測結果の比較

16

Concentration[-]Screw length [mm]図4 軸垂直断面内体積重み付け化学種濃度陰/陽解析結果平均値の比較
Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

17Concentration f [-]Flow direction

100 rpm

200 rpm
Screw length [mm]Theory:f (t ) = f (t − ∆t ) exp(−κ∆t )∆t : Suminf information
400 rpm

図5 軸垂直断面内体積重み付け化学種濃度陰解析結果平均値と理論解の比較
Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

図6 化学種濃度陰解析結果の比較

18

乖離要因分析
Suminf 情報を利用した予測結果

 κ f ( z ) = f (0) exp  −z uave uave =?

TSS output

=
>
<


 κ 
 κ 1f ( z ) = f (0)  exp  −z + exp  −z 
 u fast 2u slow  


1u fast + uslow )
(
2

Case

ufast

uslow

Δu

A

7

5

2

B

9

3

6

C

11

1

10

uave [mm/s]

f(z)/f(0)[-]

表10 ケーススタディ条件

6

TSSの解析結果には，周方向に発生する流量分
布の不均一性が反映されていることが，Suminf情報を利用した計算結果との乖離の要因

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

Screw length [mm]図7 流速差と平均分布量の関係
19

100 rpm

400 rpm
200 rpm
図8 化学種濃度陰解析結果2D展開情報の比較

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

20RTD [sec-1]

滞留時間分布(RTD : Residence Time Distribution)の計算

Df ∂f=+ u i∇f = 0
Dt ∂t

Time [sec]Tracer concentration [-]t = tmint = tmaxt

f out (t ) = fin  RTD(t )dt

fin(t)=fin0t = tmax

fout(t)1 df out (t )RTD(t ) =fin dt

tmax

RTD(t )dt = 1

tmin

t = tmin Time [sec]∵ f out (tmax ) = fintres = 

tmax

tmin

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

21

t×RTD(t )dt

RTD [sec-1]

fout(t) [-]Time [sec]Time [sec]

図9 流出口連続体トレーサー濃度予測結果の比較

図10 滞留時間分布(RTD)予測結果の比較

表11 平均滞留時間と計算所要時間の比較
300

Case Suminf res.: Vol / Qex [sec] TSS output:  ti RTD(ti )∆t [sec] CPU time. [sec]
i =1

1
2
3

65.4
55.6
50.5

63.6
53.7
48.6
(Δt :0.5 sec, 300 cycle unsteady analysis)

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

22

166
226
393

20 kg/h
100 rpm
Full packed

Standard Screw

KD(θ:22.5 °)RTD [sec-1]

Suminf Resav:27.1sec

KD(θ:45 °)KD(θ:90 °)Time [sec]Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

23

Application sample : 二軸スクリュ内における過酸化物によるポリプロピレンの分解反応解析

ポリプロピレン
Polypropylene

プロピレン
Propylene
H

CH3

C

C

H

H

・・・

H

CH3

H

CH3

H

CH3

C

C

C

C

C

C

H

H

H

H

H

H

ぺルオキシドPeroxide
2 R O

R O O R
O-Oの単結合は非常に壊れやすいCopyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

フリーラジカル（遊離基）Free radical

24

・・・

フリーラジカルによる水素原子の引き抜き

・・・

H

CH3

H

CH3

H

CH3

C

C

C

C

C

C

H

H

H

H

H

H

R

・・・

・・・

H

CH3

H

CH3

H

CH3

C

C

C

C

C

C

H

H

H

H

H

・・・

R

+
O H

H

CH3

H

CH3 H

CH3

C

C

C

C

C

C ・・・

H

H

H

H

+
O

高分子鎖のβ切断
H

CH3

H

CH3

H

CH3

C

C

C

C

C

C

H

H

H

H

H

・・・

Transfer
Termination
Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

25

+

H

Free radical reaction kinetic model 1),2)
ペルオキシド(Peroxide)の濃度変化

D[I ]
Dt

= −κ d [ I ]ペルオキシドラジカル（Peroxide radical）の濃度変化

D  R o Dt∞

= 2 fkd [ I ] − κ1  R   (n − 1) [ Pn ]on=2

ポリプロピレン高分子鎖(鎖数 n )の濃度変化

D [ Pn ]Dt= κ1  R o   2  [ Pi ] − (n − 1) [ Pn ]  i = n +1
∞

[ I ] : ペルオキシド濃度[mol/m3]
 R o  : ラジカル濃度[mol/m3][ Pn ] : 鎖数 n の高分子濃度[mol/m3]

κ d : ペルオキシド分解速度定数[s-1]κ1 : 高分子鎖分解速度定数[m3/mol/s]

1) Suwanda, D., Lew, R. and Balke, S. T.: “Reactive extrusion of polypropylene I:
Controlled degradation”, J. Appl. Polym. Sci. , 35(4),1019-1032(1988)
2) Suwanda, D., Lew, R. and Balke, S. T.: “Reactive extrusion of polypropylene II:
Degradation kinetic modeling”, J. Appl. Polym. Sci. , 35(4),1033-1048(1988)

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved26

Closure moment method 3)

MWD (Molecular Weight Distribution）

Mn ≤ Mw ≤ Mz

i番目のモーメント(ith moment)[mol/m3]：

Qi =  ni [ Pn ]n=2

数平均分子量(Number average molecular weight)[kg/mol]：Q1
M n = m0Q0

Frequency∞log (Molecular weight)

重量平均分子量(Weight average molecular weight)[kg/mol]：M w = m0Q2Q1z平均分子量(z average molecular weight)[kg/mol]：Q3
M z = m0
Q2

n : チェインの数

m0 : 繰り返し単位の分子量
ポリプロピレンの場合，プロピレンの分子量：42×10-3kg/mol

3) Tzoganakis, C., Vlachopoulos, J. and Hamielec A. E.,” Production of Controlled-Rheology Polypropylene
Resins by Peroxide Promoted Degradation During Extrusion” , Polym. Eng. Sci., 28(3),170-180(1988)
Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

27モーメント方程式の導出
0-th moment: Q0
D [ Pn ]Dt∞

D [ Pn ]

 Dt
n=2

 ∞

= κ1  R   2  [ Pi ] − (n − 1) [ Pn ]  i = n +1

o

∞
 ∞ ∞

= κ1  R   2  [ Pi ] −  (n − 1) [ Pn ] n=2
 n = 2 i = n +1
o

二重和公式：∞∞∞  [P ] = [P + P
n +1i

n = 2 i = n +1

n+2

+ Pn +3 + iii]

n=2

= [ P3 + P4 + P5 + iii ] + [ P4 + P5 + P6 + iii ] + [ P5 + P6 + P7 + iii ]∞=  (n − 2) [ Pn ] = Q1 − 2Q0
n=2

DQ0= κ1  R o  ( Q1 − 3Q0 )Dt

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved28

1-st moment: Q1
D [ Pn ]
Dt

∞


n=2

 ∞

= κ1  R   2  [ Pi ] − (n − 1) [ Pn ]  i = n +1


Dn [ Pn ]
Dt

o

∞
 ∞ ∞

= κ1  R   2 n  [ Pi ] −  n(n − 1) [ Pn ] n=2
 n = 2 i = n +1
o

二重和公式：∞∞∞

 n  [ P ] =  [ nP + nP
n +1

i

n=2

i = n +1

n+2

+ nPn +3 +iii ]n=2= [ 2 P3 + 2 P4 + 2 P5 + 2 P6 iii ] + [3P4 + 3P5 + 3P6 + 3P7 iii ] + [ 4 P5 + 4 P6 + 4 P7 +iii]
n2 − n − 2
11[ Pn ] = Q2 − Q1 − Q0
2
2
2
n=2
∞

=

DQ1
= −2κ1  R o  Q0
Dt

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved29

2-nd moment: Q2
D [ Pn ]
Dt

∞


n=2

 ∞

= κ1  R   2  [ Pi ] − (n − 1) [ Pn ]  i = n +1


Dn 2 [ Pn ]Dto∞
 ∞ 2 ∞
= κ1  R   2 n  [ Pi ] −  n 2 (n − 1) [ Pn ] n=2
 n = 2 i = n +1
o

二重和公式：∞∞∞

 n  [ P ] =  n P + n P
2

2

n=2

i = n +1

2

n +1

i

n+2

+ n 2 Pn +3 +iii n=2= [ 4 P3 + 4 P4 + 4 P5 + 4 P6 iii] + [9 P4 + 9 P5 + 9 P6 + 9 P7 iii ] + [16 P5 + 16 P6 + 16 P7 +iii ]

DQ2
1
 1

= κ1  R o   − Q3 + Q1 − 2Q0 
Dt
3
 3


Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved∞ (n − 2)(n − 1)(2n − 3)= 
+ (n − 2)(n − 1) + (n − 2)  [ Pn ]
6

n=2 
3
2
∞
 2n − 3n + n − 6 =  [ Pn ]
6
n=2 


30

Closure relation4)

2Q2
Q3 =
2Q2Q0 − Q12
Q1Q0

(

)

4) Hulburt, H. M. and Katz, S. ”Some Problems in Particle Technology : A Statistical
Mechanical Formulation”, Chem. Eng. Sci. , 19,555-574(1964)
Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved31ペルオキシドラジカル（Peroxide radical）の濃度変化

D  R o Dt∞

= 2 fkd [ I ] − κ1  R   (n − 1) [ Pn ]on=2

= 2 fkd [ I ] − κ1  R o  ( Q1 − Q0 )Steady state approximation

2 fkd [ I ] − κ1  R o  ( Q1 − Q0 ) = 0∵ Q1 − Q0 >> 1

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved32

PPの分解プロセスを表現する化学反応モデル
ペルオキシド(Peroxide)の濃度変化

D[I ]
Dt

= −κ d [ I ]モーメント方程式

 Q − 3Q0 
DQ0= 2 f κd [ I ] 1

Dt
Q
−
Q
0 
 1
 Q0 
DQ1= −2 f κ d [ I ] 

Dt
Q
−
Q
0 
 1
1
 1

−
Q
+
Q
−
2
Q
3
1
0


DQ2
3
= 2 f κd [ I ] 3

Dt
Q
−
Q
1
0




2Q2
Q3 =
2Q2Q0 − Q12 )(Q1Q0

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

33

efficiency
f : Initiator(0.6~1)

粘度と各種平均分子量との関係

42g/mol

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

34

log η0η0 ∝ M w3.4
Slope=3.4

Slope=1
log Mw

5) Fox, T. G. and Flory, P. J. ,” Viscosity—Molecular Weight and Viscosity—Temperature
Relationships for Polystyrene and Polyisobutylene“, J. Am. Chem. Soc. ,70,2384-2395(1948)

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

35

Modified Cross model6)
η=

η0
,
1− n
1 + C (η0γɺ )T 
η0 = B exp   , T0 
B = k1M w3.7，
Mz
C = k2Mw

非噛合型異方向回転二軸スクリュの
2D FAN法を利用した等温完全充満
部分領域の解析

6) Fukuoka, T. and Min, K. , “Flow Analysis a Reactive Processing in a Twin Screw
Extruder”，Seikei-Kakou, 7(8),521ー528(1995)

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved36解析モデルBarrel:200℃，HTC(1000W/K/m2)
Phead:3MPa

Qex : 10kg/h
200℃

100 rpm
u：流動情報T：温度情報C：化学成分情報
η：粘度情報表12 ケーススタディ条件

Case

Thermal condition
for chemical model1

Isothermal (200℃)

Cross

u

→ C, T=200℃

2

Non-isothermal

Cross

u,T

→ C

3

Non-isothermal

Modified Cross6)

u,T

→ C

Viscous modelη (γɺ, T , C )Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved37図11 フローカーブのCross モデルフィット（Material fit data base, PP_D3)
Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved38各種モーメント流入境界値
ρ750kg / m3

解析条件
入力基本情報ρ : 750 kg/m3

M n : 74400 kg/mol
M w : 356000 kg/mol
M z : 1176000 kg/mol

Q0 =

Mn

=

74400kg / mol= 0.010mol / m3

750kg /m3
Q1 === 1.7857×104 mol / m3
−3
m0 42×10 kg / molρ

MW Q1 356000kg / mol×1.7857×104 mol / m3
Q2 ==42×10−3 kg / mol
m0= 1.5136×1011 mol / m3

M z Q2 1176000kg / mol×1.5136×1011 mol / m3
Q3 =
=
m0
42×10−3 kg / mol= 4.2381×1018 mol / m3

計算条件

ペルオキシド流入境界値[ I0 ] :非ニュートン反復計算数：10
化学反応成分：4 種

1000 ppm=10-3 mol/m34
1.49×1012分解速度定数[s-1] kd = 1.98×10 exp  − T + 273.15 
Initiator efficiency[-]f = 1.0

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

39

化学反応モデル反復計算数：10
化学反応モデル緩和係数：1.0

計算時間：918 sec
IntelR CoreTM i7-9700
CPU @3.00 GHz

Residence time [s]Temperature [℃]Screw length [mm]
図13 温度分布予測結果の比較

Screw length [mm]
図12 滞在時間分布予測結果の比較

Case1,2Case 3
図14 充満率分布予測結果の比較
Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved40

Case1,2Case 3
完全充満領域拡大

未充満領域拡大

図15 圧力分布予測結果の比較
部分充満領域の充満率は，粘度には依存せず，押出量と回転を一定とした場合，変化しな
い．しかし，未充満領域の押出方方向への分布状態は，粘度に依存して変化する．Case3では，PPの分解に伴う粘度低下が考慮されているため，下流側では先端背圧3MPaの抵抗の影響がCase1,2と比較して高まるため，完全充満状態の領域が拡大する．逆に上流側では，隣接する完全充満領域の圧力が低下するため，その上流側の未充満領域が拡大する．

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved41化学反応解析結果の妥当性についての検討

Residence time [s]TSS:Steady analysis systems of partial
differential equations
Unsteady analysis of systems ofordinary differential equations
d [I ]dt= −κ d [ I ] Q − 3Q0 
dQ0= 2 f κd [ I ] 1
dt
 Q1 − Q0 

 Q0 
dQ1= −2 f κ d [ I ] 

dt
Q
−
Q
0 
 1
1
 1

− Q3 + Q1 − 2Q0 dQ2
3= 2 f κd [ I ] 3dt
Q1 − Q0

Screw length [mm]TSSモデルの軸方向分割数：537分割(周方向76分割，要素数：40812)
Unsteady analysis condition :Δt=0.01s, 12000 cycles

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

Q3 =

42

2Q2
2Q2Q0 − Q12
Q1Q0

(

)

Moment Q0 [mol/m3]

Peroxide concentration [mol/m3]

Screw length [mm]Moment Q1 [mol/m3]

Moment Q2 [mol/m3]

Screw length [mm]Screw length [mm]Screw length [mm]

図16 TSS解析結果と非線形連立常微分方程式の非定常解析結果の比較
Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved43妥当性検討結果
TSSを利用した定常非線形連立偏微分
方程式の反復計算結果≅ 非定常非線形連立常微分方程式の非定常解析結果

化学種の濃度分布に観られる乖離は，滞在時間の表現の差に起因：TSS；未充満部Δt短，充満部Δt長，非定常解析結果；Δt=0.01 s一定
当化学反応モデルでは，1次モーメントQ1は殆ど変化しない．
Q1 =17857 [mol/m3]∞

m0Q1 = m0  n [ Pn ] = ρ
n=2m0 : 42 ×10−3 [kg/mol]ρ = 750 [kg/m3]
PPの分解に伴って粘度は変化するが，密度は殆ど変化しない：
非圧縮粘性流体としてのモデル化が妥当
Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

44

Moment Q0 [mol/m3]

Peroxide concentration [mol/m3]

Screw length [mm]Moment Q1 [mol/m3]

Moment Q2 [mol/m3]

Screw length [mm]Screw length [mm]Screw length [mm]図17 ペルオキサイド濃度と各種モーメントの予測結果の比較
Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

45

Average molecular weight [kg/mol]

MWD parameter [-]Screw length [mm]Screw length [mm]

図18 各種平均分子量とMWDパラメータの予測結果(Case 3)
Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

46

Case2 : Cross model

Case3 : Modified Cross model
図19 粘度分布予測結果の比較
Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved47

B:100 mm

D:632 mm

C:200 mmlog(η) [Pa・s]

A:0 mmlog(γɺ ) [s-1]
図20 スクリュ内各所でのフローカーブ予測結果(Case 3,200℃)
Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved48ユーザ定義ルーチン
Initialsetforchemのユーザ定義部
当ルーチン内で利用する変数の型宣言

解析対象とする移流方程式の数＝4

化学反応計算の初期値（境界値）設定

代数的関係式の数＝6コモン変数の数＝2

ポスト項目の見出し

各種平均分子量の初期値設定
Viscalと共用する変数定義

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved49化学種要素変数の初期設定

流入口(in=1~iinletnodenum)に
濃度境界値を設定

化学種節点変数の初期設定

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved50コンターランク数値の指数表示．
テキストボックスには小数点以下のの桁数を指定．

Initialsetforchemでユーザが定義したポスト項目見出し．

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

51

chemfscalの内容
当ルーチン内で利用する変数の型宣言
Initiator efficiency
ペルオキシド分解速度定数
パラメータ

TSSが解析対象とする移流方程式：

D

chemf
(
ic
,
ie
)
+

 f (ic, ie) = chems (ic, ie)
Dt 

ペルオキシド分解速度定数の計算
tempe(ie,1):要素ieの温度

D[I ]

D

= −κ d [ I ]   κ d +[I ] = 0
Dt
Dt



 Q − 3Q0 
DQ0= 2 f κd [ I ] 1

Dt
Q
Q
−
1
0


 Q0 
DQ1= −2 f κ d [ I ] Dt Q1 − Q0 1
 1
 − 3 Q3 + 3 Q1 − 2Q0 
DQ2= 2 f κd [ I ]Dt
Q1 − Q0

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved52

変数対応関係[ I ] : chemc(1, ie)
Q0 : chemc(2, ie)
Q1 : chemc(3, ie)
Q2 : chemc(4, ie)chemvariableの内容

当ルーチン内で利用する変数の型宣言
プロピレンの分子量：4.2D-02 kg/m3

Closure relation

変数対応関係Chempar(1,ie) : Q3
Chempar(2,ie) : 数平均分子量Mn
Chempar(3,ie) : 重量平均分子量Mw
Chempar(4,ie) : Z平均分子量Mz
Chempar(5,ie) : Mw/Mn
Chempar(6,ie) : Mz/Mw

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

53

nonnew : 非ニュートン識別番号
1:Newtonian
2:Power law
3:Log polynomial
4:Carreau
5:Carreau Yasuda
6:Cross

viscalの内容

Tb T + 273.15 

η0 = a exp a→bcoef (Material fit パラメータ),
Tb→tref (Material fit パラメータ),
T→temhi:要素ieの層iの温度3.7

 M (t ) η0 (t ) =  w  η0 (0) M w (0) dc =M z (t ) / M w (t )M z (0) / M w (0)
η=

η0
1− dn

 η * gamm 
1 + dc *exp  0 rcoef 

dc→反応によるshear thinning 性の変化を表現するパラメータ
rcoef : Material fit パラメータ
gamm : 要素の層iのひずみ速度
vish(i,ie,ib) : 要素ieの層iのひずみ速度(ib=1)
omgc : 緩和パラメータ(0.5)

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved54ユーザ定義ルーチンの利用方法
1) 公開されているソースコードを用途向きに書き直す．
2) Fortran complier を利用して，実行ファイルmultiprofilesimulator.exe 作成．
3) 実行ファイルmultiprofilesimulator.exe を¥bin¥x86¥Systemフォルダーにコピー．
4) GUIに新しく設けられたCall user’s routine チェックボックスをチェック状態とする．
5) 従来通りの方法に従ってプログラムを実行．

チェック状態とするとユーザ定義ルーチンをコール
チェック状態 ： 熱流動解析との連成（各非ニュートン反復計算毎にコール）
非チェック状態 ： 熱流動解析と非連成（熱流動解析後にコール）
ユーザが定義した連立非線形偏微分方程式の反復収束計算回数（線形の場合は１）反復収束計算の緩和係数
離散化方程式の反復法マトリックスソルバー（Gauss-Seidel法）の最大反復計算回数
反復法マトリックスソルバーの反復計算終了
基準値εAx − B≤ ε → 反復計算終了
B

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved55

② ユーザ対応改良成果(1) 解析のリスタート機能(2) 解析の収束状況確認ファイル出力 (.calinf)(3) 液滴分裂/合体解析の調整パラメータ追加(4) Option Form の設定パラメータ保存

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved56

(1) 解析のリスタート機能
本機能を利用すると, 熱流動反復計算ループ終了時の解析結果を再利用することで, 計算時間の短縮を図る
ことができます.

解析条件の設定

非ニュートン熱流動反復計算ループ

熱流動場の計算

リスタート解析手順
解析開始
解析条件の設定
新機能

リスタート用の熱流動解析結果読込み
反復回数の増加 ?

y
非ニュートン反復計算ループ
回数増加分の実施
解析オプション機能／熱流動場と独立

解析開始

通常の解析手順n
熱流動反復計算を省略して,解析オプション機能（繊維破断解析,
脱揮解析,
液滴分散解析, etc.）のパラメータ調整が可能.

流速，ひずみ速度の計算
節点流量の計算
部分充満状態の計算
温度の計算
解析オプション機能／熱流動場と連成

n
新機能

収束?

y

リスタート用の熱流動解析結果出力
解析オプション機能／熱流動場と独立

解析結果出力 (.twinres など)

解析結果出力 (.twinres など)

解析終了

解析終了
Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved57

(1) 解析のリスタート機能

TSS Ver.9.0.0 Analysisタブ 初期画面

利用手順
(事前準備)
0. TSS Ver.9.0.0 で通常（従来通り）の解析を実施します.
⇒ リスタート用の熱流動解析結果を得ることが目的
解析条件ファイル名.recalc が出力される.

(リスタート解析の条件設定方法)
1. Analysis タブ下段に設置された Recalculation欄の□ Recalculate にチェックを入れる.
⇒ Select ボタンをクリックし, 0. で解析済みの解析条件ファイル.tscal を選択する.

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved58

(1) 解析のリスタート機能

TSS Ver.9.0.0 Analysisタブ / リスタート用解析条件ファイル読込後

利用手順
2. 手順1. で解析条件ファイル(.tscal)を読込むと, 該当ファイルの解析条件内容が設定されます.

5.

3. リスタート解析で得られる解析結果ファイル名(.twinres, etc.)
を別名で保存する場合は, Result File Name欄に設定された元の名前を変更します.
（変更しない場合は上書き保存されます.)4. リスタート解析で得られる解析条件ファイル名(.tscal)を
4. の解析結果ファイル名と同じにする場合（推奨）には,Calculation Control Information File Name欄 右のCopyボタンをクリックすると, 手順3. で設定した名前が反映されます.

4.

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved59

(1) 解析のリスタート機能

TSS Ver.9.0.0 Analysisタブ / 名前変更および条件変更後

利用手順
5. リスタート解析の目的となる条件変更を行います.

変更を推奨【変更可能な項目】
・非ニュートン反復計算回数（元よりも多い場合に追加計算）
・解析オプション機能の各種パラメータ- Melting/Morphology model タ ブ- Binary System Condition Set タ ブ- Forming model タ ブ- Surface renewal model タブ

変更可能【変更不可能な項目】
・材料データ(.pro), メッシュモデルデータ (.twinmsh)
・非ニュートン熱流動反復計算ループに直接関係するパラメータ- Analysisタブで読込専用（灰色）になっている項目- Thermal Boundary Condition （熱境界条件）
・Dynamic解析はリスタート機能の対象範囲外になります.

6. 手順5. で条件変更後, 通常の解析と同様に, Analysis タブ
下段のSaveボタンをクリックして解析条件ファイルを保存し,Executeボタンをクリックすると, 元の解析条件ファイルの熱流動結果(.recalc) を利用したリスタート解析が実施されます.

変更を推奨

7. 解析条件保存および解析実行（従来通り）
Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved60

(1) 解析のリスタート機能／解析事例(A) 基準解析ケース／解析条件
・樹脂データ ：HDPE_B3 (Material fit データベース)
・スクリュモデル: バレル径20mm, L/D=40 (1600mm)
Flow direction

・バレル温度: 200℃ Max. (スクリュ入口30℃)
・解析条件： 未充満解析(Un-fill analysis)
押出量 40 kg/h,
出口圧力 1.0 MPa
スクリュ回転数 100rpm,
・オプション機能： 繊維破断解析(Fiber Attrition Model)

破断頻度 Pi を決定する主要パラメータ(詳細は手引書参照)Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved61

(1) 解析のリスタート機能／解析事例(A) 基準解析ケース／解析結果○充満率分布Flow direction○熱流動解析結果の概要(fam0002.twininf)○応力分布(kPa)○重量平均繊維長(um)Flow direction

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved62

(1) 解析のリスタート機能／解析事例○熱流動解析結果の概要(fam0004.twininf)(B) リスタート解析ケース(A)の熱流動
解析結果を使用.

熱流動反復計算を省略することで, 計算時間が1,102 ⇒ 163 sec に減少～～～○重量平均繊維長(um)(B)の名前に変更(A)の結果を選択

Cb=0.004に変更

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved63

(2) 解析の収束状況確認ファイル出力 (.calinf)
本機能では, 非ニュートン熱流動反復計算回数毎の, (1) スクリュ出口の押出量 (Mass flux, kg/h)
および, (2) スクリュ出口の平均温度 (℃) を, 解析結果ファイル名.calinf ファイルに自動出力するようにしました.
.calinf ファイル（エクセルのカンマ区切りで開いた場合）(1) 押出量(2) スクリュ出口平均温度↑ p.61 の (A)基準回帰ケースの例 (fam0002.calinf)

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved64

(2) 解析の収束状況確認ファイル出力 (.calinf).calinf で出力される押出量 Mass Flux は, スクリュ出口の解析結果になりますが,
Ver.9.0.0 では熱流動解析の計算ルーチンを改良し, 解析精度（流量収支）を向上させました.
○解析事例／サイドフィード解析
副材: 5.0 kg/h
500Pa・s一定

・スクリュモデル: バレル径20mm, L/D=15 (600mm)
・スクリュ回転数: 100rpm ⇒ 平均充満率 9.9%400rpm ⇒ 平均充満率 2.6%
出口圧力
0.1 MPa

主材: 5.0 kg/h
1000Pa・s 一定
Flow direction○解析結果／押出量（解析結果）のスクリュ軸方向依存性 (.binaryinf の Flux 出力結果から算出)

理論解 10.0 kg/h
旧バージョンでは,充満率が極端に小さい場合に流量収支に誤差が見られるが,
新バージョンでは,
平均充満率が2.6% (400rpm)と極端に小さい場合でも理論解と一致した.

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved65

(3) 液滴分裂/合体解析の調整パラメータ追加
本機能では, Ver.8.0.0で改良された液滴分裂解析において, キャピラリー数の算出に利用されるひずみ速度 γɺ の抽出範囲を調整するための設定項目を追加しました.
(理論の詳細は, TwinScrewSimulatorVer.8.0.0改良成果資料.pptx の p.44～を参照ください),
キャピラリー数

Ca =

τ
σ /R

=

ηm γɺ
σ /R

Ca が大きいほど変形 or 分裂し易い．

Ca : キャピラリー数
τ : 応力σ : 表面張力係数
R:

液滴半径

ηm : マトリクス粘度γɺ : ひずみ速度

液滴初期半径計算開始Z位置
液材体積分率表面張力係数液滴平衡径
無次元時定数のモデルパラメータ⇒ 設定値以上の肉厚要素を
ひずみ速度の抽出範囲とする.
（肉厚下限値）
液滴材料データ

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

66

tb* = aλ b ,a = 91.41,
b = 0.3397

(3) 液滴分裂/合体解析の調整パラメータ追加○解析事例

・スクリュモデル: バレル径20mm, L/D=8.05 (322mm)
最小肉厚: 0.5 mmt
出口圧力
2.0 MPa

主材: HDPE
押出量 6.0 kg/h
回転数 250rpm
Flow direction○解析条件: Thickness criteria （肉厚下限値） 違いの 3条件(A) Thickness criteria: 0 mm
(デフォルト設定, スクリュ全領域の γɺ が抽出対象)(B) Thickness criteria: 1 mm
⇒ 肉厚 1mm 以上の要素が γɺ の抽出対象.
(C) Thickness criteria: 2.5 mm⇒ 肉厚 2.5mm 以上の要素が γɺ の抽出対象.
(B)の設定例

液滴材:
50 mPa・s 一定⇒ 粘度比 10-5初期半径 2000um
Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

67

(3) 液滴分裂/合体解析の調整パラメータ追加○解析結果:【 Thickness criteria vs. ひずみ速度 γɺ 】

【 Thickness criteria vs. 液滴分散径 (Droplet radius) 】

Thickness criteria が大きいほど, ひずみ速度 γɺ （周方向平均値）
が小さくなるため, Ca が減少し分散径が小さくなる.
また非ニュートン粘性流体の場合, γɺ が小さくなるとマトリクス(主材) の粘度 η m および 液滴材の粘度 η dも変化するため, 粘度比 λ が変わり, 分裂状態も変化する.

【充満率分布】

λ=Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

68

ηd
ηm

粘度比が大きいほど分裂し難い．

(4) Option Form の設定パラメータ保存
本機能では, Option Form タブで変更された設定値をファイル保存し, 次回起動時にも設定値を反映させることができます.

フォームを閉じる
変更設定値を保存するボタン
デフォルト設定値に戻すボタン

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved69

(4) Option Form の設定パラメータ保存
利用手順
1. 必要に応じて, Option Form の設定値を変更します.
(従来通り)
2. 変更した設定値を保存したい場合は,フォーム右下のSave ボタンをクリックします. 以降のTSS起動時には保存時の設定値が初期設定されます.
3. 変更前のデフォルト設定値に戻したい
場合は, フォーム右下の Default Set
ボタンをクリックします.

2.
3.

(変更例) メッシュサイズが大きく, 解析
結果描画時にメモリエラーが発生する
場合には, デフォルト値 300,000よりも大きくすると改善する.

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

70(変更例) メッシュサイズが大きく, メッシュデータ保存時に要素数または節点数がデフォルト値の50,000以上の場合,
最大要素数または最大設定数を大きくすると（ブロック数はエレメント数に応じて削減する）
メッシュデータ保存時のメモリエラーが解消する.

③ TSSカスタマイズ機能
ユーザ定義ルーチン ／ Fortran compiler の利用方法
TSS Ver.9.0.0 で利用可能な開発環境／下表2種類のいずれかを利用してください.
開発環境（compiler）

解析で使用されるマトリクスソルバ

使用PCのOS (Operating System)① Intel Fortran

SMSソルバライブラリ (TSS標準ソルバ)

Windows② GFortran (GNU Fortran)

当社自社開発ソルバ

Windows (MinGWを使用)

インストールおよびコンパイル方法（当社実施事例）は以下の別資料を参照ください.
- TwinScrewSimulatorVer9.0.0(2022)TSSカスタマイズ環境設定方法.pdf（参考）TSS Ver.9.0.0 /GPU対応版(Ver.8.0.0 から提供) はユーザ定義ルーチンのサポート対象外になりますが,
ユーザ定義ルーチン以外の解析には使用可能です.
開発環境（compiler）

解析で使用するマトリクスソルバ

使用PCのOS(operating system)

NVIDIA nvfortran

CUDAライブラリ（TSS/GPU版標準ソルバ）

Linux (PGI Fortranは2021年に販売停止)* TSS/GPU対応版が利用可能なPC環境は, 下記2点を満たす必要があります.
1. NVIDIA グラフィックボードが搭載されていること.
2. GPUのCUDA Computer Capability (グラフィックスボードのハードウェアバージョン) が7.5以下であること.
当社社開発環境 GPU：Nvidia GeForce GTX 1660 Ti, CUDA Computer Capability: 7.5, Compiler: PGI Community Edition

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

71

TSS Ver.9.0.0 フォルダ構成
TwinScrewSimulatorVer9.0.0¥bin¥x86 内にに3つのSystemフォルダが存在します.

フォルダの内容と対象ユーザ
標準仕様

Systemver.9.0.0【TSS/標準版】 : 通常解析を実施されるユーザ、①Intel Fortranを利用されるユーザ
Systemver.9.0.0_GPU 【TSS/GPU版】 : GPUを利用されるユーザ（利用可能なPC環境はp.52参照）
Systemver.9.0.0_GFORT: ②GFortranを利用されるユーザ

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved72