MATLAB: Solving PDE for Hygromechanical Coupling with f coefficient as a function of previous PDE results

Question

Hi everyone. I'm trying to solve a hygromechanical coupling using solvepde. The problem consists of a transient diffusion equation and a stationary 2D elasticity. The PDEs are:

for the diffusion where D is the coefficient of diffusion and C is the diffusing concentration and:

for the 2D elasticity where c is the elasticity tensor and alpha is the swelling-shrinkage coefficient. I solved the diffusion equation first by:

%% Geometry (a hollow cylinder)
D_e = 29; % in mm

D_i = 14; % in mm
r_e = D_e/2; % exterior radius
r_i = D_i/2; % interior radius
r1 = [1 0 0 r_e];
r2 = [1 0 0 r_i];
gdm = [r1; r2]';
g = decsg(gdm, 'R1-R2', ['R1'; 'R2']');
%% Solving diffusion
diffmodel = createpde(2);
geometryFromEdges(diffmodel,g);
mesh = generateMesh(diffmodel,'Hmax',1);
Dx = 0.55; Dy = 0.45; % Coefficient of diffusion for x and y direction
c_xy = [Dx; Dy];
specifyCoefficients(diffmodel,"a",0, ...
    "c",c_xy, ...
    "d",1, ...
    "m",0, ...
    "f",[0;0]);
applyBoundaryCondition(diffmodel,"dirichlet","Edge", ...
    1:diffmodel.Geometry.NumEdges, ...
    'u',[0.07;0.07]);
setInitialConditions(diffmodel,0.05);
tlist = 0:0.1:10;
resdiff = solvepde(diffmodel,tlist);

And then, I want f to represent the interpolated gradient of the diffusion solution (resdiff) times alpha (coefficient of swelling-shrinkage), I wrote a function in new file, myf.m, which contains:

function fout = myf(region,state,resdiff,teval)
    coefRGx = 0.25; % swelling-shrinkage coef for x direction

    coefRGy = 0.35; % swelling-shrinkage coef for y direction

    xx = (region.x);
    yy = (region.y);
    [gradx,grady] = evaluateGradient(resdiff,xx,yy,1,teval); %teval because the previous solution is transient

    fx = [gradx.*coefRGx;grady.*coefRGy];
    fout = (fx);
end

then I describe the elasticity model by:

structmodel = createpde(2);
geometryFromEdges(structmodel,g);
meshstruct = generateMesh(structmodel,'Hmax',1);
E_p = 10e6;
mu_p = 0.40;
G = E_p/(2*(1+mu_p));
mu = 2*G*mu_p/(1-mu_p);
c_k = [2*G+mu; 0; G;   0; G; mu; 0;  G; 0; 2*G+mu];    % elasticity tensor for 2D PDE
teval = 100;                                           % evaluated time step
f_k = @(region,state) myf(region,state,resdiff,teval); % f-coefficient as a function
specifyCoefficients(structmodel,'m',0,'d',0,'a',0, ...
    'c',c_k,'f',f_k);
applyBoundaryCondition(structmodel,"dirichlet","Edge", ...
    5:8,'u',[0;0]);                                    % the center part is held in place
setInitialConditions(structmodel,[0;0]);               

Everything is alright until I tried to solve it:

structres = solvepde(structmodel);

where it returns me an error:

Error using formGlobalKF2D

Coefficient evaluation function,

"@(region,state)myf(region,state,resdiff,teval)", was

requested to calculate coefficients at 1500 locations so

should have returned a matrix with 1500 columns. Instead it

returned a matrix with 1 columns.

I was a bit surprised, but I thought it was just a problem of vector size, so I transposed the fout in my myf.m file, essentially editing it into:

function fout = myf(region,state,resdiff,teval)
    coefRGx = 0.25; % swelling-shrinkage coef for x direction
    coefRGy = 0.35; % swelling-shrinkage coef for y direction
    xx = (region.x);
    yy = (region.y);
    [gradx,grady] = evaluateGradient(resdiff,xx,yy,1,teval); %teval because the previous solution is transient
    fx = [gradx.*coefRGx;grady.*coefRGy];
    fout = (fx).'; % transposed
end

but it gave me another error:

Error using formGlobalKF2D

Coefficient evaluation function,

"@(region,state)myf(region,state,resdiff,teval)", was

requested to calculate coefficients at 1 locations so should

have returned a matrix with 1 columns. Instead it returned a

matrix with 2 columns.

So I was lost here. Could you tell me where I did wrong?

Thank you in advance,

Ahmad

Answer 1

You should be all set if you replace the last two lines in your function with the following one line:

fx = [gradx.'.*coefRGx;grady.'.*coefRGy];

Note, as error message says, size of fx should be a matrix of size 2x1500, instead your code was providing 3000x1 vector.

Regards,

Ravi