Constant specification and normalization in PDEtool

I'm interested in using the PDE Tool in MATLAB to get concentration gradients through complex shapes, and thought I'd learn it by trying to reproduce a simple known result. For oxygen diffusion through a spheroid, the analytical expression is known and can be explicitly evaluated . For example, we can calculate the oxygen concentration through a spheroid of radius 500um with oxygen consumption rate a = 4.8 mmHg /s and a diffusion constant of D = 2 X 10E-9 m^2 / s, the oxygen profile through the spheroid with a constant outer boundary oxygen pressure of 100 mmHg.

Implementing this in MATLAB using PDEtool should be straight-forward as the problem has radial symmetry. Using PDETool, I draw a circle centred at (0,0) with a radius of 1. I then set the Dirchlet conditions on the circle to be 100 mmHg, and specify the PDE to be lapacian of u = a/D. So in the elliptic equation box, I set c = 1, a = 0. Now, here's my issue - the constant f required by PDETool should be proportional to a/D,but normalized so the domain -1 to 1 corresponds to -500 to + 500um. in SI units, a/D is 2.4 x 10E9 mmHg/ m^2. This is the same as 0.0024 mmHg / square micron. If the domain between 0 and 1 is 500um, 1um^2 = (1/500)^2 units, and thus we get f = -600 mmHg/ unit^2.

However, if I put this value in and solve, I get an incorrect result. If instead I change the value of f to -400 mmHg/unit^2 however, I get a result which closely matches the analytical solution, but I cannot see where this factor of 2/3 or so would come from. These are all shown below for clarity.

Does anyone know where I'm going wrong? I suspect there's two options - either I am doing something daft in my attempts at normalization OR the fact that PDEtool solves in 2D rather than 3D might be the root of my error, but as far as I am aware PDEtool is more than capable of dealing with such symmetrical problems. Any advice welcomed!

Answer 1

[Jim Joy]

Hi David,

The PDE Modeler App only allows you to model problems with 2D Geometries. Even though the Laplacian has similar forms in 2D and 3D, there are slight differences, and it is necessary to take the dimension of the problem into account.

If you are working in MATLAB R2017a or newer, you can solve the stated problem for a spherical region as shown in the code snippet below:

gm = multisphere(1);
model = createpde;
model.Geometry = gm;
applyBoundaryCondition(model,'dirichlet','Face',1,...
                     'u',100);
specifyCoefficients(model,'m',0,...
                        'd',0,...
                        'c',1,...
                        'a',0,...
                        'f',-600);
generateMesh(model);
sol = solvepde(model);

To plot the solution as a function of 'r', you can interpolate the solution along the x-axis using the "interpolateSolution" function, as shown below:

xq = (-1:0.01:1);
yq = zeros(1,length(xq));
zq = zeros(1,length(xq));
solInterp = interpolateSolution(sol,xq,yq,zq);
plot(xq,solInterp)

The output plot appears similar to the expected theoretical result you posted.

Best Regards,

Jim