MATLAB Answers

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

7 views (last 30 days)
Ahmad
Ahmad on 20 Oct 2019
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);
[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

  0 Comments

Sign in to comment.

Accepted Answer

Ravi Kumar
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:
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

  1 Comment

Ahmad
Ahmad on 21 Oct 2019
Hi Ravi,
Thanks for your help, the whole function works now! I hadn't thought of trying to transpose the gradx and grady instead, great catch! Thanks again!
Best regards,
Ahmad

Sign in to comment.

More Answers (0)

Sign in to answer this question.