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:
D_e = 29;
D_i = 14;
r_e = D_e/2;
r_i = D_i/2;
r1 = [1 0 0 r_e];
r2 = [1 0 0 r_i];
gdm = [r1; r2]';
g = decsg(gdm, 'R1-R2', ['R1'; 'R2']');
diffmodel = createpde(2);
geometryFromEdges(diffmodel,g);
mesh = generateMesh(diffmodel,'Hmax',1);
Dx = 0.55; Dy = 0.45;
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;
coefRGy = 0.35;
xx = (region.x);
yy = (region.y);
[gradx,grady] = evaluateGradient(resdiff,xx,yy,1,teval);
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];
teval = 100;
f_k = @(region,state) myf(region,state,resdiff,teval);
specifyCoefficients(structmodel,'m',0,'d',0,'a',0, ...
'c',c_k,'f',f_k);
applyBoundaryCondition(structmodel,"dirichlet","Edge", ...
5:8,'u',[0;0]);
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;
coefRGy = 0.35;
xx = (region.x);
yy = (region.y);
[gradx,grady] = evaluateGradient(resdiff,xx,yy,1,teval);
fx = [gradx.*coefRGx;grady.*coefRGy];
fout = (fx).';
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
