In[]:=
patternproperties=μ,A,γ,B-,λ;params={Dm0.03,Dc1,k1,R02,A0.001,M4,v0,tend4000000,throughPos0.65,Muldenbreite3,Muldensteigung1.5,α0.4,L10,Lminus2.2};f[m_,c_]:=c*m^2-k*m;mstat[x_]:=(2μ+A*(-Tanh[(x-(α*L-Lminus))/λ]+Tanh[(x-(α*L+Lminus))/λ]))/.patternpropertiescstat[x_]:=(γ/2+B*(-Tanh[(x-(α*L-Lminus))/λ]+Tanh[(x-(α*L+Lminus))/λ]))/.patternproperties
Dc
k
R0
2
Dm
Dc
k
R0
2
Dm
2
Dm
k
Dc
R0
Dm
k
Dc
R0
2
Dm
k
In[]:=
R1[X_]:=-0.5*A/(Muldenbreite^2)*(Tanh[(X-(throughPos*L-0.5*Muldenbreite))*Muldensteigung]-Tanh[(X-(throughPos*L+0.5*Muldenbreite))*Muldensteigung])*(Muldenbreite^2-(X-throughPos*L)^2)R1firstderiv[x_]:=-((0.5A(Muldenbreite^2-(-LthroughPos+x)^2)(-MuldensteigungSech[Muldensteigung(-0.5Muldenbreite-LthroughPos+x)]^2+MuldensteigungSech[Muldensteigung(0.5Muldenbreite-LthroughPos+x)]^2))/Muldenbreite^2)+(1.A(-LthroughPos+x)(-Tanh[Muldensteigung(-0.5Muldenbreite-LthroughPos+x)]+Tanh[Muldensteigung(0.5Muldenbreite-LthroughPos+x)]))/Muldenbreite^2R1secondderiv[x_]:=-0.5`A(-4(-LthroughPos+x)(-Muldensteigung+Muldensteigung)-2(-Tanh[Muldensteigung(-0.5`Muldenbreite-LthroughPos+x)]+Tanh[Muldensteigung(0.5`Muldenbreite-LthroughPos+x)])+(-)(2Tanh[Muldensteigung(-0.5`Muldenbreite-LthroughPos+x)]-2Tanh[Muldensteigung(0.5`Muldenbreite-LthroughPos+x)]))gamma[x_]:=Sqrt[1+R1firstderiv[x]^2]Plot[{R1[x]/.params,R1firstderiv[x]/.params,R1secondderiv[x]/.params},{x,0,L/.params},PlotRangeAll]Plot[gamma[x]/.params,{x,0,L/.params}]
1
2
Muldenbreite
2
Sech[Muldensteigung(-0.5`Muldenbreite-LthroughPos+x)]
2
Sech[Muldensteigung(0.5`Muldenbreite-LthroughPos+x)]
2
Muldenbreite
2
(-LthroughPos+x)
2
Muldensteigung
2
Sech[Muldensteigung(-0.5`Muldenbreite-LthroughPos+x)]
2
Muldensteigung
2
Sech[Muldensteigung(0.5`Muldenbreite-LthroughPos+x)]
In[]:=
Inits={m[x,0],c[x,0]}{mstat[x],cstat[x]};
In[]:=
Clear[x];simulation[f_,params_,inits_,s_]:=First@NDSolve[{(*Dynamics*)D[m[x,t],t](Dm/gamma[x]^2)*D[m[x,t],{x,2}]+(Dm/gamma[x]^2)*(R1firstderiv[x]/(R0+R1[x]))*D[m[x,t],x]-(Dm*R1firstderiv[x]*R1secondderiv[x]/gamma[x]^4)*D[m[x,t],x]+(f/.{mm[x,t],cc[x,t]}),D[c[x,t],t]Dc*D[c[x,t],{x,2}]+Dc*2*(R1firstderiv[x]/(R0+R1[x]))*D[c[x,t],x]-2*gamma[x]*(f/.{mm[x,t],cc[x,t]})/(R0+R1[x]),(*D[m[x,t],t]Dm*D[m[x,t],{x,2}]+(f/.{mm[x,t],cc[x,t]}),D[c[x,t],t]Dc*D[c[x,t],{x,2}]-2*(f/.{mm[x,t],cc[x,t]})/R0,*)(*Boundaryconditions*)(*(D[m[x,t],x]/.x0)==(D[m[x,t],x]/.x(L/.params)),(D[c[x,t],x]/.x0)==(D[c[x,t],x]/.x(L/.params)),m[0,t]m[(L/.params),t],c[0,t]c[(L/.params),t],*)(D[m[x,t],x]/.x0)0,(D[c[x,t],x]/.x0)0,(D[m[x,t],x]/.x(L/.params))0,(D[c[x,t],x]/.x(L/.params))0,(*Initialdata*)inits}/.params,{m[x,t],c[x,t]},{x,0,L/.params},{t,0,tend/.params},(*,Method"BDF"*)StartingStepSize{(L/.params)/256,0.0001},MaxStepSize{(L/.params)/256,10000},Method{"MethodOfLines","SpatialDiscretization"{"TensorProductGrid","MinPoints"256,"MaxPoints"256,"DifferenceOrder"1},"DifferentiateBoundaryConditions"{True,"ScaleFactor"2},"Method"{"IDA","ImplicitSolver""Newton"}}];
In[]:=
runSimulation=simulation[f[m,c],params,Inits,1];
In[]:=
Manipulate[Plot[{m[x,t]/.runSimulation/.{xxi,tti},(R0/.params)+1000*(R1[xi]/.params)},{xi,0,(L/.params)},EpilogLine[{{y,-100},{y,100}}],PlotRangeAll],{ti,0,tend/.params,1},{y,0,(L/.params)}]Manipulate[Plot[c[x,t]/.runSimulation/.{xxi,tti},{xi,0,(L/.params)},EpilogLine[{{y,-100},{y,100}}],PlotRangeAll],{ti,0,tend/.params,1},{y,0,(L/.params)}]
In[]:=
(*totalmassofsystemmustbeconservedwithnofluxboundarysorperiodicboundarysoneoptiontocheckvaldityofnumericalresult*)Ntotal=NIntegrate[2*Pi*((R0/.params)+(R1[x]/.params))*(gamma[x]/.params)*(mstat[x]/.params)+Pi*((R0/.params)+(R1[x]/.params))^2*(cstat[x]/.params),{x,0,L/.params}]NIntegrate[((m[x,t]/.runSimulation/.{xxi,t(tend/.params)})*2*Pi*(gamma[xi]/.params)*((R0/.params)+(R1[xi]/.params))+Pi*((R0/.params)+(R1[xi]/.params))^2*(c[x,t]/.runSimulation/.{xxi,t(tend/.params)})),{xi,0,L/.params},Method{Automatic,"SymbolicProcessing"0}]