Reverse flow not working when using dynamic flow rates

[schmoelder]

If dynamic flow rates are enabled in the connections, the column velocity (which used to determine the flow direction if the cross section area is also provided) is no longer taken into account.

Here are some examples:
dyn_flow_rate_examples.zip

[schmoelder]

After some debugging, @Immudzen and I found the issue:

CADET/src/libcadet/model/ModelSystemImpl-Residual.cpp

Line 499 in bc099ea

m->setFlowRates(_flowRateIn[i], _flowRateOut[i]);

In here, we call setFlowRates() which calculates the absolute flow rate of a unit operation but independent of direction.

CADET/src/libcadet/model/parts/ConvectionDispersionOperator.cpp

Line 220 in bc099ea

void ConvectionDispersionOperatorBase::setFlowRates(const active& in, const active& out, const active& colPorosity) CADET_NOEXCEPT

To update the direction, we would need to call the notifyDiscontinuousSectionTransition() of the ConvectionDispersionOperator.

CADET/src/libcadet/model/parts/ConvectionDispersionOperator.cpp

Line 186 in bc099ea

bool ConvectionDispersionOperatorBase::notifyDiscontinuousSectionTransition(double t, unsigned int secIdx)

However, we do not have access to secIdx inside the residual function. To solve this, we could - on section transition - store the current direction in the column object and include the direction directly in the calculation of the velocity.

[sleweke-bayer]

Sorry, I still don't understand. Please elaborate some more. Why is there no discontinuous section transition?

[schmoelder]

The _curVelocity is calculated in setFlowRates(). However, because it only takes into account the flow rate that enters the unit, it is at this point positive. Whenever notifyDiscontinuousSectionTransition() is called, it actually takes into account the velocity and makes _curVelocity negative (if required). However, for dynamic flow rates, the velocity changes not only at section times, but is recalculated for every time step and is therefore mostly never set to its negative value.

[schmoelder]

I already have an idea how to fix it; I will create a PR in case this is going somewhere. ^^

[sleweke-bayer]

Ah, I see. The _curVelocity is correct in the call to residual() after notifyDiscontinuousSectionTransition(). But it is overwritten in subsequent calls to residual(). This happens pretty much instantly (in terms of time) if dynamic flow rates are configured.

I haven't given it much thought, but what about preserving the sign in setFlowRates()?

void ConvectionDispersionOperatorBase::setFlowRates(const active& in, const active& out, const active& colPorosity) CADET_NOEXCEPT
{
	// If we have do not have cross section area, interstitial velocity is specified explicitly in the model
	if (_crossSection <= 0.0)
		return;

	if (_curVelocity < 0.0)
		_curVelocity = -in / (_crossSection * colPorosity);
	else
		_curVelocity = in / (_crossSection * colPorosity);
}

[schmoelder]

I guess that should also work, I will give it a try! Otherwise I would have suggested to simply add dir as an attribute

void ConvectionDispersionOperatorBase::setFlowRates(const active& in, const active& out, const active& colPorosity) CADET_NOEXCEPT
{
	// If we have cross section area, interstitial velocity is given by network flow rates
	if (_crossSection > 0.0)
		_curVelocity = _dir * in / (_crossSection * colPorosity);
}

and modify the notifyDiscontinuousSectionTransition() method correspondingly, s.t. it flips sign if _dir_old != _dir. Semantically it would also separate velocity (which, except for the parameter that is read) would then always mean the interstitial velocity, and direction, which would always mean its sign.

[sleweke-bayer]

Both solutions are fine with me 😄

[schmoelder]

bool ConvectionDispersionOperatorBase::notifyDiscontinuousSectionTransition(double t, unsigned int secIdx)
{
	// If we don't have cross section area, velocity is given by parameter
	double _dir_old = _dir;
	if (_crossSection <= 0.0)
		_dir = getSectionDependentScalar(_velocity, secIdx);
	else if (!_velocity.empty())
	{
		// We have both cross section area and interstitial flow rate
		// _curVelocity has already been set to the network flow rate in setFlowRates()
		// the direction of the flow (i.e., sign of _curVelocity) is given by _velocity
		_dir = static_cast<double>(getSectionDependentScalar(_velocity, secIdx));
		if (_dir_old * static_cast<double>(_dir) < 0.0)
			_curVelocity *= -1.0;
	}

	return (_dir_old * static_cast<double>(_dir) < 0.0);
}