Newer
Older
Lars Hubatsch
committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
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