General-purpose
Multi-Phase Dynamics Program

MUFFIN- MUltiFarious FIeld simulator for Non-equilibrium system -Tatsuya Yamaue
JCII, Doi projectDevelopers of MUFFIN – WG3 Doi Project Nagoya Univ.
- Takashi Taniguchi
JCII, Doi project- Tatsuya Yamaue (Japan Research Institute. Co.Ltd.)- Makoto Sasaki (Japan Research Institute. Co.Ltd.)- Naoki Kobayashi (Mitsui Chemical. Co.)
- Masahiko Noda (JSR. Co.)
Yamagata Univ.
- Akiyoshi Kuroda
Hiroshima Univ.
- Toru Okuzono

Time/Space Scales of MUFFIN
Time Scale

Multi-PhaseDynamics

0 sec-3 msec

GOURMETMUFFIN-6 μsec

SUSHI-9 nsec

PASTA-12 psec-15 fsec

COGNACSpace Scale-15fm

-12pm

-9nm

-6μm-3mm0m

Scope of Problems handled by MUFFIN
Aim:Design and Control ofMulti-Phase Structure of Polymeric Materials
Dispersion structures.
- Coagulation and breakup of droplets.
- Effect of surface modification.

Micro process control by pressure, shear, electric field, etc.
- Micro fabrication of thin films

Relation between stiffness and toughness, and morphology.
- Deformation of elastic materials.

Collaboration (zooming) with meso-simulators.
Electrolyte fluid kinetics near eletrodes.
Chemical reaction in narrow channels.
Swelling and deswelling of gels.

Features of MUFFIN
Various Continuum models◆ Various methods.
- Finite difference method (FDM) : Eulerian picture.
- Finite element method (FEM) : Eulerian and Lagrangian picture.
◆ Various Meshes.
- Structured mesh.
- Unstructured mesh (2D/3D). (Delaunay triangulation 2D/3D)◆ Various Dynamics.
- Free energy : Flory-Huggins, Electro-Chemical, Elastic, Gels,….
- Velocity : Navier-Stokes, Stokes, Oseen, Two fluid model,…..
- Concentration : Diffusion and flow of polymer, solvent, counter ions,….
- Electric potential : Poisson, Poisson-Boltzmann, ….
- Pressure………
- Polymer Stress………
- …….

T.Yamaue
T.Taniguchi

Architecture of MUFFIN

Object-oriented designing
for general multi-scale continuum fields simulation
SuperField. (base object of physical fields)- Equip various methods for physical fields. Message passing between different fields.

DynamicsManager.
- Control I/O, construct dynamics and execute simulation for given schedule.

FreeElectricVelocity Pressure Energy Potential…
SuperFieldMesh

Solver B.C. Result Parameter

Dynamics Manager

Application layer
Various simulators
MUFFIN library layer
API for various fields.

Core component
libraries (reusable).

UDF(XML) I/O

Component (Object) Structures of MUFFIN- “MUFFIN” isn’t a name of simulator but a name of architecture. -T.Yamaue
T.Taniguchi

Structure of I/O data sheets of MUFFIN
Parameter-Collaboration of different simulators.
-Proposal of a Front-end of multi-scale simulator.

Mesh data

Visualization

Plot

Partial Region data
Field data
Partial Region(Boundary) condition
Unit set on boundaryDynamics ManagerResult
Various unit sets andcurrent unit set.

cf)
XSIL (XML Scientific Interchange)
Prof. R.Williams, CACR, Caltech.
(Center for Advanced ComputingResearch)Simulation Packages of MUFFIN
Six simulation engines (three fluid simulators,
two solid simulators and one optics simulator) are now on release.
PhaseSeparation
Electrolyte
MEMFluid
Elastica
GelDyna
TURBAN

shear

PhaseSeparation : The multi-phase dynamics simulator forpolymeric fluids under the shear flow and the electric field with the GelDyna : The largehydrodynamic effect. (Fig. the phase separation dynamics of polymer deformation dynamics
blend under the shear (left), Spin coating simulation with evaporation simulator for multi-phase gelsystem by change of external
under shear (right))
stimuli, such as temperature,
pH, and electric field, andexternal forces.

TURBAN : The lighttransmittance simulator

Elastica : The linear elasticity simulator
for multi-phase solid system, including
both isotropic and anisotropic elasticity.
(Fig. the deformation (left) and strainenergy (right) of cylinder phase structure
obtained by SUSHI under the externalforce.)

Electrolyte : The

MEMFluid : The chemical
multi-component
electrolyte fluid
reaction and extraction
dynamics simulator for multidynamics simulator
under the shear flow and component electrolyte fluids
the electric field.
in the micro-fluidity chip.

Multi-Fluid Phase Dynamics Simulator“PhaseSeparation”

T.Taniguchi
M.Sasaki

The equations of PhaseSeparation
Basic model:∂ψ α= −∇ ⋅ (ψ α v ) + ∇ ⋅ ( Lα ∇µα )∂t∂v
t= −∇p + ∇ η ∇v + (∇v ) + K
ρ∂t∇⋅v = 0[{}]●Flow： incompressible, slow flow （Stokes flow∂ｖ /∂t = 0）
●Chemical potential μα , driving force K are derived from thefunctional of the free energy. Lα is the transport coefficient.

T.Taniguchi

Examples of the free energy functional●Flory-Huggins●Multi component fluid under an electric fieldT.Taniguchi

Applications of “PhaseSeparation FDM”
Phase separation
with hydrodynamics
effect.
•Flory-Huggins model.
2-dimensional, Periodic
B.C.

t=1

t=4

t=20
Coagulation of
2 droplets,
2-dimensional.

t=0

t=2

t=10

T.Taniguchi

Application of “PhaseSeparation FDM”
Phase separation under shear→→→Simulation of Spin Coating

M.Noda(An application of “PhaseSeparation FDM”)
Problems for spin coating– Surface roughness– Phase separated structure inside film
Experimental approaches– Selectivity to the solvent– Adhesion to the substrate– Viscosity
Proposed models

Macromolecules, 30, 4995 (1997)

Macromolecules, 29, 3232 (1996)

Macromolecules, 30, 4995 (1997)

Macromolecules, 29, 3232 (1996)

Polymer, 42, 1121 (2001)

M.NodaSpin Coating Simulation
Evaporation of solvent– evaporation rate
J. Electrochem. Soc., 138, 317 (1991)

E = κ ( x 0 1 − x1 ∞ )κ：mass transfer coefficientx01：initial solv. fraction
x1∞ ：equilibrium solv. fraction

Volume Fraction

1.0
0.8
0.6

Solvent
PolymeVoid0.4
0.2
0.001020
Step3040

M.NodaSpin Coating Simulation◆ Solvent Selectivity (csolv-poly)- Poly1 / Poly2 / solvent (better for Poly1)- PMMA / PS / MEK (better for PMMA)◆ Adhesion to the substrate (gS )- Poly1 / Poly2 /substrate (better for Poly2)- PMMA / PS / ODM (better for PS)◆ Evaporation rate (E )- related to spinning rate
◆ Viscosity (h )- effect of hydrodynamic flow

Electrolyte Fluid Dynamics Simulator“Electrolyte”● Simulate ion kinetics in electrolyte solutions.
● 3-dimensional FDM and FEM simulators.
● Calculate the electro-kinetic phenomena such as electroosmosis/electro-phoresis etc.)
The equations of Electrolyte●Time evolution of ion density●Flow field●Electric field

T.Taniguchi
M.SasakiT.Taniguchi

Applications of “Electrolyte FDM”
Ion density distribution

−

t=1

t=10＋

Ion distribution between electrodes
under electric field
Flow ofelectrolyte
around achargedparticle.

t=1t=10Velocity (on wall)

Application of “Electrolyte FEM”

M.Sasaki

Ion distribution around a charged spheremoving in an electrolyte solution＋＋＋＋＋＋＋＋ ＋＋＋＋＋（3-dimensional Delaunay mesh ）M.Sasaki

Application of “Electrolyte FEM” :Electro osmosis and Electrophoresis
Electrolyte flow induced by electric fieldaround a charged particle(simulation of electro-osmosis or electro-phoresis)

EEEE

T.Yamaue

Micro Fluidics Chip Simulator “MEMFluid”
Equations of motion.
Concentration of ions : ∂Cα(diffusion and reaction) ∂t = −∇ ⋅ ( vCα ) − ∇ ⋅ jα + ∑ R1αβ C β + ∑ R 2αβγ C β Cγββ ,γFlux density of ions : jα = − Lα [k BT {∇Cα + ∑ χ αβ Cα ∇C β } + eZα Cα (∇Φ − E 0 )][ {β}]

Oseen equation :t()
0=
−∇p+
∇η
∇
v+
∇v
+Kw
(laminar flow)
Volume force : K = −∑ k BT {∇Cα + ∑ χ αβ Cα ∇C β }αβ}(∇ ⋅ v = 0)
In Bulk(charge neutrality)

Laplace equation : ∇ 2 Φ = 0Cα : α-ion concentration jα : Flux density of α-ion Lα : Mobility of α-ionχαβ : Interaction parameter between α-ion and β-ion
Zα : Charge of α-ion
E0 : External electric fieldΦ : Electric potentialApplications.
MEMS (Micro Electro Mechanical System), Lab-on-a chip.
Micro reactor, μTAS (Total Analysis System), Bio chips.

T.Yamaue

Application of “MEMFluid”:Y-shaped micro-reactor
APressure field

C
40sec

60sec
16cases:
ΔP=1.0,2.5,5.0,10.0
R2=1e-3,1e-2,1e-1,1.0
Results:
For increase Pc.
R2 > 0.4:ΔP->10.0
R2 < 0.4:ΔP->1.0

80sec

100secProductivity (Pc) of C-ion
vs logR2 and logΔP

T.Yamaue

Application of “MEMFluid” :Electro-osmotic chip
Velocity on boundary is determined by electro-kinetics
Velocity at the boundary :εv eo =ςEηw
Helmholtz-Smouluchowski eq.
Calculation of ζ-potential on boundary:- 1-dim Poisson-Boltzmann eq. (1-1 electrolyte)
σσ22k BT1/ 2ς =ln[+(+1)]
1/ 2
Simulation result:
Ze(8Cεk BT )8Cεk BT- Condenser model (linear approximation)
valid for Ze | ς | / k BT << 1 ς ≤ 25( mV ){ς =σεκZα Cα e∑α2κ =(2)1 / 2εk BT
Scaling factor: Bejurem length/System size :lbl=e2( k B T εl )≈ 10 −3 − 10 − 2

ΔP=10.0, R2=1.0

Multi-Phase Linear Elasticity Simulator“Elastica”

M.Sasaki● Basic model： linear elasticity, 3-dimensional FEM

F {u i ( x )} = ∫ V d d x{ f } − ∫ V d d xρ ( x )g i u i ( x ) − ∫ St d d −1 xT i u i ( x )■ isotropic elastic material21 K ( x)(ell )2
f = G ( x ) eij − δ ij ell  +d2

1  ∂ui ∂u j eij =+2  ∂x j ∂x j 

■ anisotropic elastic material（axissymmetric）

f = D1 (eii ) + D2 (ni n j eij ) + D3(ell ⋅ ni n j eij ) + D4 nl eil ⋅ nk eik + D5 eij eij
2●2Arbitrary mixture and morphology of isotropic/anisotropic materials.

Stiffness Matrix for Anisotropic elasticity（axissymmetric）
0 0ll
σ xx  nσ   l k + m k − m 0 0
 yy  σ zz   l k − m k + m 0 0σ  = 00m 0
 yz  0
σ zx  0000 µ  000 0
σ xy  00  e xx 0  e yy 0   e zz  0  e yz 0   e zx  µ  exy f = D1 (eii ) + D2 (ni n j eij ) + D3(ell ⋅ ni n j eij ) + D4 nl eil ⋅ nk eik + D5 eij eij22

 D1 = (k − m ) 2
 D = (n + k − m ) 2 − l − µ 2 D3 = l − k + m
D = µ − m 4 D5 = m

M.SasakiApplication of “Elastica” : Elasticity of the two
T.Yamauephase structure calculated by SUSHI.
M.Sasaki【Before deformation】【Volume fraction】
Top:
apply load
・component1(green)z：hard, x,y ：soft・component2（red）
x,y：hard, z ：softｚ

Bottom: fixed【Strain energy】【after deformation】

Application of “Elastica”: Deformation of photo
M.Nodaregist pattern by capillary force

M.Sasaki

0.2µm

negative pressure induced by thecapillary force
1.0µm【before deformation】【after deformation】 Strain energy

Application of “Elastica”:Effective modulus oftwo phase system for various morphology
Total free energy and effective modulus

M.NodaM.Sasaki

Morphology (PP + Elastomer)21
1 2 F = G ∑  eij − δ ij ell   + K ∑  ekk d
2
 


Effective modulus is calculated byindependent deformation modes1 −δee∑  ij d ij ll  F
G +K =− 1 2
1 2e∑  2 kk ∑  2 ekk (a)dispersed(b)bi-continuous (C) by SUSHI

Calculated Young’s modulus

2modulus
spherebi-continuous
sushi1sushi2

E

3.753
34.67
91.84
50.99

analytic-model
seriesDaviessushi2: volume fraction is reset to 0 and 1
using a threshold value.

Gel and Elasticity Dynamics Simulator“GelDyna”

T.Yamaue

Equations of motion. (Stress-Diffusion Coupling Model for Gels)
Velocity of polymer :Velocity of solvent :incompressibility :ς (v p − v s ) = −φ ∇ ⋅ p + ∇ ⋅σς ( v s − v p ) = − (1 − φ ) ∇ ⋅ p∇ ⋅ [ φ v p + (1 − φ ) v s ] = 0v p :velocity of polymer

v s :velocity of solvent

p :pressureς :friction coefficientφ : volume fraction of polymer

σ : Cauchy stress tensor of polymer

Free energy of Gels.

φφk BT
1dF =d x [( 1 − φ ) ln( 1 − φ ) + χφ (1 − φ ) + ν 0( tr W + 2 ln)]
∫φ0φ0
v12

Stress of polymer networks

φσ ij = − [ φ f m ' ( φ ) − f m ( φ )] δ ij + ν 0( W ij − δ ij )φ0mixing free energy (Flory-Huggins): f m (φ ) = (1 − φ ) ln( 1 − φ ) + χφ (1 − φ )Wij: finger strain tensor

χ : polymer-solvent interaction

φ 0 : volume fraction of polymer in the relaxed stateν 0 : crosslink number density in the relaxed state

T.Yamaue

Applications of “GelDyna”
Applications.
Super-water absorbing gels
Temperature sensitive gels.
Drug delivery systems (DDS)
Actuators, sensors and switching devices.

Applications to swelling of NIPA gels.
2D

Bending of gels in rapid
deswelling process

3D

Light Transmittance Simulator“TURBAN”

N.Kobayashi=TURBidity ANalyzer“TURBAN” predicts thetransparency of polymeric material
having spherulite structures.
The image of the spherulite growth simulator.

“TURBAN” consists of two parts:(i) Calculate the spherulite growth
based on Avrami models.
(ii) Calculate the light transmittancebased on Maxwell equation.

The model figure about the light which
penetrates spherulite system.

N.Kobayashi“TURBAN” sample work
Light transmittance prediction of a PE film

Simulation conditions:- in the cubic cell of 1283
- a visible light wavelength = 400-700nm- non-uniform nucleation- using known properties of PE blown film
The turbidity spectra for system with thevarious average spherulite diameters.

Result:
The peak of turbidity is seen near thewavelength comparable to the meanspherulite diameter.
The peak of the turbidity is in the lower
range of a visible light.
These calculation are in good agreementwith experiments.

The turbidity integration value for thevisible light in each spherulite diameter.

T.Yamaue
A.Kuroda

Support tools of MUFFIN”Pre-processors”“MILK” (Mesh Generator)- Delaunay3D/2D auto meshing.
- Tetrahedron and Triangle.
- Internal Structures.

“NASTRAN BULK FILTER”- Import NASTRAN BULK mesh data to UDF mesh data.
FEM pre-processor（ex. HyperMesh）

GOURMETUDF convert“MeshFieldConvertor (IMPORT/EXPORT_...)”- Import/Export fields (morphology) data with SUSHI.

Support tools of MUFFIN”Post-processors”“MeshFieldShow (SHOW_...)”- Show mesh and field.
- Structured and Unstructured.
- Any cross section and Partial Regions.

“MeshFieldPlot (PLOT_...)”- Plot field (using gnuplot).
- Structured and Unstructured.
- Any cross section and Partial Regions.

“udf2avs”- Convert UDF field to AVS field.
- Structured and Unstructured.

Some Support tools for each package.
- “ModelingSupporter”, “analyze_reactor”, ……

T.Yamaue
T.Taniguchi
M.SasakiMUFFIN : Linking to other layers
Elastic modulus of
single components

Averaged elastic
modulus of multicompound

Diffusion constant of
single components
X-parameters

Surface modification

Surface tension
Stress Morphology
energyViscosity

PASTA
COGNACMUFFIN

SUSHI

Bulk
volumefraction

Strain rate
StrainMorphology

Summary
MUFFIN (MUltiFarious FIeld simulator for Non-equilibrium system)
is an architecture for simulations of spatio-temporal behavior of multiphase polymeric systems (〜0.1μm to 1mm, 〜sec).
Six simulation engines, “PhaseSeparation”, “Electrolyte”,“MEMFluid”, “Elastica”, “GelDyna” and “TURBAN” areimplemented.
MUFFIN is designed by the object-oriented modeling and serves
various reusable components libraries, such as structured andunstructured mesh, foundation of fields, I/O routines, matrix solvers,
finite elements solvers and “DynamicsManager”.
The MUFFIN component library will be useful for you to construct
your own simulator.