Hello!
I have a prolem when converting formula into fenics code.
I am implementing the CBS algrithm in book THE FINITE ELEMENT METHOD FOR FLUID DYNAMICS wrote by Zienkiewicz, here are formulas.
\Delta U_i = U_i^{n+1}-U^{n}_i=\Delta U^{**}_i-\Delta U^{*}_i
\begin{array}{l}{\int_{\Omega} N_{u}^{a} \Delta U_{i}^{*} \mathrm{d} \Omega} \\ {\qquad \begin{aligned}=&-\Delta t\left[\int_{\Omega} N_{u}^{a} \frac{\partial}{\partial x_{j}}\left(u_{j} U_{i}\right) \mathrm{d} \Omega+\int_{\Omega} \frac{\partial N_{u}^{a}}{\partial x_{j}} \tau_{i j} \mathrm{d} \Omega-\int_{\Omega} N_{u}^{a}\left(\rho g_{i}\right) \mathrm{d} \Omega\right]^{n} \\ &+\frac{\Delta t^{2}}{2}\left[\int_{\Omega} \frac{\partial}{\partial x_{k}}\left(u_{k} N_{u}^{a}\right)\left(-\frac{\partial}{\partial x_{j}}\left(u_{j} U_{i}\right)+\rho g_{i}\right) \mathrm{d} \Omega\right]^{n} \\ &+\Delta t\left[\int_{\Gamma} N_{u}^{a} \tau_{i j} n_{j} \mathrm{d} \Gamma\right]^{n} \end{aligned}}\end{array}
\begin{aligned} \int_{\Omega} N_{p}^{a} \Delta \rho \mathrm{d} \Omega=&-\Delta t \int_{\Omega} N_{p}^{a} \frac{\partial}{\partial x_{i}}\left(U_{i}^{n}+\theta_{1} \Delta U_{i}^{*}-\theta_{1} \Delta t \frac{\partial p^{n+\theta_{2}}}{\partial x_{i}}\right) \mathrm{d} \Omega \\=& \Delta t \int_{\Omega} \frac{\partial N_{p}^{a}}{\partial x_{i}}\left[U_{i}^{n}+\theta_{1}\left(\Delta U_{i}^{*}-\Delta t \frac{\partial p^{n+\theta_{2}}}{\partial x_{i}}\right)\right] \mathrm{d} \Omega \\ &-\Delta t \int_{\Gamma} N_{p}^{a}\left[U_{i}^{n}+\theta_{1}\left(\Delta U_{i}^{*}-\Delta t \frac{\partial p^{n+\theta_{2}}}{\partial x_{i}}\right)\right] n_{i} \mathrm{d} \Gamma \end{aligned}
\begin{aligned} \int_{\Omega} N_{u}^{a} \Delta U_{i}^{\star \star} \mathrm{d} \Omega=& \int_{\Omega} N_{u}^{a} \Delta U_{i} \mathrm{d} \Omega-\int_{\Omega} N_{u}^{a} \Delta U_{i}^{\star} \mathrm{d} \Omega \\=&-\Delta t \int_{\Omega} N_{u}^{a}\left(\frac{\partial p^{n}}{\partial x_{i}}+\theta_{2} \frac{\partial \Delta p}{\partial x_{i}}\right) \mathrm{d} \Omega \\ &-\frac{\Delta t^{2}}{2} \int_{\Omega} \frac{\partial}{\partial x_{j}}\left(u_{j} N_{u}^{a}\right) \frac{\partial p^{n}}{\partial x_{i}} \mathrm{d} \Omega \end{aligned}
here is my fenics code(after omiting boundary terms):
def grad_p(p_n, p, theta):
return (1-theta)*grad(p_n)+theta*grad(p)
F1 = - dot(u, v)*dx - dt*dot(dot(u_n,grad(u_n)),v)*dx - dt*mu*inner(grad(u_n),grad(v))*dx + 0.5*dt*dt*dot(dot(u_n,grad(u_n)),dot(u_n,grad(v)))*dx
F2 = dot(u_n+theta1*(u_ - dt*grad_p(p_n, p, theta2)),grad(q))*dx
F3 = dot(-u + u_n + u_, v)*dx - dt*dot(grad_p(p_n,p_,theta2),v)*dx - 0.5*dt*dt*dot(grad(p_n),dot(u_n,grad(v)))*dx
I can not get a right answer.