Non-uniform heat flux on a cube face defined by a matrix

I am trying to model the heating effects of a patterned beam (e.g. a multimode laser spot) on a large rectangular slab of material.

I can define a uniform heat flux and emissivity on the front face but I would like to be able to apply the heat flux as a 2D array of flux values (basically projecting the radiometric image). I know that this non-uniform heat flux needs to be defined as a @function handle but I just can't get it to work (because of "incompatible array sizes"). Do I need to make the heat flux map align with the mesh nodes?

% Create the PDE model
model = createpde("thermal","transient");

% Define a single cube using multi cube command
geo = multicuboid(1,1,1);

% Include the geometry in the model
model.Geometry = geo;

% Plot the geometry with face labels
pdegplot(model,"faceLabels","on")

% Generate mesh
generateMesh(model);

% Display the mesh (using default triangle size)
pdemesh(model);

% Set the thermal properties
thermalProperties(model,"ThermalConductivity",0.435, "MassDensity",1.450, "SpecificHeat",1109);

% Define this constant
model.StefanBoltzmannConstant = 5.670373E-8;

% CREATE THIS FUNCTION FOR THE NEXT STEP (A simple pattern for testing)
%function Qflux = externalHeatFluxI(region,~)
%Qflux = zeros(100);
%Qflux(:, 45:55) = 30;
%Qflux = Qflux+10;
%end

% Apply spatially varying heat flux to front face
X = zeros(100);
X = X+10;
region.x = X;
region.y = X;
region.z = X;
thermalBC(model,"Face",6, "HeatFlux",@externalHeatFluxI, "Vectorized","on")

% Set initial temperature
thermalIC(model,16);

% Specify Solution Times
tlist = 0:1:1000;

% Calculate Solution
thermalresults = solve(model,tlist);

% Plot Temperature Distribution
pdeplot3D(thermalresults.Mesh(:,end),ColorMapData=thermalresults.Temperature(:,end));

Answer 1

[Divyam]

Hi @Lee,

To set the boundary conditions for the thermal model, the function handle that is used to define the "HeatFlux" property needs to be defined. Since the function definition has not been attached with the code, here is a sample function that I defined to run the code without any issues.

function qflux = externalHeatFluxI(region, ~)
    % Initialize the heat flux array
    qflux = zeros(size(region.x));
    
    % Define your pattern (e.g., a simple Gaussian spot)
    % Adjust these parameters to match your desired pattern
    centerX = 0.5;
    centerY = 0.5;
    sigma = 0.1; % Standard deviation for Gaussian

    % Calculate the heat flux at each point
    for i = 1:length(region.x)
        % Calculate the distance from the center of the pattern
        distance = sqrt((region.x(i) - centerX)^2 + (region.y(i) - centerY)^2);
        
        % Apply a Gaussian profile
        qflux(i) = 30 * exp(-(distance^2) / (2 * sigma^2)) + 10;
    end
end

To test the entire code, directly add this function definition to file with the rest of the code and make the necessary changes to incorporate specific requirements. The attached file, "ThermalBoundaryConditions.m" can also be used to run the code on your end.