<h2>Preface</h2>
<p>I’ve realized that all my <a rel="nofollow" href="https://incidental-ideas.org/2017/05/25/regular-expression-treemaps-to-visualize-emergency-department-visits/" target="_blank">previous posts</a> were quite substantial in length and took quite a long time to create them. From this point forward I’ll be generating posts of shorter length (partially for my sanity and more for my impulsivity with ideas). A few of these posts won’t be public health related (like this one) but they will be more of a personal exploration and application of different tools and packages in R. This post will be exploring the <strong><a rel="nofollow" href="https://cran.r-project.org/web/packages/animation/index.html" target="_blank">animation</a></strong> and <strong><a rel="nofollow" href="https://github.com/dgrtwo/gganimate" target="_blank">gganimate</a></strong> package to create some (hopefully) cool animated visualizations.</p>
<h2>Laminar Flow</h2>
<p>A simple definition of <strong><a rel="nofollow" href="https://en.wikipedia.org/wiki/Laminar_flow" target="_blank">laminar flow</a></strong> is when a fluid flows through a pipe in parallel layers with no disruptions between these layers. Turbulent flow occurs when there is mixing or any disruption between these layers.  Laminar flow is a simple analysis of fluid dynamics that considers the <strong>viscosity</strong> (“thickness/gooeyness”) of a liquid. This model is based on a number of assumptions where the fluid flows through a straight cylindrical pipe of fixed diameter, there is no acceleration (only a balance between pressure and viscous/shear forces) and gravitational forces are ignored. The walls of the pipe exhibit the maximum <a rel="nofollow" href="https://en.wikipedia.org/wiki/Shear_force" target="_blank">shearing force</a> and the axial centre of the pipe has zero shear. This results in a parabolic <a rel="nofollow" href="http://hyperphysics.phy-astr.gsu.edu/hbase/pfric2.html" target="_blank">velocity profile</a> where the velocity is at it’s maximum at centre of the pipe.</p>
<figure data-shortcode="caption" class="wp-caption"><img data-orig-file="https://eugejjoh.files.wordpress.com/2018/04/sketch-laminar-flow-turbulent-flow.png" data-orig-size="842,768" data-medium-file="https://eugejjoh.files.wordpress.com/2018/04/sketch-laminar-flow-turbulent-flow.png?w=300&amp;h=274&amp;fit=300%2C274&amp;resize=300%2C274" data-attachment-id="1961" data-image-meta="{" data-image-description="" src="https://eugejjoh.files.wordpress.com/2018/04/sketch-laminar-flow-turbulent-flow.png?w=300&h=274&fit=300%2C274&resize=300%2C274" data-permalink="https://incidental-ideas.org/2018/04/05/laminar-flow/sketch-laminar-flow-turbulent-flow/" data-large-file="https://eugejjoh.files.wordpress.com/2018/04/sketch-laminar-flow-turbulent-flow.png?w=736" alt="sketch-laminar-flow-turbulent-flow" data-image-title="sketch-laminar-flow-turbulent-flow" sizes="(max-width: 300px) 100vw, 300px" width="300" data-comments-opened="1" class="wp-image-1961" data-recalc-dims="1" height="274"><figcaption class="wp-caption-text"><b>© </b><a rel="nofollow" href="https://www.cfdsupport.com/OpenFOAM-Training-by-CFD-Support/node275.html" target="_blank">CFD support</a></figcaption></figure>
<h2>MATLAB to R?</h2>
<p>First you might wonder, why is this post randomly on laminar flow? A short trip down memory lane while organizing some <a rel="nofollow" href="https://www.schulich.uwo.ca/biophysics/undergraduate/course_information.html" target="_blank">undergraduate course</a> folders resulted in my rediscovery of some lecture notes and code related to <a rel="nofollow" href="https://en.wikipedia.org/wiki/Medical_physics#Medical_biophysics_and_biomedical_physics" target="_blank">medical biophysics</a> and <a rel="nofollow" href="https://en.wikipedia.org/wiki/Hemodynamics" target="_blank">hemodynamics</a>. Not having used <a rel="nofollow" href="https://en.wikipedia.org/wiki/MATLAB" target="_blank">MATLAB</a> in over 4 years, <strong>I was curious whether I could translate some of the MATLAB code into something similar into R</strong>. Some points on the <a rel="nofollow" href="https://www.mathworks.com/discovery/matlab-vs-r.html" target="_blank">differences</a>, MATLAB is a proprietary numerical computing programming software based on matrix algorithms while R is a open-source programming language centered on statistical analyses. There are some StackOverflow <a rel="nofollow" href="https://stackoverflow.com/questions/1738087/what-can-matlab-do-that-r-cannot-do" target="_blank">discussions</a> on these comparisons if you want to explore this further. I also discovered the <a rel="nofollow" href="https://cran.r-project.org/web/packages/pracma/index.html" target="_blank"><strong>pracma</strong> package</a> for numerical analysis that uses some MATLAB function names (#tbt). Now to the actual R code!</p>
<h2>Initial Conditions</h2>
<p>First we’ll load the relevant packages. I’ve commented in the code lines on some basic descriptions on each package. Just a small note that some of these packages need to be downloaded from GitHub repositories using the <a rel="nofollow" href="https://cran.r-project.org/web/packages/devtools/index.html" target="_blank"><strong>devtools</strong> package</a>.</p>
<pre>
# Packages ####
library(pracma) #load practical math package
library(ggplot2) #data viz
library(scales) #plot scales
library(dplyr) #data wrangling
library(viridis) #data viz palette
library(animation) #animation
library(gganimate) #ggplot animation
</pre>
<p>Next, for any type of model you need to first set the <strong>initial conditions</strong> or parameters. In this case, we’ll specify the conditions for the pressure difference at the beginning at the end of the pipe, the viscous force, the dimensions of the pipe (radius and length).</p>
<pre>
## Initial Conditions ####
# Pressure difference (delta P)
Pi &lt;- 10 #initial pressure [Pa]
Pf &lt;- 50 #final pressure [Pa]

## Boundary Conditions ####
# Pipe and fluid characteristics
mu &lt;- 8.94e-4 #viscosity for water. units: [Pa*s] at temp = 25 C
R &lt;- 0.1 #radius of pipe [m]
l &lt;- 4 #length of pipe [m]

N &lt;- 100 #radial resolution
r &lt;- seq(-R,R,length.out=N) #diameter of pipe
S &lt;- pi*(R^2) #cross-sectional area of pipe
</pre>
<p>Since we're also playing with the <a rel="nofollow" href="https://cran.r-project.org/web/packages/animation/index.html" target="_blank"><strong>animation</strong> package</a>, we’ll also input a vector of time units to loop and animate the visualization.</p>
<pre>
# time dependent
time &lt;- seq(0,19) #length 20
Pit &lt;- 0.2*time*Pi
Pft &lt;- 0.2*time*Pf
Pdiff &lt;- abs(Pft-Pit) #pressure difference over time
</pre>
<h2>Data Setup</h2>
<p>Now what we’ll do is run a loop over the time sequence we’ve specified (20 arbitrary units of time, seconds for example) to create the output data. The pressure difference will be increased as a function of time, simulating increased flow of the fluid (water in this case) over time. The equation in the code below is the flow velocity as a function of the pipe radius derived from <strong><a rel="nofollow" href="https://en.wikipedia.org/wiki/Hagen%E2%80%93Poiseuille_equation" target="_blank">Hagen-Poiseuille’s Law</a></strong> specified below. So after running <strong><a rel="nofollow" href="https://www.datacamp.com/community/tutorials/r-tutorial-apply-family?tap_a=5644-dce66f&tap_s=10907-287229" target="_blank">lapply()</a></strong> over the pressure differences, we then massage the list output into a data frame – first by “unlisting” the list elements, creating a matrix object where each column represents a different time point, and coercing the matrix into a data frame (and as <a rel="nofollow" href="https://simplystatistics.org/2015/07/24/stringsasfactors-an-unauthorized-biography/" target="_blank">always remembering</a> to use <b>stringsAsFactors = FALSE</b>).</p>
<p>Hagen-Poiseuille’s Law<br>
Q = πR<sup>4</sup>Δp / 8μl</p>
<pre>
# Create Data Viz ####
# plot shows V (average speed) and r (pipe radius)
# ((Pdiff)*(R^2)/(4*mu*l))*(1-(r^2)/R^2) #eqn for V
out1 &lt;- lapply(seq_along(Pdiff), function(x) {
((Pdiff[x])*(R^2)/(4*mu*l))*(1-(r^2)/R^2)
})

out1_df &lt;- data.frame(matrix(unlist(out1),ncol=length(time)),stringsAsFactors = FALSE)
names(out1_df) &lt;- seq_along(names(out1_df))
</pre>
<p>Now we have a data frame with each column containing the <strong>velocity of the fluid as a function of the radius</strong>. There are twenty time periods and 100 rows representing the resolution of the x-axis (e.g. 100 ticks). The problem now is that this data frame is not in a <a rel="nofollow" href="http://ggplot2.tidyverse.org/reference/" target="_blank"><strong>ggplot2</strong></a> friendly format. We’ll have to use the <a rel="nofollow" href="https://cran.r-project.org/web/packages/tidyr/index.html" target="_blank"><strong>tidyr</strong> package</a> to manipulate the data. The code below <a rel="nofollow" href="http://garrettgman.github.io/tidying/" target="_blank">“gathers”</a> the columns so they become grouped elements in one column and then we use the <a rel="nofollow" href="https://cran.r-project.org/web/packages/dplyr/index.html" target="_blank"><strong>dplyr</strong> package</a> to add the radius values as new column, coerce the “time” column into an integer type, and add the flowrate (a function of the velocity profile and pipe cross-sectional area) as another column. A note here that when the <strong>rbind()</strong> function is used, if the length/# of rows of the object’s (<strong>r</strong>) being “binded” is a multiple of the other object’s (<strong>df_t</strong>) length/# of rows, it will repeat based on the multiple. Since we looped based on the the radial resolution (100) specified before, <strong>rbind()</strong> will fill the column with 20 multiples (20*100) of the radius ticks.</p>
<pre>
df_t &lt;- tidyr::gather(out1_df,key=&quot;time&quot;,value=&quot;V&quot;)

df_t2 %
 mutate(time = as.integer(time)) %>% #change to integer type
 mutate(Qv = V/(1-(r^2)/R^2)) #add volumetric flowrate (Qv)
</pre>
<h2>Data Visualization</h2>
<p>Now we have our data in a nice format to visualize and loop over to create our GIF. We specify the data frame we want to use, map the aesthetics and specifically the <b>frame</b> and <b>cumulative</b> arguments. The <b>frame</b> aesthetic tells the <strong><a rel="nofollow" href="https://github.com/dgrtwo/gganimate" target="_blank">gganimate</a></strong><a rel="nofollow" href="https://github.com/dgrtwo/gganimate" target="_blank"> package</a><strong><a rel="nofollow" href="https://github.com/dgrtwo/gganimate" target="_blank"></a></strong> to loop over the values defined in <b>frame</b> or time in our case. We’ll use the amazing <a rel="nofollow" href="https://cran.r-project.org/web/packages/viridis/index.html" target="_blank"><strong>viridis</strong> package</a> to add some time colours, use the built-in <strong>ggplot2</strong> dark theme, tinker with the text sizes, and provide a good title and axis labels (yay for plotting best practices). The represented here displays the flow rate through the pipe going from purple (low flowrate) to yellow (high flowrate).</p>
<pre>
p1 &lt;- ggplot(data=df_t2, aes(y=r,x=V,col=Qv,frame=time,cumulative=TRUE)) +
 geom_path() + geom_point(alpha=0.05) +
 scale_colour_viridis(option = &quot;C&quot;, discrete = FALSE) +
 theme_dark() +
 labs(title = &quot;Hagen-Poiseuille's Law: Velocity Profile t =&quot;,
 y = &quot;Pipe Radius (m)&quot;, x = &quot;Velocity (m/s)&quot;)

ggsave(p1,filename = &quot;Vprofile1.png&quot;,width = 250, height=100,units = &quot;mm&quot;)
</pre>
<p><img data-orig-file="https://eugejjoh.files.wordpress.com/2018/04/vprofile1.png?w=456" data-orig-size="2952,1181" data-medium-file="https://eugejjoh.files.wordpress.com/2018/04/vprofile1.png?w=456?w=300" data-attachment-id="1964" data-image-meta="{" data-image-description="" src="https://eugejjoh.files.wordpress.com/2018/04/vprofile1.png?w=456" data-permalink="https://incidental-ideas.org/2018/04/05/laminar-flow/vprofile1-3/" data-large-file="https://eugejjoh.files.wordpress.com/2018/04/vprofile1.png?w=456?w=736" alt="Vprofile1" data-image-title="Vprofile1" sizes="(max-width: 736px) 100vw, 736px" data-comments-opened="1" class="wp-image-1964" data-recalc-dims="1"></p>
<p>The title might look strange but the empty “t = ” will be useful to indicate which time point each frame describes when we create the GIF (look at the <b>title_frame</b> argument in the <a rel="nofollow" href="https://www.rdocumentation.org/packages/gganimate/versions/0.1.0.9000?tap_a=5644-dce66f&tap_s=10907-287229" target="_blank"><strong>gganimate()</strong></a> function below. This will add the value provided to the frame aesthetic and include it into the plot title. We also specify the intervals for each frame (0.2 seconds) and the dimensions of the output.</p>
<pre>
gganimate(p1,filename = &quot;Vprofile1.gif&quot;, title_frame = TRUE,
 interval=0.2, ani.width = 900, ani.height = 350)
</pre>
<p><img data-orig-file="https://eugejjoh.files.wordpress.com/2018/04/vprofile11.gif?w=456" data-orig-size="600,300" data-medium-file="https://eugejjoh.files.wordpress.com/2018/04/vprofile11.gif?w=456?w=300" data-attachment-id="1959" data-image-meta="{" data-image-description="" src="https://eugejjoh.files.wordpress.com/2018/04/vprofile11.gif?w=456" data-permalink="https://incidental-ideas.org/2018/04/05/laminar-flow/vprofile1-2/" data-large-file="https://eugejjoh.files.wordpress.com/2018/04/vprofile11.gif?w=456?w=600" alt="Vprofile1.gif" data-image-title="Vprofile1" data-comments-opened="1" class="wp-image-1959" data-recalc-dims="1"></p>
<p>Now we have this nice animation visualizing the laminar flow velocity profile of a pipe at different time points where the pressure difference increases, in turn increasing the maximum velocity of the fluid. Another way to see it, is that at time 0, there is no pressure difference, therefore no flow. But once the pressure difference increases (decreased lower pressure at the end of pipe e.g. increased flow of fluid out of the end of the pipe) the velocity profile or flow increases.</p>
<p>The entire script containing the code is found in my <a rel="nofollow" href="https://github.com/eugejoh/explore_ani" target="_blank">GitHub repo</a> if you would want to fork it or tinker with it yourself. Again any feedback or comments are welcome, please check out my <a rel="nofollow" href="https://incidental-ideas.org/2017/12/09/watermain-breaks-in-the-city-of-toronto/" target="_blank">other posts</a>, and hope to write another post soon!</p>
	<img alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/6v6DPG7TTJg" height="1">