第154回 講演会「混合・混練評価にCAE技術はどこまで使えるのか？」
2016/08/04 プラスチック成形加工学会

二軸スクリュ押出機の2.5D FEM熱流動シミュレーションと実験検証
谷藤 眞一郎株式会社HASL〒179-0041 東京都練馬区石神井町3-30-23
石神井ウエスト201
TEL 03-5923-6988 FAX 03-5923-6987
E-mail:tanifuji@hasl.co.jp
URL:http://www.hasl.co.jp
Copyright© 2010- Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

AGENDA◇技術開発の背景

Page
2◇二軸スクリュ押出機内成形現象の定量化法
・二軸スクリュモデリング技術(Template & UV mapping technology)

5・Generalized Hele-Shaw flow 2.5D FEM formulation

10・未充満状態の定量化法 (Pressure downstream update scheme)

15・溶融可塑化モデル(Concentrated suspension model)

22◇2.5D FEM解析技術の適用事例
・繊維破断解析(Fiber attrition model)

33

・二軸スクリュ押出機内充満状態の実験検証解析

39

◇まとめ
Copyright© 2010- Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

431

◇技術開発の背景
Commercially available simulation software for the twin screw extruder

Copyright© 2010- Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved2数値解析技術の特徴.

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

3

HASL Twin Screw Simulator 技術開発目標
・Robust to use, Easy to use, User-friendly (2D FANメリット踏襲)
・ Short computational time :数分～数時間以内(2D FANメリット踏襲)
・ Integrated simulation: 全域解析(2D FANメリット踏襲)
・ Detail & User customized analysis : (close in on 3D FEM,FVM)

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

◇二軸スクリュ押出機内成形現象の定量化法
二軸スクリュモデリング技術
Self Wiping Twin Screw & Kneading disk*)

 CL 
 2 Rs ψ s = cos −1 

α=r(θ)

πnt− 2ψ sr (θ ) = C L2 − Rs2 sin 2 θ − Rs cos θθz=L pθ
2π

C L : 軸間距離Rs : スクリュ半径
α : チップ角*)参考文献：
M. L. Booy, Polym. Eng. Sci., 18,973(1978)
Copyright© 2010- Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved2ψ s : 噛合角
nt : 条数
L p : スクリュピッチ
z : スクリュ長5Template

スクリュ形状パラメータの数値
入力による容易なモデリング

Copyright© 2010- Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved62.5D Hele-Shaw FEM mesh

Co rotating intermeshingself wipe type twin screw

2.5D Hele-Shaw FEM meshes

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

7Thickness information
mapped on 2.5D HeleShaw FEM meshes

UV mapping technology

部分的にピンミキシングエレメントを配置

編集
2.5D 解析モデルのUV展開モデル（Unwound model)のインポート
Copyright© 2010- Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

解析

要素単位で肉厚（流路深さ）を任意編集（当情報を2.5D解析モデルに写像）
8

Flow direction

Flow direction
Velocity distribution in a twin screw extruder with pin mixing element

Stream lines around pin mixing element(Unwound display)
Copyright© 2010- Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

9Generalized Hele-Shaw flow 2.5D FEM formulation
流体支配方程式（円柱座標系）
vr<<vθ, vzを仮定
連続方程式：

1 ∂vθ ∂v z+=0r ∂θ ∂z
運動方程式θ成分：1 ∂  2  ∂vθ vθ   1 ∂p
 r η−  =r 2 ∂r   ∂r r   r ∂θ
運動方程式 z 成分：
Fluid element between cylinders

1 ∂  ∂vz  ∂p
 rη
=r ∂r  ∂r  ∂z
Copyright© 2010- Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved10圧力変数/座標要素内補間補間関数

4p = ∑ φα pα ,11(1 − ξ L )(1 − η L ), φ 2 = (1 + ξ L )(1 − η L ),44
11
φ3 = (1 + ξ L )(1 + η L ), φ 4 = (1 − ξ L )(1 + η L )44

φ1 =α =144

α =1α =1

θ = ∑ φα θ α , z = ∑ φα zα , r = (Rb − R s )ζ L + R s

2.5D FEM Hele-Shaw flow element for the twin screw extruder
Copyright© 2010- Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

11

(eθ)

Qα = − S αβ + S αβ p β + Dαe

節点を通過する流量ez

圧力勾配流れの流量寄与

Qαe = ∫ φα ( nθ qθ + n z q z ) dΓ,∂φα ∂φ β∂φ ∂φ βez
dS , Sαβ= Sz ∫ αdS ,∂θ ∂θ∂z ∂zSeSeeθ
Sαβ= Sθ ∫牽引流れの流量寄与
β c2 α c2 11z, S =  δ c −,S =  α c −γc βc 4
4
Dθ =∂φαdS∂θSeDαe = Dθ ∫3Rb 1
Rb r
1dr , γ c = ∫dr , δ c = ∫dr
α c = ∫ dr , β c = ∫Rs ηRs η r
Rs η r 3
Rs ηr

RbCopyright© 2010- Hyper Advanced Simulation Laboratory Co., Ltd. All Rights ReservedeθΓeRbe12Ω  2 βc  Rb − 2 γc 当定式化の特長
・連立方程式により求めるべき変数が圧力と温度変数になるため計算負荷が小さい。
・圧力と温度が求められれば、肉厚方向の線積分で流速、ひずみ速度、粘度分布などの諸量が容易に算出可能。
・流入出収支（連続方程式）を高精度に表現。

Q4Q3
要素

どのような圧力計算値を用いても、要素内節点流量の総和は０44Qα = 0∑αQ ∑ φα = 1α =1=1Q2

Q1q4q3●q1q2

連立方程式より求められる圧力計算値を用いれば、内部節点流量の総和は０●内部節点Copyright© 2010- Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

134qα = 0∑α=1Theoretical verification for a concentric annular flow

11   11 vθ = Ωr  2 − 2  /  2 − 2 Rs   Rb Rs r
1 ∂p 2r − Rs24η ∂z(log r − log Rs )1 ∂p 2( Rb − Rs2 )−4η ∂z(log Rb − log Rs )vz =Copyright© 2010- Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

14()

未充満状態の定量化法 (Pressure downstream update scheme)
充満解析
大気圧0MPa（指定）

スクリュ背圧10MPa（指定）押出量20kg/h(計算値)
スクリュ背圧5MPa（指定）

大気圧0MPa（指定）押出量40kg/h(計算値)

未充満解析
大気圧0MPa（指定）？
押出量20kg/h(指定)Copyright© 2010- Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

15

スクリュ背圧5MPa（指定）スクリュ内流量バランス式
WVb H
WH 3  dp 
Q=sin( 2θ ) − 412η  dx 
牽引流量QdQ:チャネル内流量 p:圧力
W:チャネル幅
η:流体粘度H：チャネル深さ
x:スクリュ長さ
Vb：バレル相対速度
θ：螺旋角

圧力勾配流量Qp

・牽引流量Qdは、スクリュ形状パラメータと回転数のみによって決定される。
・ dp/dx>0 → Q<Qd , dp/dx<0 → Q>Qd・ スクリュ形状、回転数、流量、材料物性に変化が無ければ、当バランス方程式は、p→p+Δp(一定値）の変換に対して不変。

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

・ スクリュ形状、回転数、流量、材料物性に変化が無ければ、当バランス方程式は、p→p+Δp(一定値）の変換に対して不変。
10MPa:20kg/hの特性は、Δp=-5MPaとすることで5MPa：20kg/hの特性と一致

Pressure (MPa)

10MPa:20kg/h

10MPa

5MPa:40kg/h
5MPa
5MPa:20kg/h

負圧領域

Screw length (L/D)

スクリュ内圧力分布

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

17

未充満状態判定法
圧力勾配

流量バランス圧力

dp/dx<0Q>Qd

必然的にp>0充満1)

dp/dx>0Q<Qdp>0充満2)dp/dx>0Q<Qd

p<0をp=0に補正

未充満3)

状態

1) Q>Qdで未充満状態は有り得ない。
2) Q<Qd , p>0 : 充満状態にあり、圧力勾配流量（バックフロー）が、流量を牽引流量から減じている状態。
3) Q<Qd , p<0 : 未充満状態と判定し、p=0,W=fWに補正、f=Q/Qd:充満率と定義（FAN法の場合）。
Copyright© 2010- Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved18

Z

1st step pressure distribution
Screw head position

Screw inlet position

Z：圧力更新の着目要素
Corrected pressure distribution

：風下要素

f=1.0
f=0Zp ( z − ∆z , θ − ∆ θ ) = p ( z , θ ) −Filling ratio distribution∂p∂p∆z −∆θ∂z∂θ

Pressure downstream update scheme

2.5D FEM

FAN

未充満領域を想定した圧力計算法
Copyright© 2010- Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved19Pressure (MPa)

10MPa

5MPa

Down stream distance (mm)Comparison of predicted pressure distribution between filled and unfilled analysis
Copyright© 2010- Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

20

Pressure (MPa)
0.0

1.0

2.0

3.0

4.0

5.0

Flow
direction

Filling ratio
0.0

0.5

Flow
direction

Predicted pressure and filling ratio distribution
Copyright© 2010- Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

211.0溶融可塑化モデル(Concentrated suspension model)
未溶融ペレット（球体にモデル化）
溶融体

未溶融ペレットを含む溶融体を粘度ηの濃厚サスペンションとしてモデル化η = η rηl
サスペンション粘度

相対粘度

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

溶融体粘度

Einstein model (希薄系)η r = 1 + [η ]φKrieger-Dougherty model*)(濃厚系)φ 
η r = 1 −  φm [η ] : 固有粘度（球の場合2.5)
φ : 粒子体積濃度φm : 最密充填体積濃度φm =π18− [η ]φm≅ 0.74

Face centered cubic

ケプラー予想：球に対してこの配置が最密充填（1611年）、2014年8月10日にヘイルズ（ミシガン大学）らによって証明完結

φm ≅ 0.64

Random close packing*)参考文献:I. M. Krieger and T. J. Dougherty, Trans.Soc.Theol.,3,137-152(1959)
Copyright© 2010- Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved23単位体積V内に半径Rの球がn個含まれていることを想定すると、

ρb43φ (0) = n πR(0) / V =3ρs

ρ s : 固体密度ρ b : 嵩密度１）二軸スクリュ内に投入されるペレットは、下流側に搬送される過程で溶融が促進されて大きさが変化するが、単位体積V内の個数は変化しないと仮定する。
33

 R (t )  R (t )  ρ b φ (0) = φ (t ) = 
 R ( 0)  R ( 0)  ρ s2)ペレットは、周囲の溶融樹脂（温度Tl）からの熱流束の影響を受けて溶融が促進される。

Copyright© 2010- Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved24球対称非定常熱伝導方程式

ρ s : 密度Ts : 温度

Tav:平均温度t : 時間
κ s : 熱伝導率 r : 半径

(ペレットの温度状態を決定する支配方程式) C ps : 比熱DTs κ s ∂  2 ∂Ts ρ s C ps= 2rDtr ∂r  ∂r hsl：ペレット/溶融樹脂間熱伝達係数
Tl：溶融樹脂温度 Tm:溶融温度 λ：潜熱

昇温フェーズ（R=R0：初期半径)

初期時刻t=00<t<tm

溶融フェーズ(DR/Dt<0)

溶融開始時刻tmtm<tT ( R, t m ) = Tm

温度境界条件：∂T∂Ts
κs s= hsl (T f − Ts ) ,=0r=R∂r r = R∂r r =0

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

温度境界条件： Ts r = R = Tm , ∂Ts = 0∂rMelting front 運動学的条件式：
κs25∂T
DR= hsl (T f − Tm ) + ρ s (C ps (Tm − Tav ) + λ )∂r r = RDt

計算例
ペレット平均温度

ペレット中心温度
ペレット径

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

出典：S.Taniguchi,M.Ohmi and
S.Ishiura, Trans. ISIJ, 571,23(1983)
出典：

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

非定常熱伝導方程式(円柱座標系)(溶融体の温度状態を決定する支配方程式)

固体ペレット／溶融体間の熱交換∂Tl  κ l ∂  ∂Tl n ∂Tl vθ ∂Tl22&ρ l C pl ++ vz
=r + ηγ + hsl (Ts − Tl )4πRr ∂θ∂z  r ∂r  ∂r V
 ∂t
熱対流

熱伝導

二軸スクリュの特徴に着目した仮定

vr << v z , vθ28

3φ (t )
hsl (Ts − Tl )R(t )ρ l : 密度

Tl : 温度C pl : 比熱t : 時間η : 粘度γ& : ひずみ速度

κ l : 熱伝導率

∂Tl ∂Tl∂T,<< l∂z ∂θ∂rCopyright© 2010- Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

粘性発熱 熱伝達二軸スクリュ内ペレット溶融状態の定量化例
0 MPa

0 MPa
Flow direction
Barrel radius:20mm
Screw speed : 100 rpm
Bulk density : 595 kg/m3
Initial pellet radius : 1.5 mm

Viscosity: Power law modelm0=5000,b=0.01,n=0.45
Physical properties

Case studies
Hsl (W/m2/K)

State

1

100

Filled2

1000

Filled3

2000

Filled4

1000

Un-filled(30kg/h)

Case

170℃

80℃
Barrel temperature condition

Copyright© 2010- Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved29溶融樹脂温度分布解析結果
30℃

180℃

Case1(hsl:100kW/m2/K)

Case2(hsl:1000kW/m2/K)

Case3(hsl:2000kW/m2/K)
Copyright© 2010- Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved30ペレット半径分布解析結果
0.0 mm

1.5 mmCase1(hsl:100kW/m2/K)

Case2(hsl:1000kW/m2/K)

Case3(hsl:2000kW/m2/K)
Copyright© 2010- Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved31

未充満＆溶融モデルCase4(hsl:1000kW/m2/K,30kg/h)解析結果
30℃

0 MPa

180℃

Temperature

3 MPa

Pressure
0 mm

1.5 mmPellet radius
Copyright© 2010- Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved32◇2.5D FEM解析技術の適用事例
繊維破断解析

繊維複合材の強度は、繊維長に依存して大きく変化する。一般的に繊維長の増加に伴って、剛性、強度、耐衝撃性能は何れも向上する。

出典：Copyright© 2010- Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved33

Fiber attrition model*)
繊維長を離散化モデルで表現：

li = i∆l (i = 1 ~ i max, ∆l : 繊維長刻み幅)

単位体積あたりに含まれる繊維長liの繊維数：

N i(i = 1 ~ i max)i max

dN i= − Pi N i + ∑ Rik N kdt
k =i +1
繊維長liの繊維数の変化率

繊維長liの繊維の破断に伴う減少率

別の繊維長(k>i)の破断に伴う繊維長liの繊維の増加率*)参考文献：J.H. Phelps (2009). “Processing-Microstructural Models for Short- and Long-Fiber ThermoplasticComposites,” PhD Thesis, University of Illinois at Urbana-Champaign, Urbana, IL 61801.

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

for γˆ& ≥ 1,

ζ：無次元抗力係数

Pi = 0 for γˆ& < 1,

ηm：流体粘度4&ςηγ4lγˆ& = 3 m i4π Ef d f

Ef :繊維ヤング率df：繊維断面直径

 ( x − µ )2 ,
Rij = C exp −22σlj
x = li , µ = , σ = Sl j2
Copyright© 2010- Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

Cb：無次元破断頻度係数

S:繊維長分布の調整パラメータC:規格化パラメータ35繊維数計算例

繊維長（μm)Copyright© 2010- Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

36Fiber number

二軸スクリュ内繊維長分布の定量化例
1000

600

900

500

800
700

400600

900

450

800

400

700

350

600

300

500

250300

500

500

400

200

400

200

300
200

345678910000210050

100
1200

100

1000300

150

1234567891102345678

910

12345

Fiber length segment
40 kg/h
100 rpmi max

6000μmInitial fiber length

重量平均繊維長分布 Lw =∑N Li =1i max∑N Li =1Copyright© 2010- Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

37i2iii

3000μm678910繊維長分布解析結果
0860

Segment 1 (Fiber length:0~600μm)
1000Segment 8 (Fiber length:4200~4800μm)
1000150Segment 10 (Fiber length:5400~6000μm)
Copyright© 2010- Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved38

二軸スクリュ押出機内充満状態の実験検証解析

プラスチック成形加工学会第27回年次大会発表資料（A-201,2016)

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

39

Flow direction

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

41

Copyright© 2010- Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved42◇まとめ成果
・2.5D FEM解析技術は、既往FAN法の強みである未充満状態や溶融状態の評価が可能であることに加えて、1次元次数の高いレベルで、各種解析結果の評価が可能である。
・当解析技術の運用は、既往FANソフトと同様、容易であり、計算負荷は、3DFEM/FVMと比較して圧倒的に小さい。
・当解析技術では、解析に際して負荷の大きいメッシュ生成作業を軽減するとともに、計算要素情報の任意編集を許容し、定型から外れた特殊エレメントの分析も守備範囲とする。
・連続体のトレース技術(非線形輸送方程式の解法）を併用することで、複合材の製造プロセスや繊維破断などの付加価値成形の分析に対応可能である。

今後の課題
・大学や企業（可視化装置）との連携を計り、解析結果の妥当性の検討や新規定量化モデルの構築などを目的とした共同研究（実験検証解析）を推進する。

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