Solving PDE of electrothermal coupling with f coefficient in stationary

I am trying to solve a electrothermal system(joule heating) in stationary by using solvepde. The PDEs are shown below.

Equation(1) is for the E-field in stationary, where phi is voltage potential and sigma is the electrical conductivity of the material. Equation(2) is the relation between temperature and voltage potential in stationary, where T is temperature and k is thermal conductivity.

 %%E field for equation (1)
Emodel = createpde;
geometryFromEdges(Emodel,g);
Emesh = generateMesh(Emodel,'Hmax',0.8); % generate mesh 
% dirichlet boundary condition for E-domain
applyBoundaryCondition(Emodel,'dirichlet','Edge',8,'u',20);
applyBoundaryCondition(Emodel,'dirichlet','Edge',7,'u',0);
applyBoundaryCondition(Emodel,'neumann','Edge',1:6,'q',0,'g',0);
% laplace equation
specifyCoefficients(Emodel,'m',0,...
                          'd',0,...
                          'c',1,...
                          'a',0,...
                          'f',0);
Eresults = solvepde(Emodel);

ElectroThermal coupling

ETmodel = createpde;
geometryFromEdges(ETmodel,g);
generateMesh(ETmodel,'Hmax',0.8);  % generate mesh 
applyBoundaryCondition(ETmodel,'dirichlet','Edge',8,'u',273.15);
applyBoundaryCondition(ETmodel,'dirichlet','Edge',7,'u',273.15);
applyBoundaryCondition(ETmodel,'neumann','Edge',1:6,'q',1/kappa,'g',h*273.15/kappa);

 f = @(region,state) sigma.*(state.ux.^2 + state.uy.^2); %  <--------------u is for T here, 
                                                             but what should I do for letting this u represent
                                                            the potential voltage according to result of Emodel above? 
 %  region.ux means partial u/ partial x in this case. 
 specifyCoefficients(ETmodel,'m',0,...
                          'd',0,...
                          'c',kappa,...
                          'a',0,...
                          'f',f);

 ETresults = solvepde(ETmodel);

I want to let f be the right part of equation 2 (let u represent "phi" instead of T). or is there some essential way to represent the coupled PDE to solve?

Thank you in advance, Xing An.

Answer 1

[Alan Weiss]

If I understand your question, you would like f to represent the interpolated squared gradient of the previous solution to Emodel. So write a function to use the interpolated gradients of the solution Eresults, such as

function fout = myf(region,state,Eresults,sigma)
[gradx,grady] = evaluateGradient(Eresults,region.x,region.y)
fout = sigma.*(gradx.^2 + grady.^2);

Set this function as your f coefficient after getting Eresults and sigma into your workspace:

f = @(region,state)myf(region,state,Eresults,sigma)

Alan Weiss

MATLAB mathematical toolbox documentation

[Xing An]

Hi Alan, Thanks for reply and your answer. The code should work, but it gives me some errors as shown below.

Error using formGlobalKF2D Coefficient evaluation function, "@(region,state)fcoefETcoupled(region,state,Eresults,sigma)", was requested to calculate coefficients at 117 locations so should have returned a matrix with 117 columns. Instead it returned a matrix with 1 columns.

Error in pde.DiscretizedPDEModel/getStationaryFEMatrices (line 80) [K, F] = formGlobalKF2D(self.emptyPDEModel, self.p, self.t, self.coefstruct,u,time);

Error in pde.DiscretizedPDEModel/checkSpatialCoefsForUDependence (line 45) [K0,A0,F0,Q0,G0,H0,R0] = self.getStationaryFEMatrices(u0,tdummy);

Error in pde.DiscretizedPDEModel (line 34) obj = obj.checkSpatialCoefsForUDependence(u0,tdummy);

Error in pde.PDEModel/solvepde (line 72) femodel = pde.DiscretizedPDEModel(self,p,e,t,coefstruct,u0);

Error in geo (line 62) ETresults = solvepde(ETmodel);

Would you like to help me further?

[Alan Weiss]

Sorry, I wasn't able to run the code to debug it. I suspect that if you take a transpose somewhere you will get a reasonable result.

Turn on the debugger and set a break point in the code. See whether region.x is a row vector or column vector. It seems that the coefficient function wants a row vector as output, but perhaps region.x is a column vector. If so, take a transpose of region.x and region.y before calculating the value of fout, or take the transpose of fout before returning it.

Alan Weiss

MATLAB mathematical toolbox documentation

[Xing An]

Alan, thank you very much. The code works well.

[ADSW121365]

Alan, is it possible given identical meshes for the two problems to avoid the interpolation step in this approach by using results.XGradients etc?

I have a question related to this here if you have any thoughts: https://uk.mathworks.com/matlabcentral/answers/547154-using-stationary-results-as-a-source-function-pde-toolbox