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

7 views (last 30 days)
Commented: Ahmad on 21 Oct 2019
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);
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);
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?

Ravi Kumar on 21 Oct 2019
You should be all set if you replace the last two lines in your function with the following one line:
Note, as error message says, size of fx should be a matrix of size 2x1500, instead your code was providing 3000x1 vector.
Regards,
Ravi

Hi Ravi,
Best regards,