CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
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

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   November 1, 2022, 09:03
Default 1D simulation diffusion term, polar coordinates, moving boundary condition
  #1
New Member
 
Andrea
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.



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
File Type: jpg model.jpg (145.5 KB, 19 views)
exnihilo is offline   Reply With Quote

Old   June 22, 2023, 05:32
Default
  #2
Senior Member
 
Matthew
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

Reply

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.