Phase-contrast MRI (PC-MRI) data have been used in the literature to perform MRI-constrained computational fluid dynamics (CFD) simulations using the finite volume method (FVM) with SIMPLER algorithm^[1]-[4]. However, the finite element method (FEM) is more flexible on dealing with complex geometries and boundary conditions, also allowing higher order approximations than FVM^[5]. The aim of this work is to show that it is feasible to construct a more robust, hybrid MRI-CFD solver, based on FEM, that enforces PC-MRI measurements directly into the CFD routine. The hybrid solver can be used to regularize PC-MRI data, either reducing noise, or obtaining simulations that are closer to PC-MRI measurements than CFD alone. The feasibility of this approach is demonstrated using a modified 2D implementation of the Navier-Stokes-continuity equations, using the Garlekin FEM. In this demonstration, the CFD solution was constrained over a restricted 2D domain by two velocity components, taken from a 4D PC-MRI dataset.

4DFT FGRE PC-MRI data were acquired from a flow phantom (Fig.1). PC imaging was performed on a 3T GE Discovery MR750 system (50mT/m, 200T/m/s), using a 32-channel head coil (resolution = 1.0×0.5×0.5mm³; FOV = 7.5×16×12cm³; Venc = 50cm/s; TR = 11.4ms; flip angle = 8.5^o; temporal resolution = 91.2ms; scan time = 40 minutes; 9 NEX; pulse cycle 60bpm). The phantom’s blood-mimicking fluid was modelled as a Newtonian and incompressible fluid. Hence, the steady incompressible 2D Navier-Stokes-continuity system of equations^[6], $$\mathbf{u}\cdot\nabla\mathbf{u}=\nabla\cdot\boldsymbol{\sigma}\,\,\,\mbox{and}\,\,\,\nabla\cdot\mathbf{u}=0,$$ were discretized using FEM^[5], where Re is the Reynolds number^[6], $$$\mathbf{u}=u_x\mathbf{i}+u_y\mathbf{j}\in\mathbb{R}^2$$$ is the velocity field, $$$p$$$ is pressure, and $$$\boldsymbol{\sigma}=-p\mathbf{I}+Re^{-1}[\nabla\mathbf{u}+\nabla\mathbf{u}^T]$$$ is the Newtonian stress tensor^[6]. The discretization uses the residues functions given by the governing equations' weak form^[5], $$R_m(\mathbf{u},p)=\int_\Omega(\mathbf{u}\cdot\nabla\mathbf{u})\cdot\boldsymbol{\psi}d\Omega+\int_\Omega\boldsymbol{\sigma}:\nabla\boldsymbol{\psi}d\Omega-\int_\Gamma(\mathbf{n}\cdot\boldsymbol{\sigma})\cdot\boldsymbol{\psi}d\Gamma$$ and $$R_c(\mathbf{u})=\int_\Omega(\nabla\cdot\mathbf{u}){\phi}d\Omega,$$ where $$$\phi\in\mathbb{R},\,\,\boldsymbol{\psi}\in\mathbb{R}^2$$$ are weight functions. The discretized equations were written as a linear system $$$\mathbf{J}\mathbf{c}=-\mathbf{R}$$$^[5], where $$$\mathbf{J}$$$ is the Jacobian matrix of the residues, $$$\mathbf{R}$$$ is the residues vector, and $$$\mathbf{c}$$$ is the solution vector, containing velocity and pressure. To incorporate PC-MRI measurements into the CFD solver, the minimization $$\min_{\mathbf{c}}{\frac{1}{2}\|\mathbf{J}\mathbf{c}+\mathbf{R}\|_2^2+\frac{\lambda}{2}\|\mathbf{S}\mathbf{c}-\mathbf{u}_{\mathrm{mri}}\|_2^2}$$ is performed using Newton's method, where $$$\mathbf{S}$$$ is a matrix used to compare velocities on solution vector $$$\mathbf{c}$$$ with the measured PC-MRI signal $$$\mathbf{u}_{\mathrm{mri}}$$$, since these vectors may have different sizes. The lumen was manually outlined to define the computational mesh. Simulation grid was designed with 1.0×0.5mm² element resolution using Q₂P_-1 elements^[5]. The Reynolds number for the problem was calculated as Re = 110. The PC-MRI velocity profile was set at the inlet of the flow domain together with no-slip boundary condition. MRI-constrained CFD solutions were obtained using two MRI-measured velocity components for different λ. A denoising experiment was performed, by adding 7.5cm/s standard deviation Gaussian noise to the measured PC-MRI data. Finally, the results were quantitatively and qualitatively compared with the "ground truth" 2D PC-MRI data*.

Figure 2 shows the u_x velocity component for: PC-MRI (Fig.2(a)), pure CFD solution, λ = 0 (Fig.2(b)), and PC-constrained CFD solution, λ = 10^-2 (Fig.2(c)). The constrained solution (Fig.2(c)) is closer to the measured PC-MRI data (Fig.2(a)) than pure CFD (Fig.2(b)). Moreover, this behavior is quantitatively confirmed using the signal-to-error ratio (SER), considering the measured PC-MRI as "ground truth" signal. While the pure CFD solution provided 4.53dB SER, the MRI-constrained CFD solution provided 8.17dB SER. This shows that constrained CFD solution can produce more realistic velocity fields. In Fig.2(b), there is a noticeable misbehavior of the solution in the center of the domain, near the bifurcation. For this problem’s Reynolds number, the Navier-Stokes equation’s diffusive term $$$(\Delta^2\mathbf{u})$$$ is much smaller than the convective term $$$(\mathbf{u}\cdot\nabla\mathbf{u})$$$ causing amplification of solution errors. Then, this solver implementation requires numerical stabilization treatment^[7]. Figure 3 shows the velocity fields from the denoising experiment for: PC-MRI (Fig.3(a)), noise corrupted PC-MRI, (Fig.3(b)) and CFD solution constrained by noisy PC-MRI, λ = 10^-3 (Fig.3(c)). The constrained solution (Fig.3(c)) corrected the noisy PC flow (Fig.3(b)), leading to a solution that is very similar to original PC-MRI. This was quantitatively confirmed using the signal-to-noise ratio (SNR), considering PC-MRI as "ground truth" signal. The noisy PC-MRI has 5.39dB SNR, whereas the CFD solution constrained by noisy PC-MRI provided 6.86dB SNR. This 1.47dB SNR increase shows that MRI-constrained CFD can also be used for denoising PC-MRI velocity fields.

This work shows that the proposed hybrid CFD-MRI finite-element solver is viable and can be used to obtain numerical solutions that are similar to the measured flow data, while satisfying the fluid physics equations; and also for reducing noise. In future works, a streamline-upwind stabilization term^[7], for high Reynolds number, will be implemented on a 3D steady version of this FEM based solver.