3D Multi-Domain PDE Modelling

The associated geometry for my problem is given below, where cell 1 is the current containing region and cell 2 is the surrounding space. The mesh is produced via gmsh and imported to MATLAB to deal with the limitations noted here:

An example code to do this is given below, please excuse the terrible use of global variables in this piece of demo code:

clear,close all; clc;

N = 3;  %Number of Equations
model = createpde(N);   %Value = #Variables 

%% Import gmsh mesh:
Total_Smoothed_Hollow_TwoCell
elements = msh.TETS; 
geom2 = geometryFromMesh(model,nodes,elements);

%% Variables:
global I_mag N_wire L_inner L_outer Z_height mu_o
I_mag = 20;                              
N_wire = 861;                            
Z_height = 8.89e-2;                      
L_inner = 10.42e-2;                      
L_outer = 15.24e-2;                      
mu_o = 4*pi*10^(-7);

%% Region 1:
c_coeff = zeros(9*N^2,1);
c_coeff = c_coeffassign_2v(c_coeff,N);       
specifyCoefficients(model,'m',0,'d',0,'c',c_coeff,'a',0,'f',@f_coeffunction,'Cell',1); %Outer

%%         UNCOMMENT ME FOR NON-SINGULAR 'RESULTS'        %%
%applyBoundaryCondition(model,'dirichlet','Face',[1,2,3,4,5,6,7,8,9,10],'u',[0;0;0]);

%% Region 2: 

c_coeff = zeros(9*N^2,1);
c_coeff = c_coeffassign_2v(c_coeff,N);      
specifyCoefficients(model,'m',0,'d',0,'c',c_coeff,'a',0,'f',[0;0;0],'Cell',2);
applyBoundaryCondition(model,'dirichlet','Face',[11,12,13,14,15,16],'u',[0;0;0]);

%% Solving: 

results = solvepde(model); 
vals = results.NodalSolution;

The code as given provides error message:

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  7.084396e-18. 
> In pde.EquationModel/solveStationary (line 49)
  In pde.PDEModel/solvepde (line 74)
  In MATHWORKS_PDETOOL_SIMPIF (line 54)

Other Attempts and Related Work:

The addition of a second boundary condition on this (the commented out one in the above) removes this singular warning however the solution is then nearly all zeros.

I have also tried using neumann conditions in various formats for the two regions to no avail. Combinations of dirichlet for one region and Neumann for the other have also been considered.

In the context of the PDE Toolbox, what does the singular warning actually mean? I'm assuming there is some error in my setup somewhere.

Answer 1

[Ravi Kumar]

Thanks for taking time to explain the problem with such a detail and for providing the reproduction code. I tried to understand and reproduce the behavior. However, one function you have used in f_coeffunction_2v.m, quad_find(), is missing from the attachment. It would help if you can attach that.

While I still don’t fully understand the setup of f-coefficient, I don’t have access to the paper, I think I sort have a clue what is going on. First, we need to understand how the two volumes (subdomains C1 and C2) that form the problem domain are related. As per the current setup, I don’t see any coupling between C1 and C2. That is C1 will have its own set of elements and C2 will have its own. Refer to this mesh plot below that I obtained using the commands:

elementsSection = model.Mesh.findElements('box',[0.0,0.5],[-0.6,0.6],[-0.6,0.6]);
figure
pdemesh(model.Mesh.Nodes,model.Mesh.Elements(:,elementsSection))
figure
pdemesh(model.Mesh.Nodes,model.Mesh.Elements(:,elementsSection),'FaceAlpha',0.4)

I guess this is not what you intend to model, right? I am guessing you want model the conductor at the center and volume of air around it to solve for A in that air domain. If so, then I think you need subdomains that don’t overlap but share interface boundaries. Would you have able to generate the mesh in gmsh such that two volumes stay separate?

Regards,

Ravi

[ADSW121365]

I think you're correct in this. Recreating the mesh as you mentioned is something that is possible (and I thought I was doing) so I will do some research into what I'm doing wrong and post an updated mesh when possible.

The coupling between the two regions is definitely an important step I'm missing out and I'm not really sure what coefficients/mathematics would be necessary to introduce such coupling.

[ADSW121365]

An updated conformal mesh has now been attached.

[Ravi Kumar]

Thanks. I think they are still seperate cells.

I did a thermal simulation with a constant temperature on one of the square's faces. In the steadystate, the whole domain needs to be at the assigned temperature. This is not happening as you can see from the code and figure below:

model = createpde('thermal','steadystate');   %Value = #Variables 
%% Import gmsh mesh:
Total_Smoothed_Hollow_TwoCell
nodes = msh.POS'; nodes = nodes./1000; %Unit Conversion from STEP MESH in mm
elements = msh.TETS; elements = elements(:,1:4); elements =elements';
geom2 = geometryFromMesh(model,nodes,elements);

thermalProperties(model,'ThermalConductivity',80);
thermalBC(model,'Face',12,'Temperature',100)
%internalHeatSource(model,1000,'Cell',1);
% thermalIC(model,0,'Cell',1)
% thermalIC(model,10,'Cell',2)
R = solve(model);

xq = linspace(-0.5E-3, 0.5E-3,100);
yq = linspace(-0.5E-3, 0.5E-3,100);
[XQ,YQ] = meshgrid(xq,yq);
ZQ = zeros(size(XQ));
iT = interpolateTemperature(R,XQ,YQ,ZQ);
contourf(XQ,YQ,reshape(iT,size(XQ))); colorbar;

Regards,

Ravi

[ADSW121365]

I have now solved this issue, you were correct and it was a mesh issue. Thanks!

[JJC]

How was the mesh issue solved?

[ADSW121365]

The Coherence; command at the start of my GMSH script.