Pom-Pom model, and with drag-strain coupling, by Matlab

Ref:

McLeish, T. and R. Larson, "Molecular constitutive equations for a class of branched polymers: The pom-pom polymer," J. Rheol. 42, 81-110 (1998).

Blackwell, R. J., T. C. B. McLeish and O. G. Harlen, "Molecular drag-strain coupling in branched polymer melts," J. Rheol. 44, 121-136 (2000).

pom-Drag-coupling-shear

function dx=pomDragCoupleshear(t,x)

global Shearr;

global tau_b;

global tau_s;

global q;

global v_star;

K12=Shearr;

dx=[

K12*x(2)+x(4)*K12-(1/tau_b)*(x(1)-1/3);

K12*x(5)-(1/tau_b)*(x(4));

K12*x(8)-(1/tau_b)*(x(7));

K12*x(5)-(1/tau_b)*(x(2));

-(1/tau_b)*(x(5)-1/3);

-(1/tau_b)*(x(8));

x(6)*K12-(1/tau_b)*(x(3));

-(1/tau_b)*(x(6));

-(1/tau_b)*(x(9)-1/3);

(K12*x(4))/(x(1)+x(5)+x(9))*x(10)-(x(10)-1)/tau_s*exp(v_star*(x(10)-1));

];

end

pom-Drag-coupling-ext

function dx=pomDragCoupleext(t,x)

global Henckyr;

global tau_b;

global tau_s;

global q;

global v_star;

k(1)=-Henckyr/2;

k(2)=-Henckyr/2;

k(3)=Henckyr;

dx=[2*x(1)*k(1)-(x(1)-1/3)/tau_b;

2*x(2)*k(2)-(x(2)-1/3)/tau_b;

2*x(3)*k(3)-(x(3)-1/3)/tau_b;

(k(1)*x(1)+k(2)*x(2)+k(3)*x(3))/(x(1)+x(2)+x(3))*x(4)-(x(4)-1)/tau_s*exp(v_star*(x(4)-1));

];

end

pom-shear

function dx=pomshear(t,x)

global Shearr;

global tau_b;

global tau_s;

global q;

K12=Shearr;

if x(10)>q

x(10)=q;

end

dx=[

K12*x(2)+x(4)*K12-(1/tau_b)*(x(1)-1/3);

K12*x(5)-(1/tau_b)*(x(4));

K12*x(8)-(1/tau_b)*(x(7));

K12*x(5)-(1/tau_b)*(x(2));

-(1/tau_b)*(x(5)-1/3);

-(1/tau_b)*(x(8));

x(6)*K12-(1/tau_b)*(x(3));

-(1/tau_b)*(x(6));

-(1/tau_b)*(x(9)-1/3);

(K12*x(4))/(x(1)+x(5)+x(9))*x(10)-(x(10)-1)/tau_s;

];

end

pom-ext

function dx=pomext(t,x)

global Henckyr;

global tau_b;

global tau_s;

global q;

k(1)=-Henckyr/2;

k(2)=-Henckyr/2;

k(3)=Henckyr;

if x(4)>q

x(4)=q;

end

dx=[2*x(1)*k(1)-(x(1)-1/3)/tau_b;

2*x(2)*k(2)-(x(2)-1/3)/tau_b;

2*x(3)*k(3)-(x(3)-1/3)/tau_b;

(k(1)*x(1)+k(2)*x(2)+k(3)*x(3))/(x(1)+x(2)+x(3))*x(4)-(x(4)-1)/tau_s;

];

end

main body:

x0=[1/3;1/3;1/3;1];

tspan = logspace(-5,5,9000);

tspan=tspan';

q=14;s_a=1;

s_b=28.7;

tau_a=0.03;

phy=s_b/(2*q*s_a+s_b);

global tau_b;

tau_b=4/(pi)^2*s_b^2*phy*tau_a*q;

global tau_s;

tau_s=s_b*tau_a*q;

G0p=1560;

A0=[1/3;0;0;0;1/3;0;0;0;1/3;1];

global Shearr;

for i=-12:9,

Shearr=10^(i/3);

shearrate(i+13)=Shearr;

[t,A]=ode45(@pomshear,tspan,A0);

a = find(A(:,10)>q);

A(a,10)=q;

StreS(:,i+13)=15/4*G0p*phy^2*A(:,10).*A(:,10).*A(:,4)./(A(:,1)+A(:,5)+A(:,9));

viscS(:,i+13)=StreS(:,i+13)/Shearr;

S(:,i+13)=A(:,4)./(A(:,1)+A(:,5)+A(:,9));

lamd(:,i+13)=A(:,10);

end

[stressmax,rowmax]=max(StreS);

tmax=(tspan(rowmax))';

strainmax=tmax.*shearrate;

plotyy(shearrate,strainmax,shearrate,stressmax,'loglog')

shearrate=shearrate';strainmax=strainmax';stressmax=stressmax';

tmax=tmax';

shearrate=shearrate';strainmax=strainmax';stressmax=stressmax';

tmax=tmax';

beta=2;

global v_star;

v_star=2/(q-1);

for i=-12:9,

Shearr=10^(i/3);

shearrate(i+13)=Shearr;

[t,A]=ode45(@pomDragCoupleshear,tspan,A0);

StreS(:,i+13)=15/4*G0p*phy^beta*A(:,10).*A(:,10).*A(:,4)./(A(:,1)+A(:,5)+A(:,9));

viscS(:,i+13)=StreS(:,i+13)/Shearr;

end