If you have access to MATLAB you can use the FEATool Multiphysics GUI to easily accomplish this; as it comes with built-in support for 4-series NACA wing profiles, automatic mesh generation, and FEniCS solver integration (you can use FEATool to export FEniCS py scripts and call the solver, or just export the mesh as a Dolfin/FEniCS xml grid data if you prefer).
You could possibly use the MATLAB FEA Toolbox and model the whole cross section as a full 3D linear elastic problem. Otherwise a 1D Euler-Bernoulli beam simplification might be suitable.
Possibly it could be due to the 'k' coeffcient being zero (then the heat equation isn't precisely defined).
If the PDE Toolbox still doesn't work it should be possible to simulate it with the Matlab FEA Toolbox. The following code should roughly be what you need:
fea.geom.objects = { gobj_circle([0 0],0.3,'R1') gobj_circle([0 0],0.35,'R2'), gobj_circle([0 0],1,'R3') }; fea.grid = gridgen( fea, 'hmax', 0.05 );
subplot(2,2,1) plotgeom(fea) subplot(2,2,2) plotgrid(fea) subplot(2,2,3) plotsubd(fea) subplot(2,2,4) plotbdr(fea)
fea.sdim = {'x' 'y'}; fea = addphys( fea, @heattransfer );
% Parameters for the human body rho_body = 985; C_body = 3470; k_body = 0.5; Q_body = 120; h_body = 4.5;
% Parameters for the insulation rho_insul = 1; C_insul = 1000; k_insul = 0+sqrt(eps); Q_insul = 0; h_insul = 26.3;
% Parameters for the air rho_air = 1; C_air = 1000; k_air = 0+sqrt(eps); Q_air = 0; h_air = 50;
fea.phys.ht.eqn.coef{1,end} = {rho_air rho_insul rho_body}; fea.phys.ht.eqn.coef{2,end} = {C_air C_insul C_body}; fea.phys.ht.eqn.coef{3,end} = {k_air k_insul k_body}; fea.phys.ht.eqn.coef{6,end} = {Q_air Q_insul Q_body};
T_amb = 293; fea.phys.ht.bdr.sel = [1 1 1 1]; fea.phys.ht.bdr.coef{1,end} = {T_amb T_amb T_amb T_amb};
fea = parsephys(fea); fea = parseprob(fea); fea.sol.u = solvestat(fea);
figure postplot( fea, 'surfexpr', 'T' )
I would recommend the classic DFG cylinder benchmarks where the drag and lift coefficients have been computed up to machine precision. Model set up https://www.featool.com/model-showcase/04_fluid_dynamics_03_flow_around_cylinder1 and references:
[1] John V, Matthies G. Higher-order finite element discretizations in a benchmark problem for incompressible flows. International Journal for Numerical Methods in Fluids 2001.
[2] Nabh G. On higher order methods for the stationary incompressible Navier-Stokes equations. PhD Thesis, Universitaet Heidelberg, 1998. Flow Around a Cylinder
Depending if your goal/task is to do it yourself, or if you can use toolboxes then both the PDE Toolbox and FEATool Multiphysics supports heat transfer simulations in Matlab, see for example the following heat exchanger tutorial:
https://www.featool.com/model-showcase/01_quickstart_02_heat_exchanger1
If you are allowed to use a toolbox and not have to code it yourself, here is a tutorial for porous media simulation with the Brinkman equations that might help, and Darcy's law is also possible.
Although it should work with PDE toolbox, you might have found a bug in which case you can maybe report it to the Mathworks and directly get help from them.
Alternatively, heat transfer with different material properties in different subdomains can also be simulated with the FEATool Multiphysics toolbox, for example as done in this example https://www.featool.com/model-showcase/02_heat_transfer_03_heat_transfer3
This seems to be almost exactly what you want. If you have your fem element data structure you can import it and make something similar.
MATLAB postprocessing and visualization of an axi-symmetric disc brake
You could also use the functions available with FEATool Multiphysics Matlab FEM toolbox to create unstructured plots. For example if have grid points and cell connectivities, you can import the grid as described here
https://www.featool.com/doc/grid.html#grid_imp
that is
fea.grid.p = % your imported grid points size[n_sdim,n_p]. fea.grid.c = % your imported cell connectivities size[n_nodes_per_cell,n_c]. fea.grid.a = gridadj(fea.grid.c,2); fea.grid.b = grid.bdr(fea.grid.p,fea.grid.c,fea.grid.a); fea.grid.s = ones(1,size(fea.grid.c,2));
Then if your solution corresponds to the grid points / nodes you can import it into the fea.sol.u solution field
fea.sdim = {'x' 'y'}; fea.dvar = {'u'}; fea.sfun = {'sflag1'}; fea.sol.u = % your solution data as a column vector. fea = parseprob(fea);
and then use either the Gui or the Matlab functions such as
postplot(fea,'surfexpr','u','isoexpr','u')
to plot surface plot and iso contour lines. Other plot options and Plotly output is also available as described in https://www.featool.com/doc/postplot_8m.html