# 1D simulation diffusion term, polar coordinates, moving boundary condition

 User Name Remember Me Password
 Register Blogs Members List Search Today's Posts Mark Forums Read

 LinkBack Thread Tools Search this Thread Display Modes
November 1, 2022, 09:03
1D simulation diffusion term, polar coordinates, moving boundary condition
#1
New Member

Andrea
Join Date: Oct 2022
Posts: 3
Rep Power: 3
I am trying to integrate the system I am showing in the page 3 of the attached pdf but I am having some problems with Cw and the radius because they should reach a minimum and a maximum respectively but in my simulation they go towards +inf and -inf, and Cw is decreasing too fast. Cw0, K, Cm are parameters, I also attach the whole paper if something is not clear enough. (I am using matlab). I am using forward euler and finite difference method. It is probably better to integrate with a variable space step but I cant figure out how to do it, so I just made a big enough domain to allow the radius to increase. I have started recently cfd so be kind ahahah.

Code:
```clc; clear all; close all
%% DATA
K = 1;
Diff = 4e-9;
Cw0 = 55e-3; %mol/m3
r0 = 1e-3;
L = 12e-3;
N = 500;
tf = 150;
Cm = 30e-3; %mol/m3
Cinf = 0;
%% preprocessing
h = L/N;
grid1D = linspace(0,L,N+1);
%% SOLUTION
index = find(grid1D >= r0);
dt = 0.005;
tsteps = tf/dt;
%initial sol
C = [Cw0.*ones(1,index(1)-1)./K,Cm*ones(1,index(end)-index(1))]';
%mem allocation for plot
Cplot = zeros(size(grid1D));
%loop
for t = 1:tsteps
C0 = C;
% bc 6: eq 6 integration
if t ~= 1
indexrd = find(grid1D >= rd);
for counter_bc6 = index(1):indexrd(1) - 1
int6(counter_bc6) = (-4*pi*C(counter_bc6)*grid1D(counter_bc6)^2 -4*pi*C(counter_bc6 + 1)*grid1D(counter_bc6+1)^2)*h/2;
end
Cw = sum(int6)/(4/3*pi*r0^3) + Cw0;
else
Cw = Cw0;
end
% eq 3
C(index(1)) = Cw/K;
% eq 5
if t == 1
rd = r0 + h;
else
d = -Cm + 8*Cm -8*C(indexrd(1)- 2) + C(indexrd(1)- 3);
rd = -Diff*(d)/12/h/Cm*dt + rd;
end
if (mod(rd,h)~=0)
rd = rd + (h - mod(rd,h));
end
%eq 4
indexrd_new = find(grid1D>=rd);
%expl
for i = index(1):N-1
C(i) = C0(i) + dt*(Diff/h^2*(C0(i+1) + C0(i-1) - 2*C0(i)) + 2*Diff/i/h*(C0(i+1)-C0(i-1))/h );
end
for i = indexrd_new(1) :length(C)
C(i) = Cm;
end
for i = 1:index
C(i) = Cw;
end
if (mod(t,100)==0)   % => Every 50 time steps
for i=1:N-1
Cplot(i) = 0.5*(C(i+1) + C(i));
end
plot(grid1D,Cplot);
hold on
plot(grid1D(indexrd_new(1))*ones(2,1),[0.028 0.052])
%scatter(t*dt,Cw)
hold off
ylim([2.8e-2 5.2e-2])
xlim([0 L/2])
xlabel("length [m]")
ylabel("Concentration [mol/m3]")
message = sprintf('t=%d', ceil(t*dt));
time = annotation('textbox',[0.15 0.8 0.1 0.1],'String',message,'EdgeColor','k','BackgroundColor','w');
drawnow;
end
end```
.rtcContent { padding: 30px; }.lineNode {font-size: 12pt; font-family: "JetBrainsMono Nerd Font", Menlo, Monaco, Consolas, "Courier New", monospace; font-style: normal; font-weight: normal; }
Attached Images
 model.jpg (145.5 KB, 19 views)

 June 22, 2023, 05:32 #2 Senior Member   Matthew Join Date: Mar 2022 Location: United Kingdom Posts: 177 Rep Power: 4 The way to include Neumann boundary conditions into the diffusion is to use ghost cells at each point. You can then correctly fit in the ghost points that you get from the basic finite difference, that's how you do this normally. You have a moving boundary condition, so what you could do is fix the boundary with the transformation: r is mapped to r/R(t), then this fixes the boundary. The resulting equations you get will be a coupled set of equations for the nodes and the moving boundary. You can then use Newton-Raphson to solve these equations if they're nonlinear, or a matrix method if they're linear.

 Thread Tools Search this Thread Search this Thread: Advanced Search Display Modes Linear Mode

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are Off Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post charan3007 SU2 0 October 21, 2021 08:27 CFXer CFX 4 July 16, 2020 23:44 mpeppels CFX 11 August 22, 2019 07:30 Abo Anas CFX 27 May 11, 2018 01:44 amodpanthee CFX 31 April 19, 2018 18:02

All times are GMT -4. The time now is 05:10.

 Contact Us - CFD Online - Privacy Statement - Top