Hello everyone,
I am working on a simulation of non-isothermal viscoelastic multiphase flow, and I am trying to implement the log-conformation formulation of the Oldroyd-B model. However, I am facing errors related to solver convergence, and I am not sure whether my implementation is correct. I am particularly confused about whether I have properly implemented the log-conformation transformation and its evolution equation. I have attached my code below. I would really appreciate it if someone could review it and point out any mistakes or suggest improvements. // ==============================================================================
// NON_ISOTHERMAL VISCOELASTIC LOG CONFORMATION OLDROYED B MULTIPHASE FLOW
// ==============================================================================
load “msh3”
load “iovtk”
// GEOMETRY (UNCHANGED) =========================================================
real Db=2.0, Dc=0.4, Lb=5.0, Ln=1.5, angle=90.0;
real Lair=5.0, Wair=2.0, Hb=Db/2.0, Hc=Dc/2.0;
real Lc=(angle>=89.9)?0.0:(Hb-Hc)/tan(angle*pi/180.0);
real Y0=0, Y1=-Lb, Y2=Y1-Lc, Y3=Y2-Ln, Y4=Y3-Lair;
int n=15;
border Inlet(t=Hb,-Hb){x=t;y=Y0;label=1;}
border Barrel(t=Y0,Y1){x=-Hb;y=t;label=2;}
border ContraL(t=0,1){x=-Hb+t*(Hb-Hc);y=Y1-t*Lc;label=3;}
border NozzleL(t=Y2,Y3){x=-Hc;y=t;label=6;}
border Step(t=-Hc,-Wair/2){x=t;y=Y3;label=4;}
border AirSide(t=Y3,Y4){x=-Wair/2;y=t;label=4;}
border AirBot(t=-Wair/2,Wair/2){x=t;y=Y4;label=5;}
border AirR(t=Y4,Y3){x=Wair/2;y=t;label=4;}
border StepR(t=Wair/2,Hc){x=t;y=Y3;label=4;}
border NozzleR(t=Y3,Y2){x=Hc;y=t;label=6;}
border ContraR(t=0,1){x=Hc+t*(Hb-Hc);y=Y2+t*Lc;label=3;}
border BarrelR(t=Y1,Y0){x=Hb;y=t;label=2;}
mesh Th=buildmesh(Inlet(n)+Barrel(2n)+ContraL(n)+NozzleL(2n)+Step(n)+AirSide(3n)
+AirBot(3n)+AirR(3n)+StepR(n)+NozzleR(2n)+ContraR(n)+BarrelR(2*n)
);
// PARAMETERS ===============================================================
real Re=0.01, Pr=50.0, Br=0.1, Fr=1.0, Eo=2.0;
real Wi0=0.5;
real beta=1./9., alpha=1-beta;
real C1=15, C2s=50;
real sigma=1.0;
real dt=0.001;
int Nt=300;
// CREATE SPACES =============================================================
fespace Vh(Th,P2);
fespace Qh(Th,P1);
Vh u,v,uold,vold,wu,wv;
Qh p,q;
Qh Temp,Told;
Qh phi,phiold;
Qh s11,s12,s22;
Qh sold11,sold12,sold22;
Qh G11,G12,G22;
// FUNCTIONS HAVISIDE AND WLF =================================================
func real H(real ph){return 0.5*(1+tanh(ph/0.02));}
func real WLF(real T){return exp((-C1*T)/(C2s+T+1e-8));}
// INITIAL ==================================================================
phi = (y > Y0 - 0.2) ? 1 : -1;
u=0; v=0; Temp=0; s11=0; s12=0; s22=0;
// TIME LOOP ==================================================================
for(int it=0; it<Nt; it++)
{
uold=u; vold=v; Told=Temp; phiold=phi;
sold11=s11; sold12=s12; sold22=s22;
// NON_ISOTHERMAL AND PHASE PROPERTIES =====================================
Qh ftheta=WLF(Told);
Qh etas=beta*ftheta*H(phiold)+(1-H(phiold))*0.01;
Qh etap=alpha*ftheta*H(phiold);
Qh Wi=Wi0*ftheta;
Qh rho=H(phiold)+0.1*(1-H(phiold));
// DEVSS ===================================================================
solve G11p(G11,q)=int2d(Th)(G11*q)-int2d(Th)(dx(uold)*q);
solve G22p(G22,q)=int2d(Th)(G22*q)-int2d(Th)(dy(vold)*q);
solve G12p(G12,q)=int2d(Th)(G12*q)-int2d(Th)(0.5*(dy(uold)+dx(vold))*q);
// LOG-CONFORMATION → A ======================================================
Qh trS=0.5*(sold11+sold22);
Qh delS=sqrt(max(0.25*(sold11-sold22)^2+sold12^2,1e-12));
Qh e1=exp(trS+delS), e2=exp(trS-delS);
Qh A11=0.5*(e1+e2)+0.5*(e1-e2)*(sold11-sold22)/(2*delS);
Qh A22=0.5*(e1+e2)-0.5*(e1-e2)*(sold11-sold22)/(2*delS);
Qh A12=0.5*(e1-e2)*sold12/delS;
Qh ArrInv=0.5*(1./e1+1./e2)+0.5*(1./e1-1./e2)*(sold11-sold22)/(2*delS);
Qh AzzInv=0.5*(1./e1+1./e2)-0.5*(1./e1-1./e2)*(sold11-sold22)/(2*delS);
Qh ArzInv=0.5*(1./e1-1./e2)*sold12/delS;
// CURVATURE + SURFACE TENSION ================================================
Qh phix=dx(phiold), phiy=dy(phiold);
Qh norm=sqrt(phix^2+phiy^2+1e-10);
Qh nx=phix/norm;
Qh ny=phiy/norm;
Qh kappa=dx(nx)+dy(ny);
Qh Fx = -sigma * kappa * phix;
Qh Fy = -sigma * kappa * phiy;
// NAVIER-STOKES ==============================================================
real tauDEV=0.1;
solve NS([u,v,p],[wu,wv,q])= int2d(Th)( rho*(u*wu+v*wv)/dt
+(2*etas/Re)*(dx(u)*dx(wu)+dy(u)*dy(wu) +dx(v)*dx(wv)+dy(v)*dy(wv))
+tauDEV*(dx(u)*dx(wu)+dy(u)*dy(wu)+dx(v)*dx(wv)+dy(v)*dy(wv))
-p*(dx(wu)+dy(wv)) -q*(dx(u)+dy(v))
)
-int2d(Th)( rho*(convect([uold,vold],-dt,uold)*wu +convect([uold,vold],-dt,vold)*wv)/dt
+tauDEV*(G11*dx(wu)+G12*dy(wu)+G12*dx(wv)+G22*dy(wv))
+(etap/(Re*max(Wi,1e-6)))*((A11-1)*dx(wu)+A12*dy(wu)+A12*dx(wv)+(A22-1)*dy(wv))
+Fx*wu + Fy*wv -(rho/(Fr^2))*wv
)
+on(1,u=0,v=1) +on(2,3,4,6,u=0,v=0) +on(5,p=0);
// LOG-CONFORMATION =========================================================
Qh srr=dx(u), szz=dy(v), srz=0.5*(dy(u)+dx(v)), omg=0.5*(dx(v)-dy(u));
solve S11(s11,q)=int2d(Th)(s11*q/dt)-int2d(Th)(convect([u,v],-dt,sold11)*q/dt
+(2*srr+2*omg*sold12-(1./max(Wi,1e-6))*(1-ArrInv))*q);
solve S22(s22,q)= int2d(Th)(s22*q/dt)-int2d(Th)(convect([u,v],-dt,sold22)*q/dt
+(2*szz-2*omg*sold12-(1./max(Wi,1e-6))*(1-AzzInv))*q);
solve S12(s12,q)=int2d(Th)(s12*q/dt) -int2d(Th)(convect([u,v],-dt,sold12)*q/dt
+(2*srz+omg*(sold22-sold11)-(1./max(Wi,1e-6))*(-ArzInv))*q);
// ENERGY ====================================================================
Qh d11=dx(u), d22=dy(v), d12=0.5*(dy(u)+dx(v));
Qh DD=d11*d11+d22*d22+2*d12*d12;
Qh AID=(A11-1)*d11+(A22-1)*d22+2*A12*d12;
Qh Diss=(Br/(Re*Pr))*(2*etas*DD+(etap/max(Wi,1e-6))*AID);
solve Energy(Temp,q)= int2d(Th)(rho* Temp*q/dt+(1/(Re*Pr))*(dx(Temp)*dx(q)+dy(Temp)*dy(q)))
-int2d(Th)(convect([uold,vold],-dt,Told)*q/dt+Diss*q);
// LEVEL SET ===================================================================
phi=convect([uold,vold],-dt,phiold);
// OUTPUT ======================================================================
if(it%5==0)
// savevtk("out"+it+".vtk",Th,u,v,p,Temp,phi);
plot(phi, fill=1, value=1, WindowIndex =0, cmm="Phase: "+it);
plot(s11, fill=1, value=1, WindowIndex =1, cmm="stress: "+it);
plot(Temp, fill=1, value=1, WindowIndex =2, cmm="isotherm: "+it);
cout<<"Step "<<it<<endl;
}
cout<<“Simulation complete”<<endl;