MATLAB: Solving PDE of electrothermal coupling with f coefficient in stationary

couplefcoefficientjoule heatingpdesolvepdestationary

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.

Best Answer

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
Related Question