CFD Online Logo CFD Online URL
Home > Forums > General Forums > System Analysis

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

Register Blogs Members List Search Today's Posts Mark Forums Read

LinkBack Thread Tools Search this Thread Display Modes
Old   November 1, 2022, 09:03
Default 1D simulation diffusion term, polar coordinates, moving boundary condition
New Member
Join Date: Oct 2022
Posts: 3
Rep Power: 3
exnihilo is on a distinguished road
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.

clc; clear all; close all
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);
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));
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;
        Cw = sum(int6)/(4/3*pi*r0^3) + Cw0;
        Cw = Cw0;
    % eq 3
    C(index(1)) = Cw/K;
    % eq 5
    if t == 1
        rd = r0 + h;
        d = -Cm + 8*Cm -8*C(indexrd(1)- 2) + C(indexrd(1)- 3);
        rd = -Diff*(d)/12/h/Cm*dt + rd;
    if (mod(rd,h)~=0)
        rd = rd + (h - mod(rd,h));
    %eq 4
    indexrd_new = find(grid1D>=rd);
    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 );
    for i = indexrd_new(1) :length(C)
        C(i) = Cm;
    for i = 1:index
        C(i) = Cw;
    if (mod(t,100)==0)   % => Every 50 time steps
        for i=1:N-1
            Cplot(i) = 0.5*(C(i+1) + C(i));
        hold on
        plot(grid1D(indexrd_new(1))*ones(2,1),[0.028 0.052])
        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');
.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
File Type: jpg model.jpg (145.5 KB, 19 views)
exnihilo is offline   Reply With Quote

Old   June 22, 2023, 05:32
Senior Member
Join Date: Mar 2022
Location: United Kingdom
Posts: 177
Rep Power: 4
hunt_mat is on a distinguished road
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.
hunt_mat is offline   Reply With Quote


Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

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 Off
Trackbacks are Off
Pingbacks are On
Refbacks are On

Similar Threads
Thread Thread Starter Forum Replies Last Post
Wind tunnel flow simulation boundary condition issue charan3007 SU2 0 October 21, 2021 08:27
Table bounds warnings at: END OF TIME STEP CFXer CFX 4 July 16, 2020 23:44
Radiation in semi-transparent media with surface-to-surface model? mpeppels CFX 11 August 22, 2019 07:30
My radial inflow turbine Abo Anas CFX 27 May 11, 2018 01:44
CFD analaysis of Pelton turbine amodpanthee CFX 31 April 19, 2018 18:02

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