lead_advec_diff.Rmd

---
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()
```