single mode Rolie-Poly model to calculate: continuous extension, and shear, step extension and shear /stress relaxation, and switch rates experiments.

Post date: Nov 1, 2013 6:49:16 PM

As in:

Likhtman, A. E. and R. S. Graham, "Simple constitutive equation for linear polymer melts derived from molecular theory: Rolie–Poly equation," J. Non-Newtonian Fluid Mech. 114, 1-12 (2003).

The following codes can calculate: continuous extension, and shear, step extension and shear /stress relaxation, and switch rates experiments.

Calculate by Matlab with ODE4.

function dx=RoliePoly(t,x)

dx = zeros(6,1);

global Tau_d;

global Shearr;

global Shearr2;

global t1stop;

global t2start;

global beta;

global delta;

global Tau_R;

global switchtime;

if t<=t1stop

K12=Shearr;

switchtime=t;

elseif (t<=(t2start+t1stop) )&& (t >t1stop)

K12=0;

elseif t> (t2start+t1stop)

K12=Shearr2;

end

dx(1) = 2*x(4)*K12+1/Tau_d+beta*( ( x(1)+x(2)+x(3) )/3 )^delta * 2*(1- sqrt( 3/(x(1)+x(2)+x(3)) ) )/Tau_R -(1/Tau_d + 2*(1- sqrt( 3/(x(1)+x(2)+x(3)) ) ) * ( beta*( ( x(1)+x(2)+x(3) )/3 )^delta +1)/Tau_R)*x(1);

dx(2) = 1/Tau_d+beta*( ( x(1)+x(2)+x(3) )/3 )^delta * 2*(1- sqrt( 3/(x(1)+x(2)+x(3)) ) )/Tau_R -(1/Tau_d + 2*(1- sqrt( 3/(x(1)+x(2)+x(3)) ) ) * ( beta*( ( x(1)+x(2)+x(3) )/3 )^delta +1)/Tau_R)*x(2);

dx(3) = 1/Tau_d+beta*( ( x(1)+x(2)+x(3) )/3 )^delta * 2*(1- sqrt( 3/(x(1)+x(2)+x(3)) ) )/Tau_R -(1/Tau_d + 2*(1- sqrt( 3/(x(1)+x(2)+x(3)) ) ) * ( beta*( ( x(1)+x(2)+x(3) )/3 )^delta +1)/Tau_R)*x(3);

dx(4) = x(2)*K12 -(1/Tau_d + 2*(1- sqrt( 3/(x(1)+x(2)+x(3)) ) ) * ( beta*( ( x(1)+x(2)+x(3) )/3 )^delta +1)/Tau_R)*x(4);

dx(5) = -(1/Tau_d + 2*(1- sqrt( 3/(x(1)+x(2)+x(3)) ) ) * ( beta*( ( x(1)+x(2)+x(3) )/3 )^delta +1)/Tau_R)*x(5);

dx(6) = x(5)*K12 -(1/Tau_d + 2*(1- sqrt( 3/(x(1)+x(2)+x(3)) ) ) * ( beta*( ( x(1)+x(2)+x(3) )/3 )^delta +1)/Tau_R)*x(6);

end

function dx=RoliePolyExt(t,x)

dx = zeros(2,1);

global Tau_d;

global Extenr;

global Extenr2;

global t1stop;

global t2start;

global beta;

global delta;

global Tau_R;

global switchtime;

if t<=t1stop

K33=Extenr;

switchtime=t;

elseif (t<=(t2start+t1stop) )&& (t >t1stop)

K33=0;

elseif t> (t2start+t1stop)

K33=Extenr2;

end

dx(1) = 1/Tau_d+beta*( ( 2*x(1)+x(2) )/3 )^delta * 2*(1- sqrt( 3/( 2*x(1)+x(2) ) ) )/Tau_R - ( K33+ 1/Tau_d + 2*(1- sqrt( 3/( 2*x(1)+x(2) ) ) ) * ( beta*( ( 2*x(1)+x(2) )/3 )^delta +1)/Tau_R)*x(1);

dx(2) = 1/Tau_d+beta*( ( 2*x(1)+x(2) )/3 )^delta * 2*(1- sqrt( 3/( 2*x(1)+x(2) ) ) )/Tau_R - (-2*K33+ 1/Tau_d + 2*(1- sqrt( 3/( 2*x(1)+x(2) ) ) ) * ( beta*( ( 2*x(1)+x(2) )/3 )^delta +1)/Tau_R)*x(2);

end

Main code begins here:

global Tau_d;

Tau_d=1340;

global Z;

Z=33;

global Tau_R;

Tau_R=Tau_d/3/Z;

global Shearr;

global stepStrain;

global beta;

global delta;

beta=0.1;

delta=-0.5;

global Shearr2;

global t1stop;

global t2start;

t1stop=0.04;

t2start=3000;

Shearr2=2;

global switchtime;

global Extenr;

global Extenr2;

%%----------------------------------

G0p=6800;

global Henckyr;

Henckyr=1;

Shearr2=285;

t1stop=1.5/Shearr;

t2start=0.2;

Z=33;

beta=1;

Tau_d=445.5;

tspan = logspace(-6,1,1000);

tspan = linspace(0,10,1000);

%%-----------------------------------------------

%%-step-strain extension relaxation Modulus 1.2 1.3 1.5 1.7lambda

%%-----------------------------

%%-----------------------------------------------

x0=[1;1];

tspan = logspace(-5,3.4,10000);

tspan=tspan';

stepSratio=[1.2,1.3,1.5,1.7];

stepStrain=log(stepSratio);

relaxModulus = zeros(4,1);

i=1;

Stress=zeros(10000,4);

relaxModulus=zeros(10000,4);

l2l=zeros(10000,4);

lambda=zeros(10000,4);

HenckyS=zeros(10000,4);

engStress=zeros(10000,4);

Extenr=10;

for i=1:4,

t1stop=stepStrain(i)/Extenr;

x=ode4(@RoliePolyExt,tspan,x0);

Stress(:,i)=x(:,2)-x(:,1);

HenckyS(:,i)=(tspan - ((tspan - t1stop) + abs (tspan- t1stop) )/ 2 )*Extenr;

lambda(:,i)=exp(HenckyS(:,i));

l2l(:,i)=lambda(:,i).*lambda(:,i)-1./lambda(:,i);

engStress(:,i)=Stress(:,i)./ lambda(:,i);

relaxModulus(:,i)=Stress(:,i)./ l2l(:,i);

end

figure(5)

loglog(tspan,relaxModulus(:,1),tspan,relaxModulus(:,2),tspan,relaxModulus(:,3),tspan,relaxModulus(:,4));

axis([0.01 3000 0.1 2])

loglog(tspan,Stress(:,1),tspan,Stress(:,2),tspan,Stress(:,3),tspan,Stress(:,4));

axis([0.01 3000 0.01 1.2])

figure(3);

loglog(tspan,relaxModulus(:,1),tspan,relaxModulus(:,2),tspan,relaxModulus(:,3),tspan,relaxModulus(:,4));

axis([0.01 3000 0.01 2])

%%-----------------------------

%%continuous extension different rates

%%-------------------------------------------------

Stress=zeros(10000,5);

relaxModulus=zeros(10000,5);

viscosity=zeros(10000,5);

engStress=zeros(10000,5);

x0=[1;1];

tspan = logspace(-5,3.4,10000);

conRate=[0.001, 0.003,0.01,0.03,0.1,0.3];

for i=1:6,

t1stop=3/conRate(i);

Extenr=conRate(i);

x=ode4(@RoliePolyExt,tspan,x0);

Stress(:,i)=x(:,2)-x(:,1);

relaxModulus(:,i)=Stress(:,i)./conRate(i);

end

loglog(tspan,Stress(:,1),tspan,Stress(:,2),tspan,Stress(:,3),tspan,Stress(:,4) ,tspan,Stress(:,4) );

axis([0.01 3000 0.01 1.2])

plot(tspan,Stress(:,1),tspan,Stress(:,2),tspan,Stress(:,3),tspan,Stress(:,4) ,tspan,Stress(:,4) );

axis([0.01 30 0 300])

for i=1:5,

t1stop=3/conRate(i);

viscosity(:,i)=Stress(:,i)/conRate(i);

lambda(:,i)=(tspan - ((tspan - t1stop) + abs (tspan- t1stop) )/ 2 )*conRate(i);

engStress(:,i)=Stress(:,i)./ lambda(:,i);

end

loglog(tspan,viscosity(:,1),tspan,viscosity(:,2),tspan,viscosity(:,3),tspan,viscosity(:,4) ,tspan,viscosity(:,4) ,tspan,viscosity(:,5),tspan,viscosity(:,6)) ;

loglog(tspan,engStress(:,1),tspan,engStress(:,2),tspan,engStress(:,3),tspan,engStress(:,4) ,tspan,engStress(:,4) ,tspan,engStress(:,5),tspan,engStress(:,6)) ;

axis([0.1 3000 1 10])

loglog(lambda(:,1),engStress(:,1),lambda(:,2),engStress(:,2),lambda(:,3),engStress(:,3),lambda(:,4),engStress(:,4),lambda(:,5),engStress(:,5),lambda(:,6),engStress(:,6)) ;

axis([0.1 30 1 10])

t1stop=3000;

for i=1:6,

Extenr=conRate(i);

x=ode4(@RoliePolyExt,tspan,x0);

Stress(:,i)=x(:,2)-x(:,1);

viscosity(:,i)=Stress(:,i)/conRate(i);

HenckyS(:,i)=(tspan - ((tspan - t1stop) + abs (tspan- t1stop) )/ 2 )*conRate(i);

engStress(:,i)=Stress(:,i)./ lambda(:,i);

relaxModulus(:,i)=Stress(:,i)./ log(lambda(:,i));

lambda(:,i)=exp(HenckyS(:,i));

l2l(:,i)=lambda(:,i).*lambda(:,i)-1./lambda(:,i);

end

loglog(tspan,viscosity(:,1),tspan,viscosity(:,2),tspan,viscosity(:,3),tspan,viscosity(:,4) ,tspan,viscosity(:,4) ,tspan,viscosity(:,5),tspan,viscosity(:,6)) ;

axis([0.1 3000 0.3 3000])

figure(2);

loglog(tspan,engStress(:,1),tspan,engStress(:,2),tspan,engStress(:,3),tspan,engStress(:,4) ,tspan,engStress(:,4) ,tspan,engStress(:,5),tspan,engStress(:,6)) ;

axis([0.1 3000 1 10])

figure(3);

loglog(lambda(:,1),engStress(:,1),lambda(:,2),engStress(:,2),lambda(:,3),engStress(:,3),lambda(:,4),engStress(:,4),lambda(:,5),engStress(:,5),lambda(:,6),engStress(:,6)) ;

axis([0.1 30 1 10])

figure(4)

loglog(l2l(:,1),engStress(:,1),l2l(:,2),engStress(:,2),l2l(:,3),engStress(:,3),l2l(:,4),engStress(:,4),l2l(:,5),engStress(:,5),l2l(:,6),engStress(:,6)) ;

axis([0.1 30 1 10])

plot(l2l(:,1),engStress(:,1),l2l(:,2),engStress(:,2),l2l(:,3),engStress(:,3),l2l(:,4),engStress(:,4),l2l(:,5),engStress(:,5),l2l(:,6),engStress(:,6)) ;

axis([0 20 1 10])

%%--------------------------------

%%---************Extension switch rates

%%******************************************************

%%**************************************************

x0=[1;1];

tspan = logspace(-5,3.95,100000);

tspan1 = logspace(-5,2,5000);

tspan2 = linspace(100.3,9500,10000);

tspan =[tspan1,tspan2];

tspan=tspan';

Stress=zeros(15000,11);

engStress=zeros(15000,11);

HenckyS=zeros(15000,11);

lambda=zeros(15000,11);

t2startwitch=[0.9,5,10,20,50,100,500, 1500, 3200, 4800,7500];

Extenr=1;

Extenr2=0.03;

t1stop=log(2.2)/Extenr;

i=1;

for i=1:11,

t2start=t2startwitch(i);

x=ode4(@RoliePolyExt,tspan,x0);

Stress(:,i)=x(:,2)-x(:,1);

HenckyS(:,i)=(tspan - ((tspan - t1stop) + abs (tspan- t1stop) )/ 2 )*Extenr + ((tspan-t2start-t1stop)+ abs (tspan-t2start-t1stop) )/2*Extenr2 ;

lambda(:,i)=exp(HenckyS(:,i));

engStress(:,i)=Stress(:,i)./ lambda(:,i);

end

figure(5);

plot(tspan,HenckyS(:,1));

axis([0 30 0 1])

figure(6);

plot(tspan,lambda(:,1));

axis([0 30 0 1])

figure(7);

semilogx(tspan,engStress(:,1));

axis([0.1 3000 0 2])

figure(4);

semilogx(tspan,engStress(:,1),tspan,engStress(:,2),tspan,engStress(:,3),tspan,engStress(:,4) ,tspan,engStress(:,4) ,tspan,engStress(:,5),tspan,engStress(:,6), tspan,engStress(:,7),tspan,engStress(:,8),tspan,engStress(:,9),tspan,engStress(:,10) ,tspan,engStress(:,11) );

axis([0.1 10000 0 1.8])

for i=1:11,

jj=1;

for jj=1:15000,

if tspan(jj)> (t1stop+t2startwitch(i))

relaxbeforestress(i)=engStress(jj-1,11);

[engstressmax(i),rowmax(i)]=max(engStress(jj-1:15000,i));

rowmax(i)=rowmax(i)+jj-2;

relaxingstress(i)=engStress(rowmax(i),11);

break;

end

end

end

engstressmax=engstressmax';

relaxingstress=relaxingstress';

relaxbeforestress=relaxbeforestress';

%%--Shear

%%------------------------------------------

tspan = logspace(-5,3.4,10000);

x0=[1;1;1;0;0;0];

tspan = logspace(-3,3.4,1000);

tspan=tspan';

[t,x]=ode45(@RoliePoly,tspan,x0);

plot(tspan, x(:,4));

loglog(tspan, x(:,4));

axis([0.1 10 0.01 1.2])

%%---------------------------

%%step shear, relaxation G

%%-------------------

t1stop=0.04;

t2start=3000;

Shearr2=2;

x0=[1;1;1;0;0;0];

tspan = logspace(-5,3.4,10000);

tspan=tspan';

stepStrain=[0.1,0.6,0.7,0.8,0.9,1,1.1];

Stress=zeros(10000,7);

ShearStrain=zeros(10000,7);

relaxModulus=zeros(10000,7);

for i=1:7,

Shearr=stepStrain(i)/0.04;

x=ode4(@RoliePoly,tspan,x0);

Stress(:,i)=x(:,4);

ShearStrain(:,i)=(tspan - ((tspan - t1stop) + abs (tspan- t1stop) )/ 2 )*Shearr;

relaxModulus(:,i)=Stress(:,i)./ShearStrain(:,i);

end

loglog(tspan,Stress(:,1),tspan,Stress(:,2),tspan,Stress(:,3),tspan,Stress(:,4),tspan,Stress(:,5),tspan,Stress(:,6),tspan,Stress(:,7) );

axis([0.01 3000 0.01 1.2])

loglog(tspan,relaxModulus(:,1),tspan,relaxModulus(:,2),tspan,relaxModulus(:,3),tspan,relaxModulus(:,4),tspan,relaxModulus(:,5),tspan,relaxModulus(:,6),tspan,relaxModulus(:,7) );

axis([0.01 3000 0.1 1.2])

axis([0.00001 3000 0.1 1.2])

figure(2);

plot(tspan,relaxModulus(:,1),tspan,relaxModulus(:,2),tspan,relaxModulus(:,3),tspan,relaxModulus(:,4),tspan,relaxModulus(:,5),tspan,relaxModulus(:,6),tspan,relaxModulus(:,7) );

axis([0.037 0.043 0.95 0.98])

figure(3);

loglog(tspan,relaxModulus(:,1),tspan,relaxModulus(:,2),tspan,relaxModulus(:,3),tspan,relaxModulus(:,4),tspan,relaxModulus(:,5),tspan,relaxModulus(:,6),tspan,relaxModulus(:,7) );

axis([0.1 3000 0.1 1.1])

---------------------------------------------

loglog(tspan,Stress(:,1),tspan,Stress(:,2),tspan,Stress(:,3),tspan,Stress(:,4),tspan,Stress(:,5),tspan,Stress(:,6),tspan,Stress(:,7),tspan,Stress(:,8),tspan,Stress(:,9),tspan,Stress(:,10),tspan,Stress(:,11),tspan,Stress(:,12),tspan,Stress(:,13),tspan,Stress(:,14),tspan,Stress(:,15),tspan,Stress(:,16),tspan,Stress(:,17),tspan,Stress(:,18),tspan,Stress(:,19),tspan,Stress(:,20),tspan,Stress(:,21),tspan,Stress(:,22));

axis([0.1 10 0.01 3])

axis([0.1 10 1 3])

axis([0.1 10 0.1 3])

[stressmax,rowmax]=max(Stress);

tmax=(tspan(rowmax))';

loglog(tspan,Stress(:,1),tspan,Stress(:,2),tspan,Stress(:,3),tspan,Stress(:,4),tspan,Stress(:,5),tspan,Stress(:,6),tspan,Stress(:,7),tspan,Stress(:,8),tspan,Stress(:,9),tspan,Stress(:,10),tspan,Stress(:,11),tspan,Stress(:,12),tspan,Stress(:,13),tspan,Stress(:,14),tspan,Stress(:,15),tspan,Stress(:,16),tspan,Stress(:,17),tspan,Stress(:,18),tspan,Stress(:,19),tspan,Stress(:,20),tspan,Stress(:,21),tspan,Stress(:,22),tmax,stressmax);

axis([0.1 10 0.1 3])

Shearr3=Shearr;

Shearr=Shearr2;

t1stop=10;

[t,x]=ode45(@RoliePoly,tspan,x0);

loglog(tspan, x(:,4));

axis([0.0001 10 0.01 4])