动态接触角

#includeudf.h#includesg_mphase.h#includesg_vof.h#includesg.h#includemem.h#defineVISCOSITY0.001#defineSURF_TENS0.0728#defineMYTRUE1#defineMYFALSE0#defineHoff(x)acos(1-(2.0*tanh(5.16*(pow((x/(1+(1.31*pow(x,0.99)))),0.706)))))#definestatic_Con_Ang60.#defineindex_source2/*ThisCodecomputesthenormalsoftheVOFfunction*/DEFINE_ADJUST(store_gradient,domain){Thread*t;Thread**pt;Thread*phaset;cell_tc;intphase_domain_index=1;/*SecondaryDomain*/Domain*pDomain=DOMAIN_SUB_DOMAIN(domain,phase_domain_index);voidcalc_source();Alloc_Storage_Vars(pDomain,SV_VOF_RG,SV_VOF_G,SV_NULL);Scalar_Reconstruction(pDomain,SV_VOF,-1,SV_VOF_RG,NULL);Scalar_Derivatives(pDomain,SV_VOF,-1,SV_VOF_G,SV_VOF_RG,Vof_Deriv_Accumulate);mp_thread_loop_c(t,domain,pt)if(FLUID_THREAD_P(t)){Thread*ppt=pt[phase_domain_index];begin_c_loop(c,t){calc_source(c,t);}end_c_loop(c,t)}Free_Storage_Vars(pDomain,SV_VOF_RG,SV_VOF_G,SV_NULL);}voidcalc_source(cell_tcell,Thread*thread){realVOF_Val[2],VOF_Mag,source,VOF_Norm[2];Thread*phaset;phaset=THREAD_SUB_THREAD(thread,1);if(C_VOF(cell,phaset)!=0.0&&N_TIME1){/*ThegradientsoftheVOFfunctionarefoundinthex,yandzdir.*/if(NULLP(THREAD_STORAGE(phaset,SV_VOF_G))){Message0(N_TIME=%d,....show-grad:GradientofVOFisnotavailable\n,N_TIME);Error(0);}VOF_Val[0]=-C_VOF_G(cell,phaset)[0];VOF_Val[1]=-C_VOF_G(cell,phaset)[1];/*ThemagnitudeoftheVOFgradientsisfoundsoitcanbenormalized*/VOF_Mag=NV_MAG(VOF_Val);if(VOF_Mag!=0.0){VOF_Mag=NV_MAG(VOF_Val);VOF_Norm[0]=VOF_Val[0]/VOF_Mag;VOF_Norm[1]=VOF_Val[1]/VOF_Mag;}else{/*Thisistoavoidthedividebyzerofunction*/VOF_Norm[0]=0.0;VOF_Norm[1]=0.0;}C_UDMI(cell,thread,0)=VOF_Norm[0];C_UDMI(cell,thread,1)=VOF_Norm[1];}source=0.0;C_UDMI(cell,thread,index_source)=source;}DEFINE_SOURCE(VOF_Norms,cell,thread,dS,eqn){realsource;source=C_UDMI(cell,thread,index_source);returnsource;}DEFINE_PROFILE(con_ang,t,pos){/*Firstthevariouspointervariablesarecreated*/face_tf;cell_tc;realfeta_d,vel_Val[2],cont_Line_Vel,VOF_Normal[2],cap_Num,static_Con_Rad,x_Bottom,x_Top,x_Bisect,hoff_Old,hoff_Cur,hoff_New,finish_Cond,inv_Hoff=0.0;intnotConverged,itNum;Thread*t0,*pt;/*Thiscodeisdesignedtofindthezerofortheinvertedhoffmanfunctionbyfindingthezero*//*ofthefunctionatwhichthehoffmanfunctionresultsinthestaticcontactangle*//*Firstthevariablesareassigned*/notConverged=MYTRUE;x_Bottom=0.001;x_Top=2.0;itNum=0;static_Con_Rad=((static_Con_Ang*M_PI)/180.);/*Awhileloopperformsthebisectionmethod,asimplebutverystablezerofinder*/while(notConverged){/*Thevariablesusedinthebisectionmethodareassignedandthehoffman*//*functionsareevaluated*/itNum++;hoff_Old=(Hoff(x_Bottom)-static_Con_Rad);hoff_Cur=(Hoff(x_Top)-static_Con_Rad);x_Bisect=(x_Bottom+x_Top)/2.0;hoff_New=(Hoff(x_Bisect)-static_Con_Rad);finish_Cond=fabs(1-(x_Bisect/x_Top));/*Theloopendswhentherelativeerrorislessthan1e-8andtheinverse*//*hoffmanvalueisstoredforuselater*/if(finish_Cond0.00000001||itNum10000000){inv_Hoff=x_Bisect;notConverged=MYFALSE;}/*Conditionsforthebisectionmethod*/if((hoff_Old*hoff_New)0.0){x_Top=x_Bisect;}if((hoff_Cur*hoff_New)=0.0){x_Bottom=x_Bisect;}}/*Nowthemainloopgoesthroughallthefacesintheboundary*/begin_f_loop(f,t){/*Thecellandphasethreadsareisolated*/c=F_C0(f,t);t0=THREAD_T0(t);pt=THREAD_SUB_THREAD(t0,1);/*ThemainformulationisonlyappliediftheVOFis0*/if(C_VOF(c,pt)!=0.0){/*Thevelocitiesarerecordedineachdirection*/vel_Val[0]=C_U(c,t0);vel_Val[1]=C_V(c,t0);/*TheVOFnormalsarebroughtin*/VOF_Normal[0]=C_UDMI(c,t0,0);VOF_Normal[1]=C_UDMI(c,t0,1);/*Thecontactlinevel.iscalcfromthedotproductofVOFandVel*/cont_Line_Vel=NV_DOT(vel_Val,VOF_Normal);/*Thecapillarynumberisfoundbasedoncontlinevel.*/cap_Num=fabs((VISCOSITY*cont_Line_Vel)/SURF_TENS);/*Thedynamiccontactangleisdefinedthenstoredintheprofile*/if(cap_Num+inv_Hoff0.0){cap_Num=inv_Hoff;}feta_d=((Hoff(cap_Num+inv_Hoff))*180)/M_PI;F_PROFILE(f,t,pos)=feta_d;}else{F_PROFILE(f,t,pos)=static_Con_Ang;}}end_f_loop(f,t)}