k-Wave

1. cgenes
   

   Member
   Joined: Apr '19
   Posts: 6

   Hi

   I'm doing a simulation of thermal diffusion from a skin/air surface; kwave gives very odd results (I would copy in an image but the box won't let me)

   I'll describe: it has pronounced 'stripes' in it rather than smooth continuous diffusion as when doing a homogenous medium

   also, I can't set the boundaries to be insulating, so heat doesn't 'diffuse back round' - but I believe that is the topic of another post (I'm not quite sure how to implement this code though https://www.mathworks.com/matlabcentral/fileexchange/75071-matlab-discrete-trigonometric-transform-library)

   thanks

   Charlie

   Posted 4 years ago #

2. cgenes
   

   Member
   Joined: Apr '19
   Posts: 6

   just to add to that: the way I'm setting up my air/surface media is like this

   %% specify the medium parameters
   % create the computational grid
   Nx = 65; % number of grid points in the x (row) direction
   Ny = 65; % number of grid points in the y (column) direction
   % Nz = 201;

   % calculate grid point spacing
   dx = (20 * 1e-3) / 65 % this is the height (or width / depth) of the 'world' in mm * m / the resolution aka grid points
   dy = (20 * 1e-3) / 65

   dx = dx; % grid point spacing in the x direction [m]
   dy = dy; % grid point spacing in the y direction [m]
   % dz = 1e-4;

   kgrid = kWaveGrid(Nx, dx, Ny, dy);

   % define medium properties

   airStart = 40;

   medium.density = ones(Nx, Ny);
   medium.density(:, :) = 1079; % of tissue [kg/m^3]
   % medium.density(:, airStart:end) = 1.255; % of air [kg/m^3]
   medium.density(airStart:end, :) = 1.255; % of air [kg/m^3]
   figure, imagesc(medium.density)

   medium.thermal_conductivity = ones(Nx, Ny);
   medium.thermal_conductivity(:, :) = 0.52; % tissue [W/(m.K)]
   % medium.thermal_conductivity(:, airStart:end) = 26.02; % of air [W/(m.K)]
   medium.thermal_conductivity(airStart:end,:) = 26.02; % of air [W/(m.K)]

   medium.specific_heat = ones(Nx, Ny);
   medium.specific_heat(:, :) = 3540; % of tissue [J/(kg.K)]
   % medium.specific_heat(:, airStart:end) = 0718; % of air [J/(kg.K)]
   medium.specific_heat(airStart:end, :) = 718; % of air [J/(kg.K)]

   % define medium properties related to perfusion
   medium.blood_density = ones(Nx, Ny); % [kg/m^3]
   medium.blood_density(:, :) = 1060; % [kg/m^3]
   % medium.blood_density(:, airStart:end) = 0; % [kg/m^3]
   medium.blood_density(airStart:end, :) = 0; % of air [kg/m^3]

   medium.blood_specific_heat = ones(Nx, Ny); % [J/(kg.K)]
   medium.blood_specific_heat(:, :) = 3617; % [J/(kg.K)]
   % medium.blood_specific_heat(:, airStart:end) = 0; % [J/(kg.K)]
   medium.blood_specific_heat(airStart:end, :) = 0; % of air [J/(kg.K)]

   medium.blood_perfusion_rate = ones(Nx, Ny); % [1/s]
   medium.blood_perfusion_rate(:, :) = 0.01; % [1/s]
   % medium.blood_perfusion_rate(:, airStart:end) = 0; % [1/s]
   medium.blood_perfusion_rate(airStart:end, :) = 0; % of air [1/s]

   medium.blood_ambient_temperature = ones(Nx, Ny); % [degC]
   medium.blood_ambient_temperature(:, :) = 37; % [degC]
   % medium.blood_ambient_temperature(:, airStart:end) = 22; % [degC]
   medium.blood_ambient_temperature(airStart:end, :) = 22; % of air [degC]

   Posted 4 years ago #

3. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,643

   Hi Charlie,

   In the acoustic case, the large change in material properties between tissue and air can cause some problems (see for example this post). In the thermal case, I haven't actually tried simulating this example before. Do you think that the simulated temperature field looks non-physical? What happens if you decrease the time step, or if you artificially move the properties of air (in particular the density) closer to skin?

   Regarding the thermal boundary conditions, you first need to download this toolbox and add it to the MATLAB path (try running one of the included examples to make sure the compiled functions work). Then you can set medium.boundary_condition = 'insulating' or medium.boundary_condition = 'conducting'. Unfortunately, these flags are only compatible with homogeneous material properties. It's on the to-do list to extend it to the more general case.

   Hope that helps,

   Brad

   Posted 4 years ago #

4. cgenes
   

   Member
   Joined: Apr '19
   Posts: 6

   HI Brad

   thanks very much for your reply.

   I think I might know what's going on and why it looks so weird - although I'm probably misunderstanding something LOL

   because air density and conductivity are orders of magnitude different than in tissue, heat from the tissue escapes into the air layer which then quickly diffuses to the 'boundaries' - where I think it gets stuck. so then the air layer massively heats up unrealistically, and causes all sorts of weird effects. To contradict myself here though (LOL), shouldn't the heat start appearing on the 'other side' (in my case the tissue side) of the world?

   if my first explanation does make sense, I'm guessing if your plug-in worked with heterogeneous layers I could use that to make the world boundaries conducting? (is this analogous to the perfectly matched layer [PML] the user can turn on an off with the kspace function?)

   Or as a hack, I could make the world massive in comparison to the tissue block?

   thanks

   Charlie

   Posted 4 years ago #

5. cgenes
   

   Member
   Joined: Apr '19
   Posts: 6

   ON further thought - perhaps the air/surface interface is much more complicated than just changing the coefficients for thermal diffusion, as I would have though radiative or convection processes take over?

   the Matlab thermal toolbox has this 'boundaries' function
   https://www.mathworks.com/help/pde/ug/pde.thermalmodel.thermalbc.html
   which I think accounts for convection and radiative heat transfer

   is this a much more likely reason why kWave isn't great at modelling the air/surface boundary?

   Posted 4 years ago #

6. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,643

   HI Charlie,

   Yes - eventually kWaveDiffusion will allow the world boundaries to be conducting. I'm working on this at the moment. In the meantime, you can make the domain bigger as you say.

   The thermal boundary conditions in the MATLAB PDE toolbox are analogous to the ones in kWaveDiffusion, but can also be made more general. Setting 'Temperature', 0 in MATLAB is equivalent to a conducting BC in kWaveDiffusion, and setting 'HeatFlux', 0 is equivalent to an insulating BC. You can also specify more complicated BC (and combinations) in MATLAB which are not possible in kWaveDiffusion.

   The general PDE solved using createpde('thermal', 'transient') is the same as kWaveDiffusion (ignoring the perfusion which is modelled in k-Wave). In other words, the thermal physics is the same.

   Hope that helps,

   Brad.

   Posted 4 years ago #

RSS feed for this topic

Reply

You must log in to post.