Creating and solving a general 2D axisymmetric linear single scalar PDE with Dirichlet & Neumann BCs

See attached.

Seeking Psi(r,z,t) & first derivative of Psi_z.

How would one set up this general PDE in MATLAB - I see options for structural, thermal & modal, but not clear how to set it up for a general PDE such as this.

[Joseph Thodiyil]

Typo: Should be three Neumann type BCs and one Dirichlet type BC on page 1 of attachment.

Answer 1

[Joseph Thodiyil]

% Diffusion constant
D = 70.0; % m^2/month

% Vapor concentration values
Ci = 0.03; % gm/m^3 (Initial condition)
Cs = 0.05; % gm/m^3 (@ sink, station 2)

neqs = 1;
model = createpde(neqs); % 1 equation
height = 10.0; % meter
ra = 0.90;   z1 = 0;
rb = 0.95;   z2 = height;

Rect = [4;ra;rb;ra;rb;z1;z1;z2;z2]; % Axisymmetric model
g = decsg(Rect);
geom = geometryFromEdges(model,g);  % Creating the geometry ... not generating it correctly
pdegplot(model,'EdgeLabels','on')

% Apply Boundary Conditions

DBC = Cs; % gm/m^3
% Apply Dirichlet BC to edge 4 (z = L face)
applyBoundaryCondition(model,'dirichlet','Edge',4,...
                        'r',DBC);

NBC = [0 0 0.40]; % gm/m^2/month
% Apply Neumann BC to edge 1 (r = a face)                    
applyBoundaryCondition(model,'neumann','Edge',1,...
                        'g',NBC(1));

% Apply Neumann BC to edge 2 (z = 0 face)
applyBoundaryCondition(model,'neumann','Edge',2,...
                        'g',NBC(2));

% Apply Neumann BC to edge 3 (r = b face)
applyBoundaryCondition(model,'neumann','Edge',3,...
                        'g',NBC(3));

specifyCoefficients(model, 'm',0,'d',1,'c',D,'a',0,'f', 0);

IC = Ci; % gm/m^3
setInitialConditions(model,IC);

p = generateMesh(model,'Hmax',0.01); % 1 cm mesh size

pdeplot(model);

tlist = 1:1:12; % Time from one month to 12 months
results = solvepde(model,tlist);

Psi = results.NodalSolution;

pdeplot(model,'XYData',Psi(:,tlist(1)))
pdeplot(model,'XYData',Psi(:,tlist(12)))

[Joseph Thodiyil]

The script runs fine, but the geometry is being generated incorrectly.

I suppose that for axisymmetric models, the MATLAB PDE toolbox assumes that the axis of rotation is the vertical axis passing through r = 0.

Could you offer some guidance on how to fix this geometry issue.

Thanks,

Joseph Thodiyil

[Joseph Thodiyil]

Questions:

1. How do we determine at what time steady state concentration has been reached in the problem domain ( a < r < b, 0 < z < L)?
2. Want to plot derivative of Psi w.r.t. z (flux in z-direction) in the problem domain. How do we do that?

Answer 2

[Joseph Thodiyil]

Hello fellow MATLAB users:

I have revised the setup (Rev1 - page 2 to line up with the PDE toolbox geometry edges E1 ... E4).

Rev1 MATLAB script seems to be working, and somewhere between Month 2 & 3, steady state concentration is achieved.

However, I am not certain that the Cylindrical coordinate system formulation for the PDE has been employed in the solution of this problem.

Could someone in this MATLAB community review and advise.

In addition, I would like to plot the Vapor Diffusion Flux (J) in the problem domain as well.

Thanks,

Joseph Thodiyil

Answer 3

[Suleyman Sahin]

How can I convert a method of characteristics calculation of a 2D Model to 2D axisymmetric model ?

[Torsten]

Since the equations to solve are different, you'll have to write the code anew.

Answer 4

[Filipe Soares]

Hello Joseph. Did you find a way to solve this?

I have exactly the same problem: trying to solve the Helmholtz equation in 2D-axisymmetric domain...

Answer 5

[Carlo Silano]

Hello,

You made errors in creating the geometry. I haven't checked the rest of the code, I'm only responding to the geometry creation and the issue you reported.

At the link: https://it.mathworks.com/help/pde/ug/decsg.html (read 'Input Arguments'), you can find information on how to build the geometry. I'm attaching a sample code where I plot your geometry (I used a circular domain to represent the external environment, but it can also give you an idea of how to correctly write 'ns' and 'sf' when representing multiple geometric figures). I recommend reading the documentation and trying to modify various aspects of my code to understand.

Many of the things I've written are not necessary when defining the geometry (such as parameterization, etc.), but I believe they can be useful when building complex geometries and you want to reuse the code.

I hope I have been helpful. Have a good day.

% Link MATLAB page: https://it.mathworks.com/help/pde/ug/decsg.html
% Read Input Arguments

close all
clear
clc

% ---------------------------PARAMETRIZATION-------------------------------

height = 10; % meter
ra = 0.90;   z1 = 0;
rb = 0.95;   z2 = height;

% (1) Name; (2) R [m] centroid/center of gravity;
% (3) Z[m] centroid; (4)deltaR [m]; (5) deltaZ [m];
A={
'R1' 0 0 (rb-ra) (z2-z1)
};

R=vertcat(A{:,2});
Z=vertcat(A{:,3});

deltaR=vertcat(A{:,4});
deltaZ=vertcat(A{:,5});

% Create a circular domain containing your geometry. It will represent the
% surrounding environment

% (1) circle's name (2) Xcenter, (3) Ycenter
% (4) radius
C={
    'C1' 0.0 0.0 25
    };

% Create vectors containing only certain columns of C (circle coordinates) 
Xc=vertcat(C{:,2});
Yc=vertcat(C{:,3});
Rc=vertcat(C{:,4});

% -----------------------GEOMETRY DEFINITION-------------------------------

% For a circle, the first row contains 1. The second and third rows contain 
% the x- and y-coordinates of the center. The fourth row contains the radius
% of the circle. The radius must be a positive value

% For a rectangle, the first row contains 3, and the second row contains 4.
% The next four rows contain the x-coordinates of the starting points of the
% edges, and the four rows after that contain the y-coordinates of the starting
% points of the edges.

% Geometry description
gd = [1,Xc(1),Yc(1),Rc(1),0,0,0,0,0,0;...
      3,4,R(1)-deltaR(1)/2,R(1)+deltaR(1)/2,R(1)+deltaR(1)/2,R(1)-deltaR(1)/2,...
          Z(1)-deltaZ(1)/2,Z(1)-deltaZ(1)/2,Z(1)+deltaZ(1)/2,Z(1)+deltaZ(1)/2;]';


% name-space matrix
ns = char('C1','R1')';

% SET FORMULA
sf = 'C1+R1';

% Create a geometry and remove boundary conditions on the faces
geom = decsg(gd,sf,ns);

% Set number of equations
neqs = 1;
% Create a "model container"
model = createpde(neqs);

% geometryFromEdges for 2-D
geometryFromEdges(model,geom);

% -------------------------------PLOT--------------------------------------

% Edge plot
figure
pdegplot(geom,EdgeLabels="on")
xlim([-2 2])
ylim([-6 6])
title('Edge')

% Vertex plot
figure
pdegplot(geom,VertexLabels="on")
xlim([-2 2])
ylim([-6 6])
title('Vertex')

% Domain plot
figure
pdegplot(geom,FaceLabels="on")
xlim([-2 2])
ylim([-6 6])
title('Domain')