From 6bea983ae16470aa250bef3ebf39692dd3f11373 Mon Sep 17 00:00:00 2001
From: Lars Hubatsch <hubatsch@pks.mpg.de>
Date: Mon, 18 Nov 2019 19:06:59 +0100
Subject: [PATCH] Model doesn't work yet. Mass conservation doesn't work. Steep
gradient at interface.
---
ternary_frap.m | 87 ++++++++++++++++++++++++++++++++++++++++++++++++++
1 file changed, 87 insertions(+)
create mode 100644 ternary_frap.m
diff --git a/ternary_frap.m b/ternary_frap.m
new file mode 100644
index 0000000..3e7b8ac
--- /dev/null
+++ b/ternary_frap.m
@@ -0,0 +1,87 @@
+% Numerical solution of ternary FRAP model with solvent, bleached and
+% unbleached species. Model is assumed to be equilibrated
+% (bleached+unbleached=const.=pt). Then bleached species initial
+% conditions are introduced. Integration of model via pdepe.
+
+%% Solve pde
+x = linspace(0, 40, 500);
+t = linspace(0, 5000, 50);
+% opt = odeset('RelTol',1e-14, 'AbsTol', 1e-16,'MaxStep',1e-2);
+% sol = pdepe(0,@flory_pde, @flory_ic, @flory_bc,x,t, opt);
+sol = pdepe(0,@flory_hugg_pde, @flory_ic, @flory_bc,x,t);
+
+%% Plotting
+figure(1); hold on;
+for i = 1:50
+ cla; plot(x, pt(x)); plot(x, sol(i, :)); pause(0.01);
+end
+
+%% Plot and check derivatives of pt
+figure; hold on;
+plot(x, pt(x)); plot(x, gra_pt(x, 0.001)); plot(x, lap_pt(x, 0.001));
+plot(x, gralap_pt(x, 0.001)); plot(x, laplap_pt(x, 0.001));
+
+%% Function definitions for pde solver
+function [c, f ,s] = flory_hugg_pde(x, t, u, dudx)
+% Solve with full ternary model.
+X = 2.2;
+K = 1;
+c = 1/(pt(x)+1);
+f = dudx;
+s = u*(lap_pt(x, 0.001)-2*X*(lap_pt(x, 0.001)-...
+ gra_pt(x, 0.001)^2-pt(x)*lap_pt(x, 0.001))+...
+ K*(-laplap_pt(x, 0.001)+...
+ gra_pt(x, 0.001)*gralap_pt(x, 0.001)+...
+ pt(x)*laplap_pt(x, 0.001)))/(pt(x)+1)+...
+ dudx*(-2*X*(gra_pt(x, 0.001)-pt(x)*gra_pt(x, 0.001))-...
+ K*gralap_pt(x, 0.001)+K*pt(x)*gralap_pt(x, 0.001))/...
+ (pt(x)+1);
+end
+
+function u0 = flory_ic(x)
+ if x<5
+ u0 = 0.0;
+ else
+ u0 = 0.1;
+ end
+% u0 = pt(x)*0.5;
+end
+
+function [pl,ql,pr,qr] = flory_bc(xl,ul,xr,ur,t)
+ pl = 0;
+ ql = 1;
+ pr = 0;%ur - 0.01;
+ qr = 1;
+end
+
+function [c, f ,s] = flory_pde(x, t, u, dudx)
+% Tentatively solve system where binary mixture is assumed, with variable
+% '1' (now called pt). Basically pt is fixed and now phi_u and
+% phi_b can swap places but nothing else.
+c = 1/pt(x);
+f = dudx;
+s = -u*lap_pt(x)/pt(x);
+end
+
+function p = pt(x)
+ p = (tanh(-(x-5)*2)+1)/2+0.1;
+end
+
+function gpt = gra_pt(x, delta)
+ gpt = (pt(x+delta)-pt(x-delta))/(2*delta);
+end
+
+function lpt = lap_pt(x, delta)
+ lpt = (gra_pt(x+delta, delta)-...
+ gra_pt(x-delta, delta))/(2*delta);
+end
+
+function glpt = gralap_pt(x, delta)
+ glpt = (lap_pt(x+delta, delta)-...
+ lap_pt(x-delta, delta))/(2*delta);
+end
+
+function llpt = laplap_pt(x, delta)
+ llpt = (gralap_pt(x+delta, delta)-...
+ gralap_pt(x-delta, delta))/(2*delta);
+end
\ No newline at end of file
--
GitLab