From 8f4f49d650e420d5620a1550a4617bd564507e16 Mon Sep 17 00:00:00 2001
From: Lars Hubatsch <hubatsch@pks.mpg.de>
Date: Mon, 25 Nov 2019 16:41:39 +0100
Subject: [PATCH] Works nicely. Still missing proper boundary conditions that
take into account fluxes.
---
ternary_frap.m | 92 ++++++++------------------------------------------
1 file changed, 15 insertions(+), 77 deletions(-)
diff --git a/ternary_frap.m b/ternary_frap.m
index 76a120e..66c8867 100644
--- a/ternary_frap.m
+++ b/ternary_frap.m
@@ -4,74 +4,35 @@
% conditions are introduced. Integration of model via pdepe.
%% Solve pde
-% x = [linspace(0, 40, 100), linspace(40.4, 60, 50), linspace(60.4, 100, 100)];
-x = linspace(40, 60, 300);
-% x = linspace(0.5, 100.5, 1001);
-
-t = linspace(0, 500000, 50);
-% y = linspace(0.3000001, 1.2999999, 500);
-% x = 0.5*atanh(2*y-1.6)+50;
-a = -50;
-b = 0.5;
-% x_t = linspace(0.5, 101.5, 1011);
-% global phi_tot_g
-% phi_tot_g = phi_tot(x_t, b, a);
-% global der_phi
-% der_phi = gradient_analytical(x_t, a, b);
-% global lapl_phi
-% lapl_phi = laplacian_analytical(x_t, a, b);
-
-% opt = odeset('RelTol',1e-14, 'AbsTol', 1e-16,'MaxStep',1e-2);
-% sol = pdepe(0,@flory_pde, @flory_ic, @flory_bc,x,t, opt);
+x = linspace(49, 80, 700);
+t = linspace(0, 6, 1000);
sol = pdepe(0, @flory_hugg_a, @flory_ic, @flory_bc, x, t);
-
%% Plotting
figure(1); hold on;
-for i = 1:50
- cla; plot(x, phi_tot(x, -50, 1)); plot(x, sol(i, :)); pause(0.1);
+for i = 1:500
+ cla; xlim([49, 52]); ylim([0, 1.5]);
+ plot(x, phi_tot(x, -50, 0.25)); plot(x, sol(i, :)); pause(0.01);
end
%% Plot and check derivatives of pt
figure; hold on;
-% x = linspace(40, 60, 100);
-plot(x, phi_tot(x, -50, 0.5));% plot(x, gra_pt(x, -50, 0.5, 0.001)); plot(x, lap_pt(x, -50, 0.5, 0.001));
-% plot(x, gralap_pt(x, -50, 0.5, 0.001)); plot(x, laplap_pt(x, -50, 0.5, 0.001));
-plot(x, gradient_analytical(x, -50, 0.5));
-plot(x, laplacian_analytical(x, -50, 0.5));
-%%
-phi_tot(x, -50, 1)
+x = linspace(40, 60, 100);
+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)
% Solve with full ternary model. Analytical derivatives.
% pt ... phi_tot
% gra_a ... analytical gradient of phi_tot
-% lap_a ... analytical laplacian of phi_tot
-pt = @(x) phi_tot(x, -50, 1);
-gra_a = @(x) gradient_analytical(x, -50, 1);
-lap_a = @(x) laplacian_analytical(x, -50, 1);
-% gra_a = @(x) gra_pt(x, -50, 0.5, 0.0001);
-% lap_a = @(x) lap_pt(x, -50, 0.5, 0.0001);
-c = 1/(1.3-pt(x));
-f = dudx;
-s = u/(1.3-pt(x))*(lap_a(x)+(gra_a(x)/pt(x))^2-lap_a(x)/pt(x))...
- -dudx/(1.3-pt(x))/pt(x)*gra_a(x);
+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);
+s = 0;
end
-% function [c, f ,s] = flory_hugg_a(x, t, u, dudx)
-% % Solve with full ternary model.
-% global phi_tot_g
-% global der_phi
-% global lapl_phi
-%
-% i = round(10*x);
-% % disp(i)
-% c = 1/(1-phi_tot_g(i));
-% f = dudx;
-% s = u/(1-phi_tot_g(i))*(lapl_phi(i)+(der_phi(i)/phi_tot_g(i))^2-lapl_phi(i)/phi_tot_g(i))...
-% -dudx/(1-phi_tot_g(i))/phi_tot_g(i)*der_phi(i);
-% end
-
function u0 = flory_ic(x)
if x<50
u0 = 0.0;
@@ -79,6 +40,7 @@ function u0 = flory_ic(x)
u0 = 0.3;
end
% u0 = 0.3;
+% u0 = phi_tot(x, -50, 1);
end
function [pl,ql,pr,qr] = flory_bc(xl,ul,xr,ur,t)
@@ -92,30 +54,6 @@ function p = phi_tot(x, a, b)
p = (tanh(-(x+a)/b)+1)/2+0.3;
end
-function gpt = gra_pt(x, a, b, delta)
- gpt = (phi_tot(x+delta, a, b)-...
- phi_tot(x-delta, a, b))/(2*delta);
-end
-
-function lpt = lap_pt(x, a, b, delta)
- lpt = (gra_pt(x+delta, a, b, delta)-...
- gra_pt(x-delta, a, b, delta))/(2*delta);
-end
-%
-% function glpt = gralap_pt(x, a, b, delta)
-% glpt = (lap_pt(x+delta, a, b, delta)-...
-% lap_pt(x-delta, a, b, delta))/(2*delta);
-% end
-%
-% function llpt = laplap_pt(x, a, b, delta)
-% llpt = (gralap_pt(x+delta, a, b, delta)-...
-% gralap_pt(x-delta, a, b, delta))/(2*delta);
-% end
-
function grad = gradient_analytical(x, a, b)
grad = -(1-tanh(-(x+a)/b).^2)*1/b*0.5;
-end
-
-function lap = laplacian_analytical(x, a, b)
- lap = -2*tanh(-(x+a)/b).*(1-tanh(-(x+a)/b).^2)*1/b^2*0.5;
end
\ No newline at end of file
--
GitLab