Fractional-Order Thermomechanical Modeling of Skin Tissue with Clinically Relevant Boundary Conditions: Convective–Radiative–Evaporative Cooling and Subcutaneous Elastic Foundation

Thermal therapies require precise prediction of the temperature and stress distributions in skin tissue to ensure efficacy while minimizing tissue damage. Classical bioheat models rely on oversimplified boundary assumptions, such as thermally insulated surfaces and mechanically free membranes, that fail to capture the physiological environment. This study develops a fractional-order dual-phase-lag bioheat model that incorporates clinically realistic conditions, namely simultaneous convective, radiative, and evaporative heat losses, active epidermal cooling, and subcutaneous mechanical restraint modeled through a Winkler elastic foundation. Both the Caputo and the Atangana–Baleanu (ABC) fractional derivatives are employed to represent memory effects in biological tissues. Analytical solutions in the Laplace–Fourier domain are obtained using displacement potential functions, with numerical inversion carried out via the Stehfest algorithm and Gaussian quadrature. The results show that realistic boundary conditions substantially alter the thermomechanical response: convective and evaporative cooling reduce surface temperatures and penetration depths, whereas active cooling permits deeper heating without epidermal damage. The Winkler foundation yields higher compressive stresses than traction-free models, and the ABC operator produces smoother responses than the Caputo operator. Overall, the model reveals the trade-offs between thermal efficacy and mechanical safety, thereby bridging bioheat modeling and clinical practice.