GenBC of parallel RCR network

[Hadi Wiputra]

Hi All,

I recently try to model 2 parallel RCR circuits:

parallel RCR network.PNG (9.55 KiB) Viewed 417 times

I ran this in GenBC with the following differential equations:
c index: P=1 Q1=2 P1=3 Q2=4 P2=5
f(3)= (1D0/C(1)) * (x(2) - x(3)/Rd(1))
f(5)= (1D0/C(2)) * (x(4) - x(5)/Rd(2))
f(2)=(f(1)-f(3))/Rp(1)
f(4)=(f(1)-f(5))/Rp(2)
f(1)=(1D0/(C(1)+C(2)))*(Q(1)-x(3)/Rd(1))
f(1)=f(1)+(1D0/(C(1)+C(2)))*(-x(5)/Rd(2))
f(1)=f(1)+(1D0/(C(1)+C(2)))*(C(1)*Rp(1)*f(2))
f(1)=f(1)+(1D0/(C(1)+C(2)))*(C(2)*Rp(2)*f(4))

In the equations, many of the differentials (f-es) is dependent on other f-es. I noticed that GenBC requires me to define the differentials (f-es) based on only x-es (since f-es passed to this function are just empty containers).

Is there a set way on how to generate these differential only as a function of x-es?

Thank you.

Regards,
Hadi

[Weiguang Yang]

You can try to use the "offset" function. The offset array (size= number of unknowns) will add a constant value to the corresponding unknowns before they are returned to the 3D solver. The default values for offset are 0. For an RCR, offset can be used to make the RCR ODE in GenBC simple. Here is an example. https://simvascular.github.io/docsGenBC.html


For your parallel RCR, you can set
f(1)=dP1/dt=Q1/C1
f(2)=dP2/dt=Q2/C2

Q1+Q2=Q
Q1*Rp1+P1=Q2*Rp2+P2

offset(1)=Q1*Rp1
offset(2)=Q2*Rp2

[Hadi Wiputra]

Hi Weiguang,

Thank you for your help. However, I am not sure if I understand it well enough.

If I understood it properly, we can set four variables, x(1:4), with index 1=P1, 2=P2, 3=Q1, 4=Q2. Q is given by the program.
These ones are straight forward:
f(1)= x(3)/C1
f(2)=x(4)/C2
x(3)=Q-x(4)
However, x(4)=(x(3)*Rp1+x(1)-x(2))/Rp2 -> rearranged from Q1*Rp1+P1=Q2*Rp2+P2
It means the x(3) is dependent on x(4) and x(4) is dependent on x(3), which one should be updated first? (have this same issue in the prev equation)

Also, there are 2 offsets, but there is only one Neumann boundary (P), how do I sum the two offsets and ensure that offset + P1 = offset + P2? Should I only use one offset and ignore the other?

I appreciate your help so far, I am sorry I am a bit slow.

Regards,
Hadi

[Hadi Wiputra]

Hi Weiguang,

I've also realised that the dP1/dt is not a function Q1, it is a function of Q11, where Q1=Q11+Q12, since Q1 is the first split from the parent branch and Q11 is the split from the child branch

Regards,
Hadi

[Weiguang Yang]

Sorry, the ODE needs to be corrected.

I think you can solve Q1 and Q2 first:
Q1+Q2=Q
X(1)+Q1*Rp1-(X(2)+Q2*Rp2)=0

where X(1) and X(2) are pressures over the capacitors C1 and C2 respectively.

then plug X(1),X(2), Q1 and Q2 into f(1) and f(2):

f(1)=dP1/dt=(Q1-P1/Rd1)/C1=(Q1-X(1)/Rd1)/C1
f(2)=dP2/dt=(Q2-P2/Rd2)/C2=(Q2-X(2)/Rd2)/C2

Since there are two capacitors, you just need to use 2 unknowns.

Let
offset(1)=Q1*Rp1
offset(2)=Q2*Rp2
Since you have 1 Neumann pressure, you can choose either X(1) or X(2) as the corresponding unknown for the Neumann pressure. The pressure applied to the Neumann face will be X(1)+offset(1) or X(2)+offset(2) (This is done in GenBC.f)

I hope this makes sense.
Thanks

[Hadi Wiputra]

Hi Weiguang,

I get it now, I will try it out, thank you.

Regards,
Hadi

[Hadi Wiputra]

Hi Weiguang,

I did the GenBC simulations and it runs, then i used the GUI to extract the flowrate at each inlet/outlet surfaces. I found out that the sum of the flow from the inlets and outlets didn't add up to zero. The simulation I am running is a rigid wall, so it should be balanced. Somehow it is not. I wonder if it is the issue with the extraction of data or the run itself.

Edit: the case files are attached in google drive: https://drive.google.com/drive/folders/ ... sp=sharing

Attached is the output from the GUI. I tried the case with built in RCR boundary (not GenBC) and the in and outflow rate is balanced.

Update: I looked into the Q_General file generated, summed it up (blue line) and plot it against the inflow dirichlet boundary (in red). I think the issue is with the GenBC somehow.

all_results-flows.txt

(1.4 KiB) Downloaded 8 times

Update, Update: It is solved now, the mystery is in the timestep. But not in the RK timestep, it is in the simulation time step. Since the flow rate stepping in RK is linear, changing the RK timestep does not help at all. I have to increase the timestep from 200/cycle to about 500/cycle to ensure the flow balance is maintained. I wonder if this can be improved if the value of Q passed to RK is based on non linear interpolation.. Anyway, mystery is solved, thank you


Regards,
Hadi

[Weiguang Yang]

I feel the difference between the inflow and total outflow is more likely associated with the 3D solution strategies instead of the boundary condition. The mesh size, number of nonlinear iteration and step size have impacts on the continuity constraint. I don't think the interpolation method between t and t+dt in Genbc is critical here for flow conservation as long as the time step you use is reasonable. You can also use the outlet flow data as input to numerically integrate the lumped parameter model in matlab and check if they resulting pressures are consistent with the pressure applied to the outlet pressure by GenBC.