Created
March 7, 2019 20:45
-
-
Save paleolimbot/dd49df2bc14f3a6a9bef6ab4bf1f27bb to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
--- | |
title: "R Notebook" | |
output: html_notebook | |
--- | |
```{r setup} | |
library(tidyverse) | |
``` | |
$$ | |
\frac{dC_s}{dt} = - | |
$$ | |
```{r} | |
radial_flow_velocity <- function(radius, radius_max = 1.9 / 2) { | |
# 0 - > 0 | |
# radius / radius_max -> 1 | |
radius / radius_max | |
} | |
# output: ug / L / s | |
delta_concentration <- function( | |
average_flow_velocity = 0.8, # cm / s | |
concentration = 70, # ug/L | |
diffusion_radial = 2, # cm^2 / s | |
diffusion_long = 2, # cm^2 / s | |
delta_x = 0.1, # cm | |
delta_radius = 0.2, # cm | |
radius = 0.5, # cm, | |
concentration_grad_x = 0, # ug / L / cm | |
concentration_grad_rad = 0 # ug / L / cm | |
) { | |
radial_comp <- -delta_x * average_flow_velocity * radial_flow_velocity(radius) * concentration | |
advection <- delta_x * diffusion_long * concentration_grad_x | |
radial_diff <- 1 / radius * delta_radius * radius * diffusion_radial * concentration_grad_rad | |
radial_comp + advection + radial_diff | |
} | |
``` | |
Test the responses | |
```{r} | |
# pipe wall -> middle | |
tibble( | |
r = seq(0, 0.9, length.out = 10), | |
delta_conc = delta_concentration(radius = r) | |
) | |
``` | |
2d model | |
```{r} | |
delta_t <- 1 | |
model_initial <- expand.grid( | |
radius = seq(0, 1.9 / 2, length.out = 10), x = seq(0, 100, length.out = 10) | |
) %>% | |
# initial boundary condition | |
mutate( | |
time = 0, | |
concenteration = 0, | |
average_flow_velocity = 0.8, # cm / s | |
concentration = 70, # ug/L | |
diffusion_radial = 2, # cm^2 / s | |
diffusion_long = 2, # cm^2 / s | |
delta_x = 0.1, # cm | |
delta_radius = 0.2, # cm | |
concentration_grad_x = 0, # ug / L / cm | |
concentration_grad_rad = 0, # ug / L / cm, | |
delta_concentration_time = 0 | |
) | |
model_states <- rep(list(model_initial), 100) | |
for(i in seq_along(model_states)[-1]) { | |
previous_state <- model_states[[i - 1]] | |
state <- model_states[[i]] | |
state$time <- (i - 1) * delta_t | |
# calculate some stuff, like concentration_grad_x, concentration_grad_rad | |
# based on this state and previous state | |
new_state <- state %>% | |
group_by(radius) %>% | |
arrange(x) %>% | |
mutate(concentration_grad_x = concentration - lag(concentration)) %>% | |
group_by(x) %>% | |
arrange(radius) %>% | |
mutate(concentration_grad_rad = concentration - lag(concentration)) %>% | |
ungroup() %>% | |
# do calculation! | |
mutate( | |
delta_concentration_time = delta_concentration( | |
radius = radius | |
), | |
concentration = concentration + delta_concentration_time * delta_t | |
) | |
model_states[[i]] <- new_state | |
} | |
final_df <- bind_rows(model_states) | |
final_df %>% | |
group_by(time) %>% | |
filter(radius == max(radius)) %>% | |
ggplot(aes(time, concentration)) + | |
geom_line() | |
``` |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment