How to Solve a DAE System with Integral Constraint and Moving Boundary in Mathematica?

Question

I am currently working on solving a DAE system in Mathematica, which includes an integral constraint and a moving boundary. I have attempted to use ParametricNDSolve, but it doesn't seem to work for my problem. Actually it origins from a research article(https://royalsocietypublishing.org/doi/10.1098/rspa.2013.0066). For simplicity, here are the details: If I ignore the constraint, simply let D=1 (Dw in my code), then I could get a solution directly from NDSolve:

\[Gamma]w = 1; \[Theta]G = 
 Pi/6; Dw = 1; eqs = {x'[s] == Cos[\[Theta][s]], 
  y'[s] == Sin[\[Theta][s]], \[Theta]''[s] == 
   nx[s]*Sin[\[Theta][s]] - ny[s]*Cos[\[Theta][s]], 
  nx'[s] == -(\[Gamma]w/(r[s]))*Sin[\[Theta][s]], 
  ny'[s] == (\[Gamma]w/(r[s]))*Cos[\[Theta][s]], 
  r'[s] == 0, \[Beta]'[s] == 0}; bc0 = {x[0] == 0, 
  y[0] == 0, \[Theta][0] == 0}; bc1 = {x[Dw] == 
   r[Dw]*Sin[\[Beta][Dw]], \[Theta]'[Dw] == 0, 
  nx[Dw] == -\[Gamma]w*Cos[\[Beta][Dw]], 
  ny[Dw] == \[Gamma]w*
    Sin[\[Beta][Dw]], \[Theta][Dw] + \[Beta][
     Dw] == \[Theta]G}; var = {x, y, \[Theta], nx, ny, 
  r, \[Beta]}; sol = NDSolve[Join[eqs, bc0, bc1], var, {s, 0, Dw}]

Next, I can find out appropriate D to fit with the constraint through changing the variable and repeating the calculation. However, such method seems inefficient.

I am wodering if there exists any direct methods to deal with such system. I tried ParametricNDSolve, however, I found it seems can't work because the position of D.

Any guidance or suggestions on how to approach this problem would be greatly appreciated!

Answer 1

To compute constraint, we can use code as it is with a small modification as follows

int[Dw_?NumericQ, tension_, wetting_] := 
  Module[{\[Gamma]w = tension, \[Theta]G = wetting}, 
   eqs = {x'[s] == Cos[\[Theta][s]], 
     y'[s] == Sin[\[Theta][s]], \[Theta]''[s] == 
      nx[s]*Sin[\[Theta][s]] - ny[s]*Cos[\[Theta][s]], 
     nx'[s] == -(\[Gamma]w/(r[s]))*Sin[\[Theta][s]], 
     ny'[s] == (\[Gamma]w/(r[s]))*Cos[\[Theta][s]], 
     r'[s] == 0, \[Beta]'[s] == 0}; 
   bc0 = {x[0] == 0, y[0] == 0, \[Theta][0] == 0}; 
   bc1 = {x[Dw] == r[Dw]*Sin[\[Beta][Dw]], \[Theta]'[Dw] == 0, 
     nx[Dw] == -\[Gamma]w*Cos[\[Beta][Dw]], 
     ny[Dw] == \[Gamma]w*
       Sin[\[Beta][Dw]], \[Theta][Dw] + \[Beta][Dw] == \[Theta]G}; 
   var = {x, y, \[Theta], nx, ny, r, \[Beta]}; 
   sol2 = NDSolve[Join[eqs, bc0, bc1], var, {s, 0, Dw}]; 
   ind = Quiet@
     NIntegrate[y[s] Cos[\[Theta][s]] /. sol2[[1]], {s, 0, Dw}, 
      Method -> "LocalAdaptive", AccuracyGoal -> 5]; (-1 + 
      r[Dw]^2 (\[Beta][Dw] - 1/2 Sin[2 \[Beta][Dw]]) + 
      2 x[Dw] y[Dw] - 2 ind ) /. sol2[[1]]];

Using ContourPlot we can solve equation for constraint as

 plot=ContourPlot[int[q1, q2, Pi/2] == 0, {q1, .8, 1.255}, {q2, .1,2.5}, 
 FrameLabel -> {"Dw", "\[Gamma]w"}] 

We can retrieve roots from plot in a form of list

roots=plot[[1]][[1]][[1]];

Finally, we can compute solution for every pare of parameters in the roots

sol[Dw_?NumericQ, tension_, wetting_] := 
  Module[{\[Gamma]w = tension, \[Theta]G = wetting}, 
   eqs = {x'[s] == Cos[\[Theta][s]], 
     y'[s] == Sin[\[Theta][s]], \[Theta]''[s] == 
      nx[s]*Sin[\[Theta][s]] - ny[s]*Cos[\[Theta][s]], 
     nx'[s] == -(\[Gamma]w/(r[s]))*Sin[\[Theta][s]], 
     ny'[s] == (\[Gamma]w/(r[s]))*Cos[\[Theta][s]], 
     r'[s] == 0, \[Beta]'[s] == 0}; 
   bc0 = {x[0] == 0, y[0] == 0, \[Theta][0] == 0}; 
   bc1 = {x[Dw] == r[Dw]*Sin[\[Beta][Dw]], \[Theta]'[Dw] == 0, 
     nx[Dw] == -\[Gamma]w*Cos[\[Beta][Dw]], 
     ny[Dw] == \[Gamma]w*
       Sin[\[Beta][Dw]], \[Theta][Dw] + \[Beta][Dw] == \[Theta]G}; 
   var = {x, y, \[Theta], nx, ny, r, \[Beta]}; 
   sol2 = NDSolve[Join[eqs, bc0, bc1], var, {s, 0, Dw}]; sol2[[1]]];


solr = (sol[#[[1]], #[[2]], Pi/2] & /@ roots) // 
    Quiet; 

Visualization function $\theta (Dw)$ dependent on $\gamma w$

listtetD = 
 Table[Quiet@
   Evaluate[{roots[[i, 2]], \[Theta][roots[[i, 1]]] /. 
      solr[[i]]}], {i, Length[roots]}];
ListLinePlot[listtetD, Frame -> True, 
 FrameLabel -> {"A\[Gamma]w/(YI)", "\[Theta](Dw)"}]

As we can see this is one curve shown in Figure 3 in the paper computed for the wetting angle $\theta_{\gamma}=\pi /2$.