From 19e95405b4d2af0366bd218ffc797494440c7716 Mon Sep 17 00:00:00 2001
From: Lars Hubatsch <hubatsch@pks.mpg.de>
Date: Tue, 26 Nov 2019 11:28:24 +0100
Subject: [PATCH] Boundary conditions, initial conditions, phi_tot are now all
parameterised.
---
ternary_frap.m | 42 ++++++++++++++++++++++++------------------
1 file changed, 24 insertions(+), 18 deletions(-)
diff --git a/ternary_frap.m b/ternary_frap.m
index 1119502..f36facb 100644
--- a/ternary_frap.m
+++ b/ternary_frap.m
@@ -2,16 +2,22 @@
% unbleached species. Model is assumed to be equilibrated
% (bleached+unbleached=const.=pt). Then bleached species initial
% conditions are introduced. Integration of model via pdepe.
-
+a = -50;
+b = 0.05;
+c = 1/1000000000;
+u0 = 0.1;
%% Solve pde
-x = linspace(49, 80, 700);
-t = linspace(0, 100, 1000);
-sol = pdepe(0, @flory_hugg_a, @flory_ic, @flory_bc, x, t);
+x = [linspace(49, 55, 600), linspace(55.01, 200, 300)];
+t = linspace(0, 100, 100);
+fh_ic = @(x) flory_ic(x, u0);
+fh_bc = @(xl, ul, xr, ur, t) flory_bc(xl, ul, xr, ur, t, u0);
+fh_pde = @(x, t, u, dudx) flory_hugg_pde(x, t, u, dudx, a, b, c, u0);
+sol = pdepe(0, fh_pde, fh_ic, fh_bc, x, t);
%% Plotting
figure(1); hold on;
-for i = 1:500
+for i = 1:100
cla; xlim([49, 53]); ylim([0, 1.5]);
- plot(x, phi_tot(x, -50, 0.25)); plot(x, sol(i, :)); pause(0.01);
+ plot(x, phi_tot(x, a, b, u0)); plot(x, sol(i, :)); pause(0.01);
end
%% Plot and check derivatives of pt
@@ -21,41 +27,41 @@ plot(x, phi_tot(x, -50, 0.25));
plot(x, gradient_analytical(x, -50, 0.25));
%% Function definitions for pde solver
-function [c, f ,s] = flory_hugg_a(x, t, u, dudx)
+function [c, f ,s] = flory_hugg_pde(x, t, u, dudx, a, b, c_p, u0)
% Solve with full ternary model. Analytical derivatives.
% pt ... phi_tot
% gra_a ... analytical gradient of phi_tot
-pt = phi_tot(x, -50, 0.25);
-gra_a = gradient_analytical(x, -50, 0.25);
-c = 1/100;
-f = (1.3-pt)/pt*(pt*dudx-u*gra_a);
+pt = phi_tot(x, a, b, u0);
+gra_a = gradient_analytical(x, a, b);
+c = c_p;
+f = ((u0+1)-pt)/pt*(pt*dudx-u*gra_a);
s = 0;
end
-function u0 = flory_ic(x)
+function u = flory_ic(x, u0)
if x<50
- u0 = 0.0;
+ u = 0.0;
else
- u0 = 0.3;
+ u = u0;
end
% u0 = 0.3;
% u0 = phi_tot(x, -50, 1);
end
-function [pl,ql,pr,qr] = flory_bc(xl,ul,xr,ur,t)
+function [pl,ql,pr,qr] = flory_bc(xl, ul, xr, ur, t, u0)
pl = 0;
ql = 1;
% % No flux
% pr = 0;%ur - 0.01;
% qr = 1;
% Dirichlet BC
- pr = ur- 0.3;
+ pr = ur - u0;
qr = 0;
end
-function p = phi_tot(x, a, b)
- p = (tanh(-(x+a)/b)+1)/2+0.3;
+function p = phi_tot(x, a, b, u0)
+ p = (tanh(-(x+a)/b)+1)/2+u0;
end
function grad = gradient_analytical(x, a, b)
--
GitLab