I'm doing CFD simulations for blood flow in unstructured grids. My boundary condition at the outlets is called three-element Windkessel which basically calculates the pressure by solving this ODE:

$$R_{2}C\frac{dP}{dt} + P(t) = C R_{1}R_{2}\frac{dQ}{dt} + (R_{1}+R_{2})Q(t)$$

Where $P(t)$ is the pressure at the outlet, $Q(t)$ is the flux at the outlet, and $R_{1}$, $R_{2}$, and $C$ are known constants. My initial condition is: $P(0) = 0$. I discretized this ODE as:

$$P^{t+\Delta t} = \frac{1}{2\Delta t + R_{2}C}(R_{2}CP^{t-\Delta t}+2CR_{1}R_{2}(Q^{t}-Q^{t-\Delta t})+2(R_{1}+R_{2})\Delta tQ^{t})$$

My problem is that because my initial condition is set to $P(0) = 0$, it takes very long time to reach cyclic stability (my flux is periodic due to heart beat), which is very expensive computationally or even impossible for me to reach. Is there any way of having a better guess for initial condition here in order to make sure I can reach the cyclic stability faster?

  • $\begingroup$ Are you solving for P(t), for given Q(t)? $\endgroup$ Commented Dec 14, 2020 at 16:44
  • $\begingroup$ @MaximUmansky I'm solving for $P(t)$ when in each time step of CFD simulation I know $Q^{t}$ and $Q^{t-\Delta t}$ but have no idea about $Q^{t+\Delta t}$. $\endgroup$ Commented Dec 14, 2020 at 16:48
  • $\begingroup$ I am thinking of combining P and Q into W=(P-R1*Q), so the equation becomes $ d_t W = -W/(C R_2) + Q/C$. $\endgroup$ Commented Dec 14, 2020 at 17:15
  • $\begingroup$ @MaximUmansky So basically, you are suggesting $P(0) = R_{1} Q(0)$ might be a good guess? $\endgroup$ Commented Dec 14, 2020 at 17:18
  • $\begingroup$ Writing the equation in terms of W it is a little easier to see the structure of it. If the decay rate $\nu=1/(C_2 R)$ is large then after the transients are gone it is only the term $Q/C$ that matters. So the quantity $W(t)$ is determined by the time history (integral) of $Q(t)$, isn't it. Do we know how the time scales compare for $Q(t)$ vs. $\nu$? $\endgroup$ Commented Dec 14, 2020 at 17:25

1 Answer 1


Instead of solving the ODE one can use integration here.

Let's rewrite the equation in a form more useful for analysis,

$ \frac{d}{dt} W + \nu W = \frac{Q(t)}{C}, $

where $W = P - R_1 Q$, $\nu = \frac{1}{C_2 R}$

Now, use the integrating factor $u(t) = \exp(\nu t)$, and multiply the equation by $u(t)$. Then it becomes

$ \frac{d}{dt} (u W) = \frac{u Q}{C} $

Therefore we have the solution for $W$,

$ W = \frac{1}{C} \exp(- \nu t) \int^t_{0} \exp(\nu t^{\prime}) Q(t^{\prime}) dt^{\prime}, $

so by using the available time history $Q(t)$ and taking numerically the integral we determine $W(t)$, and then $P(t)$ is found as $W(t) + R_1 Q(t)$.

  • $\begingroup$ My understanding from your answer is that: choosing this as my initial condition would be a good guess for surpassing the initial decay: $$P_{0} = R_{1}Q(\frac{1}{\nu}) + \frac{1}{C} \exp(-1) \int_{0}^{\frac{1}{\nu}} \exp(\nu \tau) Q(\tau) d\tau$$ But the problem is at the beginning of the simulation, I have no idea about $\int_{0}^{\frac{1}{\nu}} \exp(\nu \tau) Q(\tau) d\tau$, which makes me think that you are suggesting to do the simulation one time until $t = \frac{1}{\nu}$ and save the history of $Q$ and then calculate $P_{0}$ introduced above as initial condition for a new simulation. $\endgroup$ Commented Dec 16, 2020 at 16:22
  • $\begingroup$ That seems a good idea if my simulation was not extraordinary expensive but sometimes continuing the simulation even until $t = \frac{1}{\nu}$ needs much resources beyond what is available to me, which makes me think that this approach is not practical for me. Also, keeping the history of $Q$ all the time from the beginning in each time step is a huge burden added to the already my expensive computations. $\endgroup$ Commented Dec 16, 2020 at 16:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.