From b52412d8e5e2ce2794d968b775613092862818aa Mon Sep 17 00:00:00 2001
From: Lars Hubatsch <hubatsch@pks.mpg.de>
Date: Mon, 30 Dec 2019 12:13:18 +0100
Subject: [PATCH] Generalizing script. Optimizing mesh.
---
ternary_frap.m | 53 ++++++++++++++++++++++++++++++++------------------
1 file changed, 34 insertions(+), 19 deletions(-)
diff --git a/ternary_frap.m b/ternary_frap.m
index 044c34a..1d09d6c 100644
--- a/ternary_frap.m
+++ b/ternary_frap.m
@@ -4,39 +4,50 @@
% conditions are introduced. Integration of model via pdepe.
a = -50;
-b = 0.01;
+b = 0.025;
c = 1/100000;
-u0 = 0.3;
+u0 = 0.05;
+e = 0.4;
%% Solve pde
-x = [linspace(49, 55, 600), linspace(55.01, 200, 300)];
-t = linspace(0, 0.001, 1000);
+x = [linspace(49.0, 51, 3000), linspace(51.01, 200, 300)];
+t = linspace(0.001, 1, 1000);
fh_ic = @(x) flory_ic(x, a, 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);
+fh_pde = @(x, t, u, dudx) flory_hugg_pde(x, t, u, dudx, a, b, e, c, u0);
sol = pdepe(0, fh_pde, fh_ic, fh_bc, x, t);
%% Plotting
figure(1); hold on;
for i = 1:length(t)
- cla; xlim([49, 53]); ylim([0, 1.5]);
- plot(x, phi_tot(x, a, b, u0)); plot(x, sol(i, :)); pause(0.03);
+% i = 800;
+ cla; xlim([49, 53]); ylim([0, 0.7]);
+ ax = gca;
+ ax.FontSize = 16;
+ xlabel('position'); ylabel('volume fraction');
+ plot(x, phi_tot(x, a, b, e, u0), 'LineWidth', 2, 'LineStyle', '--');
+ plot(x, sol(i, :), 'LineWidth', 2); pause(0.06);
+% % print([num2str(i),'.png'],'-dpng')
end
%% Plot and check derivatives of pt
figure; hold on;
-x = linspace(40, 60, 100);
-plot(x, phi_tot(x, -50, 0.25));
-plot(x, gradient_analytical(x, -50, 0.25));
-
+x = linspace(40, 60, 10000);
+plot(x, phi_tot(x, a, b, e, u0));
+plot(x, gradient_analytical(x, a, b, e));
+plot(x(1:end-1)+mean(diff(x))/2, ...
+ diff(phi_tot(x, a, b, e, u0)/mean(diff(x))));
+plot(x, gamma0(x, a, b, e));
+max(gamma0(x, a, b, e))/min(gamma0(x, a, b, e))
%% Function definitions for pde solver
-function [c, f ,s] = flory_hugg_pde(x, t, u, dudx, a, b, c_p, u0)
+function [c, f ,s] = flory_hugg_pde(x, t, u, dudx, a, b, e, c_p, u0)
% Solve with full ternary model. Analytical derivatives.
% pt ... phi_tot
% gra_a ... analytical gradient of phi_tot
-pt = phi_tot(x, a, b, u0);
-gra_a = gradient_analytical(x, a, b);
+pt = phi_tot(x, a, b, e, u0);
+gra_a = gradient_analytical(x, a, b, e);
+g0 = gamma0(x, a, b, e);
c = c_p;
-f = ((u0+1.3)-pt)/pt*(pt*dudx-u*gra_a);
+f = g0*(1-pt)/pt*(pt*dudx-u*gra_a);
s = 0;
end
@@ -61,10 +72,14 @@ function [pl,ql,pr,qr] = flory_bc(xl, ul, xr, ur, t, u0)
qr = 0;
end
-function p = phi_tot(x, a, b, u0)
- p = (tanh(-(x+a)/b)+1)/2+u0;
+function g0 = gamma0(x, a, b, e)
+ g0 = 10*e*(tanh((x+a)/b)+1)/2+0.001;
+end
+
+function p = phi_tot(x, a, b, e, u0)
+ p = e*(tanh(-(x+a)/b)+1)/2+u0;
end
-function grad = gradient_analytical(x, a, b)
- grad = -(1-tanh(-(x+a)/b).^2)*1/b*0.5;
+function grad = gradient_analytical(x, a, b, e)
+ grad = -e*(1-tanh(-(x+a)/b).^2)/b*0.5;
end
--
GitLab