3.1. Example 1: flow over an oscillating cylinder

The performance of the proposed method is first validated using a simulated 2-D flow over an oscillating cylinder at a diameter-based Reynolds number of $Re_D$ = 185, where $D$ denotes the diameter of the cylinder. The computational domain spans $-15D\leqslant x \leqslant 45D$ and $-20D\leqslant y \leqslant 20D$ , and the cylinder is initialised at the origin (0, 0) and moves in the cross-flow (y-axis) direction according to $y_c=0.2sin(2\pi f_c t)$ . The normalised parameter $f_c$ is set to 0.156 according to Wang et al. (Reference Wang, Vanella and Balaras2019). The ground truth velocity $\mathbf {U}$ and pressure fields $p$ are obtained by solving the non-dimensional incompressible N–S equations:

(3.1) \begin{align} \frac {\partial \mathbf {U}}{\partial t}+(\mathbf {U}\cdot \nabla )\mathbf {U} &=-\nabla p+\frac {1}{Re_D}\nabla ^2\mathbf {U}, \nonumber\\ \nabla \cdot \mathbf {U} &= 0. \end{align}

All variables are normalised by the fluid density $\rho$ , the free streamwise velocity $U_{\infty }$ and the diameter $D$ . This simulation is conducted using OpenFoam with a moving grid strategy over a body-fitted and wall-refined mesh. The non-dimensional time step in the simulation is set to $2.5\times 10^{-3}$ , with a maximum Courant–Friedrichs–Lewy number of 0.7. A no-slip condition is imposed at the cylinder surface. Note that the data are stored every twenty steps, which means that the stored time interval is $5.0\times 10^{-2}$ . After the CFD data are obtained, the scattered velocity data for training are collected by interpolating the flow fields on a uniform grid to mimic the typical PIV data. The inspected region is limited to the subdomain of $-1D\leqslant x \leqslant 3D$ and $-1D\leqslant y \leqslant 1D$ . The vector space is denoted by $1/R$ , where $R$ is the number of points along the diameter. A small $R$ indicates a low spatial resolution. There are 721 boundary points distributed uniformly on the surface of the cylinder, and the time sequence is the same as that of the velocity data. In addition, evaluation of the influence of measurement noise on the performance of the proposed method is vital. Therefore, Gaussian random noise is added to the training data (including the boundary data). The noise level $\sigma_n$ is defined as the ratio of the mean noise magnitude to the velocity magnitude. The signal-to-noise ratio (SNR) can be simply estimated as $\textrm{SNR}=10{\log}_{10}{1/\sigma _n^{2}}$ . We randomly sampled $1\times 10^6$ data points and $0.5\times 10^6$ boundary points from four hundred velocity fields without pressure information as training data, $1\times 10^6$ points to estimate the residuals of the N–S equations and $0.5\times 10^6$ points to estimate the boundary conditions.

The PINN-MB network is built for pressure reconstruction in scenarios with a moving boundary. The boundary network contains nine hidden layers and 32 neurons per hidden layer, and the flow network contains 13 hidden layers and 128 neurons per hidden layer, yielding a total of 207 707 trainable parameters. This means that we used 207 707 parameters to capture the spatiotemporal evolution of the flow over the oscillating cylinder, whereas the numerical simulation used approximately $1.4\times 10^7$ grid points (400 velocity fields, each field containing 35 940 CFD points) in the subdomain to compute the flow fields. The network was trained for 4000 epochs with 30 iterations per epoch, and the batch size was set to 10 000. The initial position $\mathbf {X}_b^0$ was set to a standard circle with $D = 1$ without measurement noise. The boundary noise at the initial position can reasonably be ignored because the experimental model is accurately designed in advance. No pressure information is provided for training. An example snapshot of the velocity and pressure is shown in figure 2. The training data with $R = 20$ and $\sigma _n = 0\, \%$ , the data learned over a high-resolution mesh with $R=100$ , and the absolute error between the exact and learned data are presented from top to bottom. From left to right, the contour maps represent $u$ , $v$ and $p$ , respectively. At this instant, the cylinder moves forward along the y-axis, as indicated by the yellow arrows. The pressure is greater in the stagnation region and lower in the vortex region. Notably, the mean pressure at each time instant is subtracted to enable comparative assessments. A visual inspection of the results suggests that the predicted high-resolution data are approximately the same as the exacted data, and the maximum absolute error of the pressure coefficient is approximately $4\, \%$ of the maximum pressure. In figure 3, noise with $\sigma _n=8\, \%$ was added to both the velocity data and boundary data. The velocity, pressure and boundary of the cylinder can still be recovered from the noisy data, although the absolute errors are larger than those in the case without noise. As shown in figures 2 and 3, a large error ( $\sim$ 0.05) is always observed in the region near the boundary for high velocity gradients. This visual inspection shows that the proposed method can accurately determine the pressure fields in scenarios with a moving boundary.

Figure 2.

A representative snapshot of the velocity, pressure and corresponding absolute error. From left to right, the contour maps represent $u$ , $v$ and pressure $p$ . The dashed black lines represent negative contour levels from $-$ 1.5 to $-$ 0.1 with an interval of 0.2, and the solid black lines represent positive contour levels from 0.1 to 1.5 with an interval of 0.2. Note that the spatially averaged pressure is subtracted from both the exacted data and the predicted data for comparative assessment. For the training data, the vector space $1/R$ is set to $1/20$ , and the noise level is set to $0\, \%$ . The learned data were predicted over a mesh with an interval of $1/R=1/100$ . Therefore, the original CFD data were interpolated on the predicted mesh to compute the absolute error. The yellow arrows indicate the motion of the cylinder.

Figure 3.

A representative snapshot of the velocity, pressure and corresponding absolute error. From left to right, the contour maps represent $u$ , $v$ and pressure $p$ . The dashed black lines represent negative contour levels from $-$ 1.5 to $-$ 0.1 with an interval of 0.2, and the solid black lines represent positive contour levels from 0.1 to 1.5 with an interval of 0.2. Note that the spatially averaged pressure is subtracted from both the exacted data and the predicted data for comparative assessment. For the training data, the vector space $R$ is set to $1/20$ , and the noise level is set to $8\, \%$ (the corresponding SNR is 21.9). The learned data were predicted over a mesh with an interval of $R=1/100$ . Therefore, the original CFD data were interpolated over the predicted mesh to compute the absolute error. The yellow arrows indicate the motion of the cylinder.

To provide a quantitative comparison, we present the relative $\mathcal {l}_{2}$ norm error of the pressure as a function of $R$ and $\sigma _n$ in figure 4. The relative $\mathcal {l}_{2}$ norm error is defined as follows:

(3.2) \begin{equation} \delta _p=\frac {\Vert p_{ref}-p_{pred}\Vert _{2}}{\Vert p_{ref}\Vert _{2}}. \end{equation}

The subscripts $ref$ and $pred$ denote the reference value and the predicted value, respectively. The relative $\mathcal {l}_{2}$ norm error of the velocity components is also estimated using the same formula. All errors are averaged over space and time. Before presenting the results, we note that PINN-based pressure reconstruction has been compared with traditional approaches within pure flow fields (Fan et al. Reference Fan, Xu, Wang and Wang2023). The previous work has shown that the PINN-based method is superior to the other two methods because the former can optimise the velocity fields. In the present work, a dual-network model is used to determine the pressure from velocity with a moving boundary. This method can address not only rigid bodies but also flexible boundaries, as shown in the next three examples. Returning to figure 4, it is evident that the pressure error ( $\sim 4\, \%$ ) is larger than the velocity error, which is consistent with the results of Wang et al. (Reference Wang, Liu and Wang2022a ) and Fan et al. (Reference Fan, Xu, Wang and Wang2023). For the velocity, the relative error in the $v$ component is slightly larger than that in the $u$ component for the difference in magnitude. These figures show that the impact of the sparsity of the data on the pressure reconstruction performance is much smaller than that of the measurement noise level. The pressure error at $R = 10$ , which corresponds to only ten points along the diameter, is nearly identical to that at $R = 160$ . However, the pressure error increases from $1\, \%$ to $4\, \%$ as the noise level $\sigma _n$ increases from $1\, \%$ (the corresponding SNR is approximately 40) to $16\, \%$ (the corresponding SNR is approximately 16) in this test case. These results demonstrate that reducing measurement noise and improving data accuracy are critical for pressure reconstruction from measured velocity data. Based on the quantitative analysis, we conclude that the proposed PINN-MB method is feasible and accurate for pressure reconstruction with a moving boundary.

Figure 4.

Relative $\mathcal {l}_{2}$ norm error as a function of noise level $\sigma _n$ and vector space $1/R$ : (a) $u$ component, (b) $v$ component and (c) pressure. In each case, the same PINN-MB model is used. The boundary network contains nine hidden layers and 32 neurons per hidden layer, and the flow network contains 13 hidden layers and 128 neurons per hidden layer. The red curves represent the mean values along the noise level or vector space.

Figure 5(a) presents the pressure coefficient around the cylinder at the time when the cylinder moves to the maximum vertical position of $y_c=0.2$ . The parameter $\theta$ denotes the angle of the surface points, as shown in figure 5(a). Note that the curves are obtained by averaging the data over all the instants when the cylinder reaches the maximum vertical position. Two prediction results with $R=20, \sigma _n=1\, \%$ and $R=160, \sigma _n=8\, \%$ are displayed in this figure. Both curves show good agreement with the CFD data, but the pressure coefficient predicted with the higher noise level has a larger deviation from the CFD data in the region of $80^\circ \leqslant \theta \leqslant 120^\circ$ . Moreover, we present the power spectral density (PSD) of the pressure at $(2,0)$ in figure 5(b). The parameter $f$ denotes the non-dimensional frequency. The distributions of the spectra obtained using the proposed neural network are consistent with the reference CFD data in the region of $0.2\leqslant f\leqslant 3$ . The deviation at a low frequency of $f=0.1$ may be attributed to the larger error in the predicted pressure at the beginning and end of the corresponding time, as shown in the top plane of figure 6. Additionally, we present the cylinder motion, i.e. the y-axis centroid of the cylinder, in this figure. No evident difference among the three curves are observed. We next analyse the PSD with high frequencies $f\gt 3$ , which shows that the predicted PSD is much lower than that of the reference curve in the high-frequency region. There are two possible reasons for this result. First, the interpolation from the CFD results to a regular mesh may cause errors at high wavenumbers. Second, the network may prefer to learn low-frequency modes rather than high-frequency modes, a phenomenon known as ‘spectral bias’ (Rahaman et al. Reference Rahaman, Baratin, Arpit, Draxler, Lin, Hamprecht, Bengio, Courville, Chaudhuri and Salakhutdinov2019). There are several methods for overcoming this well-known issue. Eldan & Shamir (Reference Eldan and Shamir2016) proved that a simple function expressed by a small 3-layer feedforward neural network cannot be approximated by any 2-layer network unless its width is exponential in dimension. This implies that a deeper network is needed for learning at a high frequency. Tancik et al. (Reference Tancik, Srinivasan, Mildenhall, Fridovich-Keil, Raghavan, Singhal, Ramamoorthi, Barron and Ng2020) showed that mapping the input variables via Fourier features enables a multi-layer perceptron to learn high-frequency functions in low-dimensional problem domains. This Fourier feature mapping significantly improves the high-frequency performance in the neural radiance field for view synthesis (Mildenhall et al. Reference Mildenhall, Srinivasan, Tancik, Barron, Ramamoorthi, Ng, Vedaldi, Bischof, Brox and Frahm2020). Moreover, for the super-resolution reconstruction of turbulence, generative adversarial networks have the superior capabilities for representing small-scale structures compared with convolution neural networks (Kim et al. Reference Kim, Kim, Won and Lee2021).

Figure 5.

(a) Pressure coefficient around the cylinder at the instant when the vertical position of the cylinder reaches a maximum of $y_c=0.2$ . The parameter $\theta$ denotes the angle of the surface points. (b) PSD of the reference pressure and the predicted pressure at $(2,0)$ , where $f$ denotes the nondimensional frequency.

Figure 6.

Time series of the pressure coefficients $C_p$ at point (2,0) (the top panel) and the cylinder location $y_c$ (the bottom panel). $y_c$ denotes the y-axis centroid of the cylinder.

3.2. Example 2: fish swimming

A more complex case involves the reconstruction of the pressure around a swimming fish. Fish locomotion has attracted considerable attention from engineers and scientists because understanding the hydrodynamics of the fish is important for bionic underwater vehicle design. The fish uses the pressure signals detected along the lateral line to sense the surrounding environment (Calicchia et al. Reference Calicchia, Mittal, Seo and Ni2023), and generates propulsive forces by creating jets and vortices in the wake flow (Triantafyllou et al. Reference Triantafyllou, Triantafyllou and Yue2000; Ting & Yang Reference Ting and Yang2009; Tu et al. Reference Tu, Wang, Wang, Gao and Wei2022). A two-dimensional direct numerical simulation (DNS) of a swimming fish is carried out using an in-house CFD solver (Liu et al. Reference Liu, Wang and Ye2018, Reference Liu, Wang, Wang and He2023) to validate the accuracy of the proposed method.

A two-dimensional NACA0012 model is used to mimic the fish body. The fish axis lies on the x-axis, and the head and tail are located at (0,0) and (1,0), respectively. The lateral displacement of the fish axis is defined as the carangiform swimming style, as measured by Videler & Hess (Reference Videler and Hess1984). This motion is prescribed as the product of an envelope $A_m(x)$ and a backward travelling wave as follows:

(3.3) \begin{equation} y/L = A_m(x)\sin\!\left(\frac {2\pi }{\lambda }(x-ct)\right). \end{equation}

Here, $A_m(x)=b_0+b_1(x/L)+b_2(x/L)^2$ , $b_0=0.02$ , $b_1=-0.08$ and $b_2=0.16$ . The wavelength $\lambda$ is equal to the fish length $L$ , and the wave velocity $c$ is set to the product of the tail beat frequency $f$ and the wavelength $\lambda$ . In this case, the parameters are set to the following values: $L$ = 0.1 m, $f$ = 2 Hz and the free stream velocity $U$ = 0.1 m s^-1. Therefore, the Reynolds number, which is based on the body length and free stream velocity, $Re=LU/\nu$ , is $1\times 10^4$ , and the Reynolds number based on the body length and tail beat frequency, $Re=L^2f/\nu$ , is $2\times 10^4$ . The non-dimensional time step normalised by the characteristic time $L/U$ in the simulation is $\textrm{d}t=2\times 10^{-3}$ , whereas the time interval of the stored data is $\Delta t=2{\textrm d}t=4\times 10^{-3}$ for every two time steps to store the flow fields. Two hundred flow fields interpolated into a uniform grid in the domains $-0.3\leqslant x/L \leqslant 1.3$ and $-0.3\leqslant y/L \leqslant 0.3$ are postprocessed to train the network, and the vector space is set to $L/400$ . Boundary data with 682 points were extracted exactly from the numerical simulation, were non-uniformly distributed on the surface of the fish, and the time sequence was the same as that of the velocity data. To test the effect of the temporal resolution on the performance of the proposed method, data for five cases with time steps $\Delta t$ = $2{\textrm d}t$ , $4{\textrm d}t$ , $6{\textrm d}t$ , $8{\textrm d}t$ and $10{\textrm d}t$ are generated to train the network. The numbers of corresponding velocity fields are 200, 100, 67, 50 and 40. For each case, we randomly sample $2\times 10^6$ data points and $1\times 10^6$ boundary points without pressure information as training data, $2\times 10^6$ residual points to estimate the loss of N–S equations and $1\times 10^6$ points are used for the boundary conditions. In this example, we do not consider the measurement noise.

The configuration of the PINN-MB model is set as follows: the boundary network contains nine hidden layers and 48 neurons per hidden layer, and the flow network contains 13 hidden layers and 160 neurons per hidden layer, yielding a total of 329 371 trainable parameters. The number of trainable parameters is much less than the CFD points in this subdomain (two hundred velocity fields with 76 982 CFD points per field). This network is trained for 18 000 epochs with 30 iterations per epoch and the batch size is set to 10 000. Note that all of the velocities, pressures and boundaries are still predicted in the vector space of $L/400$ with a time separation of $2\textrm{d}t$ . We first investigate whether PINN-MB can capture fish motion under these complex conditions. Figure 7 shows a comparison of the predicted fish bodies with those extracted from the CFD results. The red and blue curves represent the predicted results when $\Delta t$ is equal to $2{\textrm d}t$ and $10{\textrm d}t$ , respectively. It is clear that the predicted curves are approximately the same as the original grey curves at different points in time, as the loss of boundary data for all the cases converges to less than $5.0\times 10^{-6}$ . We also plot the error of the predicted body surface as a function of the time steps of the training data in figure 7(b). The boundary error normalised by the body length is less than $0.1\, \%$ over the whole test region. This result indicates that PINN-MB can accurately predict the boundary for the coupling between the boundary and flow networks through the boundary conditions. Notably, the error of the predicted boundary is slightly larger at the beginning and end of the time ( $\sim 0.2\, \%$ ) because fewer constraints are imposed at these instants.

Figure 7.

(a) Exact (grey thick curves) and predicted fish surfaces with $\Delta t = 2{\textrm d}t$ (red curves) and $10{\textrm d}t$ (blue curves) at different instants. From top to bottom, the non-dimensional time separation is 0.12. (b) Predicted error of the fish body as a function of the time steps of the training data $\Delta t$ , the error bars represent the standard deviation. The error is normalised by the body length $L$ . Note that, all of the data, including velocities and boundaries, are predicted at the time interval of $\Delta t = 2{\textrm d}t$ .

Figure 8.

Representative snapshots of pressure for DNS, PINN-MB with $\Delta t = 10 {\textrm d}t$ and the absolute error. From top to bottom, the non-dimensional times are 0.23, 0.35 and 0.47, respectively. In this case, the vector intervals of both the training data and the predicted data are set to the same values of $L/400$ . No noise is added to the training data. The yellow arrows indicate the motion of the fish’s body. The mean value is subtracted from both the exacted data and the predicted data for comparative assessment.

The instantaneous pressure around the fish at different time instants is shown in figure 8. The CFD results, PINN-MB results with $\Delta t = 10 {\textrm d}t$ and absolute error are displayed from left to right. The yellow arrows indicate the motion of the fish’s body. The mean value is subtracted from the exacted and predicted data for comparative assessment. We do not show the predicted pressure with high temporal resolution $\Delta t = 2 {\textrm d}t$ because it is almost the same as that with low temporal resolution $\Delta t = 10 {\textrm d}t$ . Generally, the PINN-MB model correctly predicts the pressure around the swimming fish. In particular, in the high-pressure region near the head, the pressure distribution is consistent with that of the DNS. At the first instant $t=0.23$ , as shown in figure 8, a strong suction zone with low pressure is created when the middle of the body bends to the right side (Thandiackal & Lauder Reference Thandiackal and Lauder2020). A strong positive pressure is also consequently generated on the right side for pushing the water. PINN-MB can still correctly recover the pressure field at this instant. Visual inspection of the absolute error distribution in figure 8 shows that large errors are distributed mainly in the region of flow separation near the tail, where the flow is more complex than in the other regions. Figure 9 clearly shows the surface pressure coefficients extracted from the DNS and the predicted results. Note that PINN-MB can directly output the pressure at the boundaries without any interpolation. The predicted profiles with high temporal resolution $\Delta t = 2 {\textrm d}t$ and low temporal resolution $\Delta t = 10 {\textrm d}t$ show good agreement with the DNS results.

Figure 9.

Comparison of the surface pressure coefficients between the CFD results (thick solid curves) and the results predicted by PINN-MB with $\Delta t = 2 {\textrm d}t$ (thin solid curves) and $\Delta t = 10 {\textrm d}t$ (thin dashed curves). The different colours denote the right (red, R) and left (black, L) sides of the fish body. From left to right, the non-dimensional times are 0.23, 0.35 and 0.47, respectively.

The results of the quantitative error analysis of the pressure reconstruction results for the swimming fish are presented in figure 10. We first present the relative $\mathcal {l}_2$ norm errors for the whole field and boundary as a function of time at a high temporal resolution of $\Delta t = 10{\textrm d}t$ . Large error fluctuations of $10\, \%$ to $30\, \%$ are observed in this figure. The error reaches a maximum at $t=0.35$ and decreases to a minimum at $t=0.23$ and $0.47$ . This is because the relative error is normalised by the magnitude of the reference pressure (CFD result), which varies over time. As shown in figures 8 and 9, the pressure magnitude at $t=0.35$ is smaller than that at $t=0.23$ and $0.47$ . Therefore, the relative error is larger at $t=0.35$ . We also present the mean pressure error as a function of the time steps of training data in figure 10(b). The mean pressure errors for the fields and boundaries are $14\, \%$ and $20\, \%$ , respectively. The error profiles remain relatively flat across all the temporal resolutions tested, which is consistent with the results reported by Calicchia et al. (Reference Calicchia, Mittal, Seo and Ni2023). This result indicates that PINN-MB is insensitive to the temporal resolution of the training data for the embedding of the N–S equations.

Figure 10.

(a) Time series of the relative $\mathcal {l}_{2}$ norm error between the predicted pressure and the DNS pressure. (b) Error of the predicted pressure as a function of the time steps of the training data $\Delta t$ . Note that all the data, including the velocities and boundaries, are predicted at a time interval of $\Delta t = 2\textrm{d}t$ . The red and black curves denote the errors of fields and boundaries, respectively.

3.3. Example 3: 2-D pressure in the left ventricle

The primary function of the LV is to pump oxygenated blood into the aorta and then to every part of the body. Left ventricle dysfunction caused by structural injuries or pathological changes can induce pathological remodelling of the LV, which can lead to heart failure (HF) with significant morbidity and mortality (Stoll et al. Reference Stoll, Hess, Rodgers, Bissell, Dyverfeldt, Ebbers, Myerson, Carlhäll and Neubauer2019). Many studies have investigated cardiovascular imaging to quantitatively assess cardiac function in terms of morphology and structure. However, ventricular flow is altered in the early stages of remodelling (Stoll et al. Reference Stoll, Hess, Rodgers, Bissell, Dyverfeldt, Ebbers, Myerson, Carlhäll and Neubauer2019), which means that flow properties can be used to diagnose and prognose cardiac function (Pedrizzetti et al. 2014). The haemodynamic force, which corresponds to the global value of the intraventricular pressure gradient (IVPG), was suggested as a potential indicator for detecting early alterations in LV dysfunction and predicting disease outcome (Lapinskas et al. Reference Lapinskas, Pedrizzetti, Stoiber, Düngen, Edelmann, Pieske and Kelle2019; Mele et al. Reference Mele, Smarrazzo, Pedrizzetti, Capasso, Pepe, Severino, Luisi, Maglione and Ferrari2019; Stoll et al. Reference Stoll, Hess, Rodgers, Bissell, Dyverfeldt, Ebbers, Myerson, Carlhäll and Neubauer2019; Vallelonga et al. Reference Vallelonga, Airale, Tonti, Argulian, Milan, Narula and Pedrizzetti2021). Therefore, in vitro experiments were conducted to measure the velocity distribution in the LV and reconstruct the intraventricular pressure using the proposed PINN-MB method for the motion of the ventricular wall.

Figure 11 shows the set-up of the time-resolved 2-D PIV measurements. The pulsatile flow was simulated using an in-house LV simulator driven by a piston pump. The displacement of the piston was controlled by a user-defined signal, which is a conversion of the time evolution of the left ventricular volume extracted from Pedrizzetti & Perktold (Reference Pedrizzetti and Perktold2003). The 3-D-printed silicone LV model (refractive index of 1.41, thickness of 2.0 mm, end-systolic volume of 80 ml, base-to-apex length of 82 mm) was immersed in an acrylic hydraulic chamber filled with the same flow media. The systole and diastole of the LV were imitated by increasing and decreasing the pressure in the hydraulic chamber by using a piston pump. This hydraulic chamber was designed to have nine sides, as proposed by Saaid et al. (Reference Saaid2019), to facilitate 2-D or 3-D PIV measurements. During systole, blood flows from the LV to the aorta through the aortic valve, and then flows through the compliance chamber, resistance valve, fluid reservoir and left atrium before returning to the LV through the mitral valve. The compliance chamber downstream of the aortic valve is used to mimic arterial compliance. A resistance valve is used to adjust the systematic resistance. The fluid reservoir is open to the atmosphere to simulate the low pressure of the venous system.

Figure 11.

Set-up of the time-resolved 2-D-PIV measurements: (a) picture of the set-up, (b) schematic of the set-up. The native mitral and aortic valves, which were anatomised from a fresh porcine heart, were connected to a 3-D-printed silicone LV. As the piston pulls back, the aortic valve is closed due to flow deceleration, and the blood analogue flows through the opened mitral valve to fill the LV. As the piston moves forwards, the mitral valve is closed, and blood is pumped into the aorta, which is driven by rapidly increasing pressure. Note that a part of the flow field is sheltered by the papillary muscles.

Native mitral and aortic valves were anatomised from a fresh porcine heart weighing approximately 330 g, because the porcine heart has been used in preclinical trials and has similar haemodynamic performance to that of humans (Gooden et al. Reference Gooden, Hatoum, Zhang, Boudoulas and Dasi2022). In this experiment, a complete mitral valve consisting of the mitral annulus, leaflets, chordae tendineae and papillary muscles was sutured to a specially designed plate, and this plate was connected with a 3-D-printed silicone LV. The papillary muscles are attached to the LV wall with connecting rods, ensuring synchronous movement with the LV. However, the papillary muscles block the laser light, producing a shadow without particles, as shown in figure 11(b). A mixed solution ( $62\, \%$ water, $38\, \%$ glycerol by volume and 60 g/L NaCl by weight at a working temperature of $26\,^{\textrm {o}}\textrm {C}$ ) with a dynamic viscosity of 3.49 mPa $\cdot$ s, a density of 1173 kg/m $^3$ and a refractive index of 1.395 is used as blood analogue. The refractive index does not perfectly match that of silicone, and a better mixed solution was reported by Brindise et al. (Reference Brindise, Busse and Vlachos2018) and Chen & Dasi (Reference Chen and Dasi2023).

Because this work focused on pressure reconstruction from velocity fields, we provide only a brief introduction to the physiological conditions. The heart rate was 60 bpm and the stroke volume was 60 ml, this configuration yielded a cardiac output of 3.6 l/min. Figure 12 shows the transmitral pressure gradient (ventricular pressure–atrial pressure) and mitral flow rate during the cardiac cycle. In systole, the maximum pressure gradient is greater than 80 mmHg, and in diastole, the maximum instantaneous flow rate is approximately 8.3 l/min. The mean and maximum mitral inflow Reynolds numbers are 837 and 1930, respectively, computed as $Re=UD/\nu$ , where $U$ denotes the mean and maximum mitral velocities deduced from the cardiac output and the maximum instantaneous flow rate. The parameter $D$ denotes the diameter of the mitral valve which is equal to 26 mm in this case. The associated Womersley number of the mitral inflow, defined as $Wo = D/2\sqrt {{2\pi }/{T\nu }}$ , where $T$ is the time of the cardiac cycle, is approximately 17. The systematic vascular resistance (SVR) is approximately 777 dyn $\cdot$ s $\cdot$ cm $^{-5}$ , which is slightly less than that of a healthy human (800–1200 dyn $\cdot$ s $\cdot$ cm $^{-5}$ ).

Under these conditions, time-resolved 2-D PIV experiments were performed to measure the velocity field at the plane of the left ventricle outflow tract (LVOT). This plane is commonly used in clinical practice to evaluate blood flow and cardiac function. Moreover, the flow over the LVOT is approximately two-dimensional because this plane is a geometrically symmetric plane of the LV. The flow was seeded with 7 $\mu$ m fluorescent Rhodamine-B-coated particles at a density of 1500 kg/m $^3$ . The response frequency was $\sim$ 945 kHz, which means that the fluorescent particles can fully follow the flow motion. These fluorescent particles were excited using a double-pulsed laser (Beamtech, Vlite-Hi-50) with a wavelength of 527 nm at a repetition rate of 0.2 kHz. The laser beam was guided to the hydraulic chamber by using a laser arm and the thickness of the laser sheet was approximately 1.0 mm. Orange fluorescence with a longer wavelength of 590 nm was emitted and recorded by a high-speed CMOS camera (Lavision Imager HS CMOS cameras) equipped with a 100 mm/f2.8 lens (Tokina atx-i) with a high-pass filter. The image pairs were recorded at 200 Hz with a full image resolution of 2016 $\times$ 2016 pixels, and the time separation was set to 400 $\mu \textrm{s}$ to ensure that the maximum particle displacement was less than 8 pixels. This configuration yields a digital resolution of 20 pixels/mm. The velocity fields were calculated using in-house PIV software with a three-pass window deformation iterative multigrid scheme (Scarano Reference Scarano2002; Wang et al. Reference Wang, He and Wang2020a ). The interrogation window size was reduced to 32 $\times$ 32 pixels with $75\, \%$ overlap, and the velocity at the final pass was smoothed using Gaussian filtering with the same size as the interrogation window. This configuration yields a spatial resolution of 1.6 mm, and the vector step is 0.4 mm. The LV edge of each frame is detected according to the image intensity gradient and converted to physical coordinates based on the mapping function of the PIV data. Each edge contains 721 points with uniform distribution. In the present work, four cardiac cycles, consisting of 800 velocity fields with boundary data, were used to reconstruct the pressure fields.

Figure 12.

Flow rate (red curve) and transmitral pressure gradient (ventricular pressure–atrial pressure, blue curve) for the native mitral valve during the cardiac cycle. The time is normalised with respect to the period.

To recover the pressure fields, the boundary network consists of nine hidden layers and 64 neurons per hidden layer, and the flow network contains 13 hidden layers and 200 neurons per hidden layer. This configuration yields 517 491 trainable parameters. The location of the LV wall at the first time step is chosen as the initial input $X_b^0$ for the boundary network. Approximately $1.4\times 10^6$ velocity vectors are randomly selected from the PIV data and used as supervised data for the flow network, $1.0\times 10^6$ boundary data are randomly selected as supervised data for the boundary network and $1.1\times 10^6$ equation points with a uniform distribution are used to collect the residuals of the 2-D incompressible N–S equations. The dual-network model was trained for 10 000 epochs with 30 iterations per epoch. The corresponding loss functions are presented in figure 13. The training losses gradually converge to stable values after 1000 epochs. The final equation loss $\mathcal {L}_{eqns}$ , boundary condition loss $\mathcal {L}_{b.c.}$ , flow data loss $\mathcal {L}_{fdat}$ and boundary data loss $\mathcal {L}_{bdat}$ are approximately $2\times 10^{-3}$ , $4\times 10^{-4}$ , $2\times 10^{-3}$ and $3\times 10^{-5}$ , respectively. The equation loss is greater than that in the 3-D example as shown in figure 17(a), this is because of the two-dimensional assumption applied at the plane of LVOT. Figure 14 presents a series of instantaneous flow fields in three different phases (from left to right: early diastole, end diastole and systole) during the cardiac cycle. From top to bottom, these figures display the contour maps of the original PIV velocity magnitude and the velocity magnitude and pressure field predicted by PINN-MB. Note that the blank region within the original velocity is caused by the shadow of the papillary muscles as shown in figure 11(b). However, the proposed method has the capability to physically fill the gappy flow fields, the resulting fields are shown in figure 14. Because the ground truth pressure is unknown, we can only draw conclusions from a qualitative examination. In early diastole, when the mitral jet enters the LV, the inflow is decelerated by the increasing pressure gradient directed from the apex to the base. At the end of diastole, the pressure is high when the jet impinges onto the apex. During systole, the pressure is greater at the apex than at the base and the pressure accelerates the outgoing blood flow. However, the pressure gradient between the apex and base is much lower than that during diastole. Physically, the pressure presents a low level in the vortex ring centre induced by the mitral jet. This evolution of the IVPG is consistent with the clinical ultrasound results of a normal person (Pedrizzetti et al. 2014; Mele et al. Reference Mele, Smarrazzo, Pedrizzetti, Capasso, Pepe, Severino, Luisi, Maglione and Ferrari2019).

Figure 13.

Loss functions (equation loss $\mathcal {L}_{eqns}$ , boundary condition loss $\mathcal {L}_{b.c.}$ , flow data loss $\mathcal {L}_{fdat}$ and boundary data loss $\mathcal {L}_{bdat}$ ) as a function of epochs for the 2-D experimental case of LV.

Figure 14.

Instantaneous flow fields in the LV during different phases of the cardiac cycle: early diastole, end diastole and systole. From top to bottom, these figures present contour maps of the original PIV velocity magnitude, and the velocity magnitude and pressure fields predicted by the PINN-MB model. The predicted mesh is the same as the PIV mesh without increasing the spatial resolution. The unit of velocity is m s^-1 and the pressure is converted to mmHg.

Figure 15.

Haemodynamic force imposed on the blood during the cardiac cycle. The force is normalised by $\rho S g$ , where $\rho$ is the blood density, $S$ is the area enclosed by the LV boundary and $g$ denotes the gravitational acceleration. The time is normalised with respect to the period. The two components of the force, $F_x$ and $F_y$ , represent the forces in the septal-lateral and apical-basal directions, respectively.

Once the surface pressure is obtained, the haemodynamic force imposed on the blood can be calculated as the integral of the surface pressure distribution:

(3.4) \begin{equation} \mathbf {F}(t) = -\int _{S(t)}^{}p\mathbf {n}{\textrm d}S. \end{equation}

Here, $S(t)$ is the surface of the LV and $\mathbf {n}$ is the outward unit normal vector. Intraventricular blood flow is driven by the pressure gradient, and the haemodynamic force $\mathbf {F}(t)$ represents a global force exchanged between the blood volume and the LV wall (Pedrizzetti Reference Pedrizzetti2019; Vallelonga et al. Reference Vallelonga, Airale, Tonti, Argulian, Milan, Narula and Pedrizzetti2021). Compared with the conventional approach for the quantitative assessment of cardiac function, such as the ejection fraction, the use of haemodynamic force may provide insights into the interaction between tissue mechanics and fluid dynamics in the LV. Figure 15 displays the haemodynamic force in a period of the cardiac cycle. We note that this force is normalised with respect to the density, area (volume for the 3-D case) and gravitational acceleration for comparison (Pedrizzetti Reference Pedrizzetti2019; Vallelonga et al. Reference Vallelonga, Airale, Tonti, Argulian, Milan, Narula and Pedrizzetti2021). The two components of the force, $F_x$ and $F_y$ , represent the force in the septal-lateral and apical-basal directions, respectively. It is clear that the apical-basal force is much greater than the septal-lateral force. During systole, blood flow is accelerated towards the outflow direction because the haemodynamic force is directed from the apex to the base, as shown in figure 14. In early diastole, the force is negative and directed from the base to the apex, generating diastolic suction due to the relaxation of the LV. After this stage, the apical pressure gradually exceeds the basal pressure due to ventricular filling, which causes force to be directed from the apex to the base to decelerate the inflow. The haemodynamic force in diastole was greater than that in systole in our study; however, the magnitude of the force in these two phases was approximately the same as that reported in a previous numerical study by Pedrizzetti (Reference Pedrizzetti2019). This may because the cardiac output in our experiment (3.6 l/min) was lower than that of a normal person (5.0 l/min). The clinical assessment of intraventricular flow is beyond the scope of this article, please refer to Pedrizzetti et al. (2014), Bermejo et al. (Reference Bermejo, Martínez-Legazpi and del Álamo2015) and Vallelonga et al. (Reference Vallelonga, Airale, Tonti, Argulian, Milan, Narula and Pedrizzetti2021) for more information.

3.4. Example 4: 3-D pressure in the left ventricle

Figure 16.

Set-up of the time-resolved tomographic PIV measurements.

Thus far, the proposed method has been validated and investigated only for different types of 2-D flows. In this section, the feasibility of the proposed method for investigations of 3-D intraventricular flow is explored. Here as well, time-resolved tomographic PIV (Tomo-PIV) experiments were conducted based on the LV simulator described in § 3.3. Four high-speed cameras with 100 mm/f2.8 lenses were placed laterally around the nine-sided hydraulic chamber, as shown in figure 16. The view angle between two adjacent cameras was set to $36^{\textrm {o}}$ . A cylindrical laser with a diameter of 6 cm was expanded using a series of concave and convex lenses to illuminate all of the particles in the 3-D-printed silicone LV. Compared with the 2-D experiment, larger fluorescent Rhodamine-B coated particles with a diameter of 50 $\mu$ m were used in the 3-D experiment due to lower energy of the expanded laser. The response frequency was $\sim$ 18 kHz, which means that the fluorescent particles can still follow the flow motion. The flow conditions were the same as those in the 2-D experiment. The raw particle images were preprocessed to remove the background and normalise the particle intensity. The particle concentration after preprocessing was approximately 0.045 particles per pixel (ppp). The 3-D intensity distribution was reconstructed using the multiplicative algebraic reconstruction technique proposed by Elsinga et al. (Reference Elsinga, Scarano, Wieneke and van Oudheusden2006), based on the mapping function after self-calibration (Wieneke Reference Wieneke2008). The reconstructed volume was 70 $\times$ 70 $\times$ 50 mm $^3$ , containing 1400 $\times$ 1400 $\times$ 1000 voxels. The displacement was estimated using 3-D cross-correlation with a final interrogation volume of 48 $\times$ 48 $\times$ 48 voxels at $62.5\, \%$ overlap. This configuration yields 75 $\times$ 75 $\times$ 52 vectors with a spatial resolution of 2.4 mm and a vector pitch of 0.9 mm. Finally, all velocity fields were postprocessed using a divergence-free smoothing algorithm to clearly visualise the flow structures (Wang et al. Reference Wang, Gao, Wang, Wei, Li and Wang2016a ). The whole Tomo-PIV process was implemented using a validated in-house MATLAB code (Gao et al. Reference Gao, Wang and Wang2012; Wang et al. Reference Wang, Gao, Wei and Wang2016b ; Zhu et al. Reference Zhu, Wang, Wang and Wang2017; Tu et al. Reference Tu, Wang, Wang, Gao and Wei2022).

Three-dimensional three-component (3-D3C) measurements were obtained through Tomo-PIV; however, it is challenging to accurately and directly determine the wall of LV. In this work, our aim was to investigate the feasibility of the proposed method in the application to a 3-D moving boundary; therefore, an alternative method for obtaining a coarse boundary from the velocity is proposed. First, the mask data are identified as the region containing non-zero velocity (the velocity magnitude is larger than 0.01 m s^-1 in practice) because tracer particles are distributed only in the LV. Second, the boundary of the LV is extracted as an isosurface of the mask data. Third, the boundary points are matched by assuming that the vertices move along the wall-normal direction. Therefore, the walls at different times are described by a series of Lagrangian points with the same number and connectivity. In the present case, there are 4001 spatial points on the wall, and the time sequence is the same as that of the velocity data. The configuration of the PINN-MB model is as follows: the boundary network consists of nine hidden layers and 64 neurons per hidden layer, and the flow network contains 15 hidden layers and 200 neurons per hidden layer. This configuration yields 598 423 trainable parameters. The location of the LV wall at the first time step is chosen as the initial input $X_b^0$ for the boundary network. Approximately $6.6\times 10^6$ velocity vectors are randomly selected from three hundred Tomo-PIV fields with a time separation of 0.01 s (100 Hz) to serve as supervised data for the flow network, $2.0\times 10^6$ boundary data are randomly selected as the supervised data for the boundary network, $6.5\times 10^6$ equation points are used to collect the residuals of the 3-D incompressible N–S equations and $2.0\times 10^6$ points are used for the boundary conditions. The PINN-MB model was trained for 20 000 epochs with 30 iterations per epoch.

Figure 17.

(a) Loss functions (equation loss $\mathcal {L}_{eqns}$ , boundary condition loss $\mathcal {L}_{b.c.}$ , flow data loss $\mathcal {L}_{fdat}$ and boundary data loss $\mathcal {L}_{bdat}$ ) as a function of epoch for the 3-D experimental case of LV. (b) PDFs of the divergence for the tomographic PIV and PINN-MB results.

Figure 17(a) presents the corresponding loss functions as a function of the training epoch. The training losses gradually converge to stable values after 1000 epochs. The final equation loss $\mathcal {L}_{eqns}$ , boundary condition loss $\mathcal {L}_{b.c.}$ , flow data loss $\mathcal {L}_{fdat}$ and boundary data loss $\mathcal {L}_{bdat}$ are approximately $2\times 10^{-4}$ , $4\times 10^{-5}$ , $3\times 10^{-4}$ and $5\times 10^{-4}$ . Figure 17(b) shows the probability distribution function (PDF) of the divergence before and after optimisation. The maximum PDF of the PINN-MB model is approximately two orders of magnitude greater than that of the raw data, implying that the reconstructed velocities better satisfy the physical constraints. The instantaneous flow fields of the raw tomographic PIV (panels a–c) and the PINN-MB results (panels d–f) are displayed in figure 18, from left to right, the time separation is 0.05 s. The vortex structures are visualised as isosurfaces with $Q=500$ , where $Q$ denotes the vortex identification criterion estimated as the second invariant of the velocity gradient tensor and the contour map denotes the pressure distribution. Note that there are no pressure data in the raw tomographic PIV data. The large-scale vortex structures are clearly approximately the same in the raw data and the optimised results. Some small structures, which may be contamination by measurement noise, disappear for the correction of the N–S equations. The pressure fields of greatest concern are represented by the colourful contour map in figures 18(d)–18(f). The pressure is converted to mmHg for clinical applicability. At the first time instant when the mitral vortex enters the LV, two high-pressure regions exist at the front and back of the vortex structure. A low-pressure region is observed, which is consistent with the vortex core, to balance the centrifugal force and pressure force. The pressure gradient between the apex and base reaches a maximum when the vortex travels towards the middle of the LV. At this time, a vortex ring separates from the main vortex and continues to move ahead until it impinges on the apex of the LV, as shown in figure 18(f). The pressure acting on the wall of the LV can be directly obtained from the flow network; therefore, it is convenient to compute the haemodynamic force to describe cardiac function. From this 3-D example, we can conclude that it is feasible to reconstruct the pressure field in a scenario with a moving boundary using the proposed method.

Figure 18.

Instantaneous left ventricular flow fields of (a, b, c) the raw tomographic PIV data and (e, f, g) the present machine learning-based results. The surface pressure can be directly obtained from the flow network. The vortex structures are visualised using the Q criterion with $Q=500$ . From left to right, the time separation is 0.05 s.