# Introduction to Modeling Surface Reactions in COMSOL

High-intensity lasers incident upon a material that is partially transparent will deposit power into the material itself. If the absorption of the incident light can be described by the Beer-Lambert law, it is possible to model this power deposition using the core functionality of COMSOL Multiphysics. We will demonstrate how to model the absorption of the laser light and the resultant heating for a material with temperature-dependent absorptivity.When light is incident upon a semitransparent material, some of the energy will be absorbed by the material itself. If we can assume that the light is single wavelength, collimated (such as from a laser), and experiences very minimal refraction, reflection, or scattering within the material itself, then it is appropriate to model the light intensity via the Beer-Lambert law. This law can be written in differential form for the light intensity I as:where z is the coordinate along the beam direction and \alpha(T) is the temperature-dependent absorption coefficient of the material. Because this temperature can vary in space and time, we must also solve the governing partial differential equation for temperature distribution within the material:where the heat source term, Q, equals the absorbed light. These two equations present a bidirectionally coupled multiphysics problem that is well suited for modeling within the core architecture of COMSOL Multiphysics. Let’s find out how…We will consider the problem shown above, which depicts a solid cylinder of material (20 mm in diameter and 25 mm in length) with a laser incident on the top. To reduce the model size, we will exploit symmetry and consider only one quarter of the entire cylinder. We will also partition the domain up into two volumes. These volumes will represent the same material, but we will only solve the Beer-Lambert law on the inside domain — the only region that the beam is heating up.To implement the Beer-Lambert law, we will begin by adding the General Form PDE interface with the Dependent Variables and Units settings, as shown in the figure below. Settings for the implementing the Beer-Lambert law. Note the Units settings.Next, the equation itself is implemented via the General Form PDE interface, as illustrated in the following screenshot. Aside from the source term, f, all terms within the equation are set to zero; thus, the equation being solved is f=0. The source term is set to Iz-(50[1/m]*(1+(T-300[K])/40[K]))*I, where the partial derivative of light intensity with respect to the z-direction is Iz, and the absorption coefficient is (50[1/m]*(1+(T-300[K])/40[K])), which introduces a temperature dependency for illustrative purposes. This one line implements the Beer-Lambert law for a material with a temperature-dependent absorption coefficient, assuming that we will also solve for the temperature field, T, in our model. Implementation of the Beer-Lambert law with the General Form PDE interface.Since this equation is linear and stationary, the Initial Values do not affect the solution for the intensity variable. The Zero Flux boundary condition is the natural boundary condition and does not impose a constraint or loading term. It is appropriate on most faces, with the exception of the illuminated face. We will assume that the incident laser light intensity follows a Gaussian distribution with respect to distance from the origin. At the origin, and immediately above the material, the incident intensity is 3 W/mm2. Some of the laser light will be reflected at the dielectric interface, so the intensity of light at the surface of the material is reduced to 0.95 of the incident intensity. This condition is implemented with a Dirichlet Boundary Condition. At the face opposite to the incident face, the default Zero Flux boundary condition can be physically interpreted as meaning that any light reaching that boundary will leave the domain. The Dirichlet Boundary Condition sets the incident light intensity within the material.With these settings described above, the problem of temperature-dependent light absorption governed by the Beer-Lambert law has been implemented. It is, of course, also necessary to solve for the temperature variation in the material over time. We will consider an arbitrary material with a thermal conductivity of 2 W/m/K, a density of 2000 kg/m3, and a specific heat of 1000 J/kg/K that is initially at 300 K with a volumetric heat source.The heat source itself is simply the absorption coefficient times the intensity, or equivalently, the derivative of the intensity with respect to the propagation direction, which can be entered as shown below. The heat source term is the absorbed light.Most other boundaries can be left at the default Thermal Insulation, which will also be appropriate for implementing the symmetry of the temperature field. However, at the illuminated boundary, the temperature will rise significantly and radiative heat loss can occur. This can be modeled with the Diffuse Surface boundary condition, which takes the ambient temperature of the surroundings and the surface emissivity as inputs. Thermal radiation from the top face to the surroundings is modeled with the Diffuse Surface boundary condition.It is worth noting that using the Diffuse Surface boundary condition implies that the object radiates as a gray body. However, the gray body assumption would imply that this material is opaque. So how can we reconcile this with the fact that we are using the Beer-Lambert law, which is appropriate for semitransparent materials?We can resolve this apparent discrepancy by noting that the material absorptivity is highly wavelength-dependent. At the wavelength of incident laser light that we are considering in this example, the penetration depth is large. However, when the part heats up, it will re-radiate primarily in the long-infrared regime. At long-infrared wavelengths, we can assume that the penetration depth is very small, and thus the assumption that the material bulk is opaque for emitted radiation is valid.It is possible to solve this model either for the steady-state solution or for the transient response. The figure below shows the temperature and light intensity in the material over time, as well as the finite element mesh that is used. Although it is not necessary to use a swept mesh in the absorption direction, applying this feature provides a smooth solution for the light intensity with relatively fewer elements than a tetrahedral mesh. The plot of light intensity and temperature with respect to depth at the centerline illustrates the effect of the varying absorption coefficient due to the rise in temperature. Plot of the mesh (on the far left) and the light intensity and temperature at different times. Light intensity and temperature as a function of depth along the centerline over time.Here, we have highlighted how the General Form PDE interface, available in the core COMSOL Multiphysics package, can be used for implementing the Beer-Lambert law to model the heating of a semitransparent medium with temperature-dependent absorptivity. This approach is appropriate if the incident light is collimated and at a wavelength where the material is semitransparent.Although this approach has been presented in the context of laser heating, the incident light needs only to be collimated for this approach to be valid. The light does not need to be coherent nor single wavelength. A wide spectrum source can be broken down into a sum of several wavelength bands over which the material absorption coefficient is roughly constant, with each solved using a separate General Form PDE interface. In the approach presented here, the material itself is assumed to be completely opaque to ambient thermal radiation. It is, however, possible to model thermal re-radiation within the material using the Radiation in Participating Media physics interface available within the Heat Transfer Module.The Beer-Lambert law does assume that the incident laser light is perfectly collimated and propagates in a single direction. If you are instead modeling a focused laser beam with gradual variations in the intensity along the optical path then the Beam Envelopes interface in the Wave Optics Module is more appropriate.In future blog posts, we will introduce these as well as alternate approaches for modeling laser-material interactions. Stay tuned!Merci beaucoup. Could I use this for 2 layers of material ? or I need special boundaries condition at the interface. Thank you for your article.Dear Emma, Yes, you can use this approach for two layers of materials, simply add another “General Form PDE” interface and enter a different absorption coefficient in the different layers. Best Regards,Is there a way to extend this approach to model beer-lambert law in x, y and z directions simultaneously?Dear Natan, If you have multiple incident beams, then you can add multiple General Form interfaces for different beam directions. If you are more generally trying to model reflection and scattering, then we suggest that you contact your COMSOL Support Team for a personalized discussion about your modeling needs.Thank you for your response. It does not work adding multiple General Form interfaces since a single domain is modeled. However, I will try another strategies.Dear Dr. Walter, Just a remark: We got a solution to our question. Thank you! Best regards, NatanDear All, I am a beginner on COMSOL multiphysics, and concerning this tutorial I have some question , first of all I think in the Beer law is Iz=-alpha(T)* I which will give Iz+alpha(T)*I not minus; the second question is Iz is always negative so why Q=Iz and not Q=-Izfinnally my last question why you used two different domains inside where the beam heats up and outside. I have tried it using only one domaine (the bigger cylinder) and the results were terriblethank you for your fast reply Best Regards, Nabil BELGHACHEMDear Nabil, Yes, you are correct. One term that was assumed here that is that the beam is propagating in the negative z-direction. Regarding your model: You would want to contact your COMSOL Support Team about this.Dear Dr Walter, thank you very much for your fast reply. Now I understand, this make sens .Best Regards, Nabil BELGHACHEMDear all, I did not begin to use Comsol yet because I don’t know if this could be the right software to solve my problem. I would need to simulate the ultrashort laser pulses interaction with semicondutcor. To do it, I have two coupled equations for electron and lattice temperature evolution, a generalized Lambert-Beer law and an equation describing the evolution od electron density due to absorption processes. My problem is that the parameters in these equations depend on the electron and lattice temperature and on electron density (which are also dynamic variables). Do you suggest to use Comsol for this kind of problem and how? Thank you. Best regards.CaterinaDear Caterina, What you’re describing certainly sounds like an interesting and doable problem. You may want to search on the usage of the two-temperature model in COMSOL. If you’re looking to engage other COMSOL users, we would encourage you to post on the discussion forum: http://www.comsol.com/community/forums/ And if you’d like to contact a COMSOL representative about purchasing a license for this type of work: https://www.comsol.com/contactDear Caterina, Briefly speaking, yes, what you are describing does all sound like a reasonable (albeit possibly quite complex) problem to solve. We would of course suggest that you investigate using COMSOL for this type of problem, yes.Walter how did you make those plots side by side. I was able to replicate your model easily enough, but not the post processing. Can you point me in the right direction?Hello Peter, You can add a “Deformation” subfeature to any plot and specify that the plot be moved by a fixed quantity to pattern several different cases in one image. You may also find these resources helpful: https://www.comsol.com/offers/Postprocessing-and-Visualization-Handbook-Part-1 https://www.comsol.com/offers/Postprocessing-and-Visualization-Handbook-Part-2Hello,thanks for the tutorial, this was very helpful. I have implemented now the Beer-Lambert-law, but now I would like to move the laser beam along a certain distance with a certain speed. How can I implement that in the existing model?Thanks in advance, BastianHello Bastian,Yes, you can move the location of the laser around on the surface of the part. Within the Dirichelet boundary condition simply make the expression for incident intensity a function of position and time. For an example model that shows you how to set up such time and space-dependent expressions, please see: https://www.comsol.com/model/laser-heating-of-a-silicon-wafer-13835Best Regards,Hello Walter,thanks for your help, I was able to implement the moving laser. Now I have another question:I have two layers of material, the upper one has a low absorption coefficient, the lower one has a high absorption coefficient. The two parts have an ideal thermal contact. So basically a small amount (about 10-20%) of the laser energy is absorbed in the upper part (which has a thickness of 2 mm), the rest reaches the lower part and is absorbed in the first 100 microns of the material. Both material interactions are with beer-lambert law.Is it possible to do this with one heat source? At the moment I implemented simply two heat sources at the top of each layer and two dirichlet boundary conditions with the different absorptions coefficients. It would be nice to implement it in a way that if the heat reaches the lower material the beer-lambert law is calculated with the new absorption coefficient.Thanks in advance! Best regards, BastianHello Bastian, What you could do, in the situation that you’re describing, is assume that the incoming light in incoherent (probably fair to do if your first layer is 2mm, and if we assume an effective wavelength of around 1um) Under those assumptions, you could implement a bi-directional Beer-Lambert law approach, solving for both forward- and backward-propagating light. Dear all, by the way is the equation for the gaussian intensity distribution correct? Because in literature I found the following equation: I_Gauss = I_max*exp((-2*(x^2+y^2))/r_g^2), where I_max is the maximum intensity at r = 0, and r_g is the radius of the gaussian beam. I_max is then defined as I_max = (2*P)/(pi*r_g^2), where P is the power of the laser. Thansk in advanceHi Walter, I am a beginner with comsol multiphonics. I am trying to model sunlight absorption in side a 50 micro layer of material using beer Lambert’s law. Is it appropriate to use beer Lambert’s law in my case? My second question is about the heat source. Is it always Iz, the differential form? regardless of having the general form PDE or not? I have created a 2 D model , I don’t have the general form PDE, can the software still solve the differential form equation? Thanks SenHello, Thank you for the blog! The model appears listed in the Application Gallery but cannot be downloaded due to an html error which redirects to the current web page. Please, could you make it availbale there? Cheers!Hello Ciprian, This should now be addressed. Best Regards,Hello again Walter, and thank you for making the model available. It is not clear for me yet where does the 40[K] come from in the expression of the absroption coefficient. Would you clarify it for us, please? Many thanks! CiprianHello Ciprian, The 40[K] term is merely there to demonstrate a temperature dependence of the absorption coefficient. It is arbitrary, and meant just for demonstration purposes. I have a question. Downloaded the model from a given link, but when I look at heat source equation, I only can choose Q=Q0, but not the one that in on the screenshot. What don’t I understand? Thanks in advance, George AloyanHi George, Thanks for your comment! For questions related to your modeling, please contact our support team.Online support center: https://www.comsol.com/support Email: support@comsol.comBest, CatyThank you for response. And question about this model – why Iz is used, but not I? Thanks in advance, George AloyanHi George,Thank you for your additional comment. When you contact our support team, they can help address this question as well.Best, CatyHello, This model is provided in the application gallery but the version is 5.2a. Can you please provide it in older version?Hi Saurabh,Thanks for your comment.For your question, please contact our Support team.Online Support Center: https://www.comsol.com/support Email: support@comsol.comBrazilChinaDenmarkFinlandFranceGermanyIndiaItalyJapanNetherlandsNorwayPortugalRussiaSpainSwedenSwitzerlandUnited Kingdom

Say you are working on a modeling case where loads are moving in such a way that they cross over different mesh elements and boundaries during the simulation. In these cases, among other instances, you may want to apply a boundary condition to only part of the geometrical boundary or only under certain conditions. In this blog post, we’ll discuss how you can utilize the flexibility of COMSOL Multiphysics to handle such situations.In the mathematical treatment of partial differential equations, you will encounter boundary conditions of the Dirichlet, Neumann, and Robin types. With a Dirichlet condition, you prescribe the variable for which you are solving. A Neumann condition, meanwhile, is used to prescribe a flux, that is, a gradient of the dependent variable. A Robin condition is a mixture of the two previous boundary condition types, where a relation between the variable and its gradient is prescribed.The following table features some examples from various physics fields that show the corresponding physical interpretation.Within the context of the finite element method, these types of boundary conditions will have different influences on the structure of the problem that is being solved.The Neumann conditions are “loads” and appear in the right-hand side of the system of equations. In COMSOL Multiphysics, you can see them as weak contributions in the Equation View. As the Neumann conditions are purely additive contributions to the right-hand side, they can contain any function of variables: time, coordinates, or parameter values.Let’s consider a heat transfer problem where a circular heat source with a radius R is traveling in the x direction with a velocity v. Its intensity has a parabolic distribution with a peak value q_0. A mathematical description of this load could beFor a traveling load, it is obviously not possible to have domain boundaries, or even a mesh, that fits the load distribution at all times.The load distribution itself can be entered in a straightforward manner. Since the variable for the radial coordinate, r, will be used in two places, it is a good idea to define it as a variable. The entire input for the moving heat source is shown below. Parameters describing the moving heat source. The variable describing the local radial coordinate from the current center of the traveling heat source. Input of the heat flux.The results from a time-dependent simulation using such data are depicted in the following animation. Symmetry is assumed in the yz-plane, so the load is actually applied on a moving semicircle. Animation of the temperature distribution as the heat source travels along the bar.Where a Dirichlet condition is given, the dependent variable is prescribed, so there is no need to solve for it. Equations for such degrees of freedom can thus be eliminated from the problem. Dirichlet conditions therefore change the structure of the stiffness matrix. When looking in the Equation View of COMSOL Multiphysics, these conditions will appear as constraints.Assume that you want the traveling spot to prescribe the temperature as exactly 450 K. This may be a bit artificial, but it displays an important difference between the Neumann condition and the Dirichlet condition. If you were to add a Temperature node and enter a similar expression ( if(r < R,450[K],0)), it would mean setting the temperature to absolute zero on the part of the boundary that is not covered by the hot spot.The intention is, however, to switch off the Dirichlet condition outside of the hot spot. There’s a small trick for doing so. If you instead enter if(r < R,450[K],ht.Tvar) as the prescribed value, you will get the intended behavior (shown in the following animation). Settings for a conditional Dirichlet condition. Animation of the temperature distribution as the prescribed temperature spot travels along the bar. In order to understand how this works, enable the Equation View, and look at the implementation of the Dirichlet condition (in this case, a prescribed temperature): Enabling the Equation View. The Equation View for the Temperature node.The constraint is formulated as ht.T0-ht.Tvar, which implicitly means ht.T0-ht.Tvar = 0. The first term is the prescribed temperature, which you enter as input. The second term is just the temperature degree of freedom cast into a variable. This constrains the temperature to be equal to the given value, unless the given value happens to be the string ht.Tvar. In that case, the symbolic algebra during assembly will reduce the expression to ht.Tvar-ht.Tvar, and further to zero. And with the constraint expression being 0, there is no constraint.In COMSOL Multiphysics, there are actually two possible implementations of a Dirichlet condition. The default case is the pointwise constraint, as referenced above, but you can also use a weak constraint. In the weak constraint, equations are added rather than removed. The heat fluxes needed to enforce the prescribed values of the temperature are then added as extra degrees of freedom (Lagrange multipliers).You can use essentially the same trick to make a weak constraint conditional, with just a small twist weaved into the mix. Using weak constraints requires that you first enable the Advanced Physics Options. Enabling the Advanced Physics Options.When weak constraints have been selected in a node within a physics interface, there will be extra degrees of freedom for the Lagrange multiplier. In our case, the degree of freedom has the name T_lm.If the same expression for the temperature as that shown above is applied, the extra degree of freedom will not get any stiffness matrix contribution on the part of the boundary where the Dirichlet condition is switched off. The stiffness matrix will thus become singular. To avoid this situation, change if(r < R,450[K],ht.Tvar) to if(r < R,450[K],ht.Tvar-T_lm*1e-2). The multiplier used for T_lm differs between varying models and physics fields, and it may require some tuning.The reason this works, as a textbook might note, is “left as an exercise to the reader”. Settings for a conditional Dirichlet condition when using weak constraints.The Robin conditions generally contribute to both the stiffness matrix and the right-hand side. The structure of the stiffness matrix is not affected, but values are added to existing positions. The Robin conditions also appear as weak contributions in the Equation View. Turning these conditions into functions of time, space, and other variables is no different than doing so for Neumann conditions.It is, however, interesting that by selecting appropriate values, you can actually morph Robin conditions into acting as approximative Dirichlet or Neumann conditions. This is especially important for cases where you want to switch between the two boundary condition types during a simulation.To create a Dirichlet condition, you assign a high value of the “stiffness”, for instance, a spring constant or heat transfer coefficient. In mathematical terms, this is actually a penalty implementation of the Dirichlet condition. The higher the stiffness, the greater the accuracy of the prescribed value for the degree of freedom. But there is a caveat: A very high stiffness will harm the numerical conditioning of the stiffness matrix. For a heat transfer problem, a starting point for choosing a “high” heat transfer coefficient \alpha could bewhere \lambda is the thermal conductivity and h is a characteristic element size. The same idea can be applied to other physics fields by replacing \lambda with the appropriate material property (i.e., Young’s modulus in solid mechanics). The factor 1000 is just a suggestion and can be replaced by 104 or 105.If you were to use convection to model the moving 450 K spot from the previous example, you could utilize the settings shown below. The built-in variable h for the element size is applied in the expression. Using a convection condition to prescribe the temperature.In the good old days, when I first began using finite element ysis, it was sometimes not possible to prescribe nonzero displacements in finite element programs for structural mechanics. The limitation was imposed by the added complexity of the programming. If this was the case, the best option was to use the penalty method by adding a predeformed stiff spring. You wouldn’t want it to be too stiff though; in those days, single precision arithmetic was still in use!Let’s turn our attention toward approximating a Neumann condition. We want a heat flux that is independent of the surface temperature. In the case of heat transfer, the Robin condition states that the inward heat flux q iswhere \alpha is the heat transfer coefficient, T is the temperature at the boundary, and T_{\textrm{ext}} is the external temperature. So if T_{\textrm{ext}} is much larger than the expected temperature on the surface, then q \approx \alpha T_{\textrm{ext}}. The strategy is then to select an arbitrary, very large T_{\textrm{ext}} and compute a suitable heat transfer coefficient, as highlighted below. Using a convection condition to prescribe the heat flux.Designers actually use this principle to introduce a fixed force into real-life mechanical components: a prestressed long weak spring. If the predefomation of the spring is much larger than the displacement of the parts to which the spring is connected, the force in the spring will be almost constant.When a boundary condition is limited by a Boolean expression like if(r < R,..., then it is more than likely that the border of the region to which it is applied will not follow the edges of the mesh elements. This will introduce discretization errors.For a Neumann or Robin condition, a numerical integration is performed over each finite element. The value of the function is evaluated at a number of discrete Gauss points in the element. If the size of the mesh elements is large in comparison to the geometrical size of the load, then the exact number of Gauss points covered by the load can significantly affect the total load. As such, there should be several elements on the patch covered by the load at any instant. A small change in the location of the load may alter the number of contributing integration points. (In reality, the number of integration points is larger.)Dirichlet conditions, at least in the pointwise sense, are instead applied to the mesh nodes. In the figure below, the temperature distribution and heat flux are shown for a certain time when simulating the moving circular spot with a prescribed temperature of 450 K. In front of the hot spot, a darker shade with 260 K is visible. Since the initial and ambient temperatures in the simulation are 293 K, this is not expected. It is a numerical artifact that is related to the fact that not all nodes on each element have a Dirichlet condition. At a discontinuity in Dirichlet conditions, there will be singularities. This is a topic of discussion in a previous blog post. Refining the mesh will reduce such an effect.The green arrows in the following figure represent the nodes at which an influx of heat is created as a reaction to prescribing the temperature. With the mesh density in the model, the approximation of a semicircle will be rather rough. Temperature distribution and heat flux around the semicircular prescribed temperature.There are many ways in which the solution can enter your boundary conditions. This will generally introduce nonlinearities, which are automatically detected by COMSOL Multiphysics.As an example, let’s look at a beam featuring a support that is placed slightly below it, inhibiting further movement after a certain deflection. This can be implemented with a conditional Dirichlet condition via a Prescribed Displacement/Rotation node in the Beam interface. Beam with a deflection, controlling support and distributed load. Settings prescribing that the beam should stop at a deflection of 2 cm.The ysis shows the expected behavior. At lower loads, the deflection shape is symmetric, whereas at higher load levels, the point on the beam where the extra support is located will stop moving. At the final load level, the beam will even undergo a change of sign in the curvature. This is visible in the deformation plot, but it is shown more clearly in a bending moment graph. The beam displacement at the support point stops at 2 cm. Bending moment along the beam for various load levels.The approach highlighted here is rather crude and the iterative solution may not have good convergence properties. A more stable implementation is to use a highly nonlinear spring at the support point, so that the reaction force is a continuously differentiable function of the displacement. This is actually similar to how penalty contact is implemented in the Solid Mechanics interface.COMSOL Multiphysics gives you access to very powerful mechanisms for prescribing nonstandard boundary conditions. Today, we have provided you with a few examples of what you can do with these conditions.For those who are interested in yzing a model with a traveling load, take a look at the Traveling Load tutorial model, available in our Application Gallery. If you have additional questions on how to prescribe nonstandard boundary conditions within your own modeling processes, please contact us today.Hi HenrikThanks for an excellent BLOG on very useful “Tricks” to handle complex physics cases in COMSOL. I was familiar with some of them, but had not really thought of all other combinations you are showing here You have such a powerful tool and there is so many engineers and scientist not even knowing about COMSOL, but too bad for them Sincerely IvarThank you! it is such useful BLOG.Hi Henrik,I have a question about solution dependent boundary condition (BC). I am solving a Poisson-continuity fully coupled equation involving two variables, which are electric potential (phi) and space charge density (rho). I am able to obtain converged results when using Dirichlet BC for phi. However, when I change the Dirichlet BC to a variable BC, which depends on the variable rho and gradient of phi, the stationary solver does not lead to convergency. I guess the model is not well posed after using the variable BC.I have no prior experience in implementing a variable BC, so any suggestion will be greatly appreciated. EricEric,Please send questions of this type to support, or post them on the user forum. I think it would be necessary to attach a model to get help in either case.Regards, Henrikthank you, it was so useful for me.In electric current, there is no term like fixed current for Neumann boundary conditionSagar – The name of the Neumann boundary condition in the Electric Currents interface is ‘Normal Current Density’. The table in the beginning of the blog post is more conceptual.Regards, Henrikthank you, it was so important for meIndiaBrazilChinaDenmarkFinlandFranceGermanyIndiaItalyJapanNetherlandsNorwayPortugalRussiaSpainSwedenSwitzerlandUnited Kingdom