Note: This discussion is about an older version of the COMSOL Multiphysics^{®} software. The information provided may be out of date.

**Discussion Closed** This discussion was created more than 6 months ago and has been closed. To start a new discussion with a link back to this one, click here.

## Maxwell stress tensor integration

Posted May 24, 2010, 8:46 PM EDT AC/DC & Electromagnetics, RF & Microwave Engineering Version 4.3 38 Replies

I'm still trying to compare the electrostatic force calculated by the expression (F=0.5*epsilon*V^2*A/g^2) with the result coming by Maxwell stress tensor in COMSOL.

I used the suggestion of Ivar KJELBERG (in the Model Exchange part): mapped mesh, axial symmetry. But still the results don't match (3 orders of magnitude).

Once I set the variable Fes in the force table, I get the result from Post processing-->data display-->global-->Fes_z

If I try to integrate the stress tensor through the same bonderies (surrounding the Fes box) I get a N/m value. What is strange is that the relation between these two forces is the perimeter of the plates (in my cases 2*radius*pi.greco since I'm using circular plates with the axial symmetry).

Is there anyone that already tried to computer the electrostatic force by COMSOL with good results?

Has anyone an idea about the behavior of the boundary integration of the Maxwell stress tensor?

Thank you very much to everyone will give some suggestions

Riccardo

Were you able to find the solution to your problem? I am having a similar problem with differences between my Maxwell's stress tensor integration and my analytical calculations for electrostatic force.

My model is attached; am I doing anything clearly wrong?

Thanks,

Kristen

Attachments:

Indeed maxwell Stress tensor calcuations are slightly tricky, as they are besd on a few hypothesis and impliesintegration of steep gradients (often).

1) check that you have "air" or vacuum all around your part

2) that the (dense) mesh is symmetric (as far as possile aroud the integration edges/boundaries)

3) check edge by edge (boundary by boundary) and look at the results w.r.t ie.e symmetries, does it look coherent

In case of doubts, use the other variatoal methods

And its worth to read carefully the 2-3 examples in the doc

Good luck

Have fun Comsoling

Ivar

What version are you using?

There is a confimed bug in the Maxwell stress tensor integration for axial symmetry in version 4.0.

Due to the following change (copied from the release notes):

"Axisymmetric Models

In version 3.5a equations in axisymmetric application modes use the independent variable for the radius, r, to account for axisymmetry. In version 4.0, the equations are compensated by using the factor 2?r.

If you have manually multiplied expressions by 2? for a version 3.5a model, please note that these may be incorrect when the model is opened in version 4.0."

Johan

I cannot open your file, your version is more recent, anyway I'm using this configuration. I send you this guide I found on internet. So, since the capacitance value matches I assume that also the Force value should match, at least the order of magnitude, instead it isn't.

I've also tried to plot the values on excel in order to find the law behind the force computed from comsol..noway.

Riccardo

Attachments:

I used all the suggestions you gave but still no match.

The force on the upper plate (force integrated on the boundary) has a different value from the one on the lower plate.

Should I instead expect the same value?

in that case, is a problem related with the mesh?

If you have time just to take a look to the file I will be really grateful.

Riccardo

Attachments:

version 3.5, so is not the bug, damn.

thank you anyway for the suggestion!

Riccardo

I am running into the same problem- the Maxwell stress tensor has a different value for top and bottom plates. This is for the most simplified case of a parallel plate capacitor, so I am quite surprised to see this result. I have air surrounding both plates and refined the mesh. Can you offer any suggestions as to why this might happen?

Thanks,

Kristen

The 2*pi*r is not really a bug, it's a decided change in the notations of the physics for V4 and its announced in the release notes, but I agree the reading in of a 3.5 file into 4 should check this

Ivar

Yes, it is a change in notation but the Maxwell Stress Tensor integration that is used by the "Force Calculation"-node has not been updated to reflect this change. This requires you to manually integrate the stress tensor to get the correct result. Perhaps the attached example files can illustrate what I mean.

Regards,

Johan

I need some time to study the files, but if you are sure about a bug, pls send the info to "Support", these people are efficient to correct whatever is wrong, and you might well find it corrected in V4.0a ;)

Have fun Comsoling

Ivar

Yes, I am sure that it will soon be fixed. I informed the support about the bug as soon as I found it and they confirmed it and are working on correcting it.

Johan

I have had a very quick look at "panic", and I have a few comments:

use the "Compute surface integrals for axisymmetric modes" "on" in 2D axi to add the 2*pi*r to the integration and you get total forces in Newton and not N/m (pls note also the change of notation how to use the 2*pi, but not of the r between V3.5 and V4)

Then OK you get different values for your line integrals on the top and bottom boundaries, and different values if you take the interiour forces (due to the er <> 1 and the way the maxwell stress tensor is calculated).

My only explanation (the quick one) is that you are getting effects due to the E field edge concentration at the external vertex of your electrode boundary. The electric field expells the electrons to concentrate on the (shrp) external rim of your electrode. You see that the electric field is concentrating at the vertex making a numerical singularity (always avoid sharp edges in HV they tend to create arcs!).

I have tried t cut your high "er" volume by a line close to the edge and you see that the forces are closer equal on top and bottom edges up to the region where the non grounded electrode starts to see high E gradients.

This is my only explanation, have you tried the other methods proposed ? for such simple models the variation principle should be easy to implement and should give good results. Checkthe ACDC doc

Have fun Comsoling

Ivar

Attachments:

I was thinking maybe the problem was in defining moving mesh properties, but after reading this thread it seems like the problem is with linking two physics using the maxwell stress tensor in version 4.0. Because I able to this simulation in V 3.5a before.

I am attaching my file here. Maybe some body could give me some insight on what might be wrong. I am fairly new with comsol, specially v 4.0, so any help would be really appreciated.

Attachments:

There are definitely a few mistakes...

--------Solid Mechanics physics domain--------

**The air domain should not be active here in the structural mechanics domain. Delete them.

**Add the bottom boundary of the bottom plate as a fixed constraint (or else it will be interpreted as being free)

**The ONLY boundary that should be added as an edge load is the bottom boundary of the top plate (unless you are trying to place the electrode at the top of the plate, like I currently am having trouble with)

**The equation for the electric potential force to be used in the structural domain is not Fes.nTx.emes (as it was in 3.5a)... For 4.0, the "Maxwell downward Surface Stress Tensor" to be used is denoted as "es.nTx_Fes".

----------Electrostatics physics domain---------

**Delete the "Charge Conservation 2" menu item.

**Zero charge should probably be applied even to the fixed side edges of the bottom plate

**The "Force Calculation 1" item may not be needed, since we already defined the Maxwell surface stress tensor.

**The Ground item may only need one boundary; the upper surface of the bottom plate.

Try adding in an ALE mesh.

Hope this helps a little!!

You were adding the force in which seems right. However, the naming convention is definitely different in v4.0

Thank you fro your help. I got the simulation to work now. I think making the "air" inactive in the solid mechanics physics was important to get the ALE working right. If its active, I only see deformation in the air , not in the cantilever itself.

But at pull-in voltage, approx. 16.7 V it runs into trouble. It cannot finish computation. I think its due to the gap becoming zero. Also I believe, at pull -in voltage there should be a sharp drop in displacement right away. But I see that it starts to bend faster after 16 V. but doesn't collapse. I was wondering whether defining a contact pair between the two surfaces would help to get a better picture of the cantilever collapsing at pull-in voltage. Has anyone done any modeling of contact pair with electrostatics? Is there any example I can look at?

Thanks again!

Have you implemented the contact pair scheme? Has it solved your problem near pull-in voltage? IÂ´m having a similar problem, when the parameter V_in gets near to the pull in voltage, the simulation stops with an error.

IÂ´m currently working on a parallel plate accelerometer model and I was playing around with it by increasing the voltage of the botton plate.

I used a 1d-graph to compare elastic force Fel (k*vertical displacement w) and eletrostatic force Fes (emes.nTz_Fes (N/m2) * Plate Area). But Fes is still far away from Fel when the simulation stopped working. I was expecting that both values should be at least close to each other near the Pull in.

Am I plotting wrong values?

Thanks!

My research is in a similar area involving electrostatics with a MEMS device for an ultrasonic transducer. From my understanding, the ALE Moving Mesh stops at the pull-in voltage. If you really need to see the membrane deflect all the way onto the bottom substrate, then I doubt this information I provided will help you visualize that. Then again, I don't see why it would be necessary (unless your specific project requires this result).

But you originally said that you were trying to recreate the "2D MEMS ALE Beam" from v3.5a to v4.0. The v3.5a also does not model the deflection past the pull-in voltage (because of the fact that the ALE Moving Mesh doesn't calculate past the pull-in, as noted in the COMSOL 2D ALE Beam modeling guide). Besides, just past pull-in all that happens is that the membrane/beam collapses onto the bottom substrate (i.e. the unstable region).

Salwa, you also stated that you think the model isn't solving because the gap is becoming zero. I'm do not remember how you were trying to model your device, but an equation that I have been going by for the initial gap height between the top and bottom electrodes has been:

g.0 = (d.membrane / ?.membrane) + (d.vacuum) + (d.bottom/ ?.membrane)

where (d.membrane) and (?.membrane) are the thickness of the domain and the relative permittivity of the material, respectively. The vacuum height is the height between the beam and the substrate (in your case, maybe it isn't a vacuum).

Depending on your setup of materials and thicknesses with respect to the equation above, it may be possible that your device is deflecting at a height more than the vacuum height. For example, if you are placing the V_in electrode at the top of your membrane/beam and the ground at the bottom of an insulation material, both enclosing a vacuum layer in between, then based on your thicknesses and material properties, it may be possible that your device collapses at a height greater than the vacuum layer, which obviously doesn't make sense. Would a contact pair solve this issue? I'm guessing yes but have not had a need to investigate.

Again, as far as actually visualizing the collapsed structure I doubt this helps. But hopefully this information was a bit informative.

Best,

Sadat

I don't know if this helps much, but I had a trouble with my Fes.nTx.emes force from the electrostatic domain in my model.

It turned out for me that the placement of the force on the boundary was not fully correct. I added my Fes force to both the top and bottom of the membrane, and things worked out for me. Let me know if this seems right for you..

Thank you very much for your reply.

I am interested in the behavior of the beam after it collapses. Practically it is seen that past pull-in voltage and collapse the cantilever beams tend to flatten on the bottom surface and then bends upwards.

Here's what I did with the simulation. I put a thin layer ( 1 um) of insulating Silica glass above the grounded bottom pad. So when the beam collapses past pull in voltage the electrical potential induced cantilever beam doesn't come in contact with the ground. This kind of worked as it was able to calculate past pull-in voltage showing huge displacement of the beam ( which is impossible, as it showed deflection (20um) past the gap (4um).)

Next I introduced a contact pair, where the top surface of the silica glass and bottom surface of the beam ( V induced) are the destination and source pair. This simulation result could calculate up to contact point showing 4um displacement exactly but not beyond that. It may be due to defining the contact pair's manual scaling options. So we'll have to play a little bit to find the suitable values of manual scaling for each variable.

Talking about Maxwell tensor, anybody know how I could transfer the maxwell tensor from the boundry layer to the subdomain.

The main reason is I need to couple the physics of acoutics and stuctural and it seem that from the loudspeaker example they use the sub-domain as the loading to the system, and have the acoustics/stuctural coupled through acceleration and pressure.

The main thing here is, I have the electrostatic force FES - Fes_nTx_emes or Fes_nTr_emes for axial symm for my force, but the tensor in mainly for the boundry. How do we transfer them to the subdomain.

Currently, I am using the boundry integration to get the full results and multiply by area and divide by volume. That does not seem to get the results I expected.

Anybody, have performed a good tranfer for the force.

Or should I go though to the first principle and enter F = 0.5 * Epsilon * V^2 * A / g^2 directly to the force, and will that work?

the maxwell stress tensor is only defined on the material - air boundary so it applies to the boundary, and is "tranported" to the domain via the structural model.

Note that if you integrate the value first you get an average value that you distribute, and not the detailed local value, this might change the moments applied

--

Good luck

Ivar

Thanks for your insight, now I understand what I did wrong on the boundry - subdomain transfer. If that is the case, what is the best way that I could put a localised load within the subdomain rather than through the boundary. I tried to search whether anybody have done this before, but it seem that from my understanding, the force need to be inserted through 1st principle. Any way, that we could transfer the localised force between the air-solid boundary to the solid domain directly?

I have recheck my example on the loudspeaker but it seems that they do not use any moving mesh in that sense, so I presume it is a bit different to me then.

Still searching for a suitable way. Anyway thanks for you help again Ivar.

Syamsul

Read more: Maxwell stress tensor integration - COMSOL www.comsol.com/community/forums/general/thread/5840/#ixzz0z0HSrpPk

__________________

moviesonlineworld.com

I'm not sure if I completely understand your problem, but from what I've noticed (and have been informed by COMSOL support), much of the Maxwell Stress Tensor contribution occurs on the bottom-portion of the plate.

For instance, if you constrain the bottom plate, and allow the top plate to fluctuate, stress placed at the bottom portion of the top plate will contribute more to the deflection than the upper portion of the plate.

Again, I'm not sure if I understood you, but this has been my experience. Hope that bit helps somewhat!

the Maxwell stress tensor integration is very sensitive to mesh to gradient ratis, you need a rather fine and if possible regular mesh t get low errors, one should also look at the local values and its fluctuations, not only the integrated final value, as this gives a clue on the accuracy of the tensor derivation.

Have you tested the results with different mesh densities, or looked at an arrow plot ofthe maxwell stress tensors ? Do it with a raw mesh overlay, so you see the true raw wireframe mesh grid density and not the "smoothed" views

--

Good luck

Ivar

I realize this is an old thread but I am now trying to do some Maxwell stress tensor calculations in COMSOL 4.3 and I am wondering if anyone knows if these issues still arise. As a simple test case, I've modeled two point charges and I'm trying to integrate to find the force one exerts on the other. For two charges +/- q, separated by a distance 2d, one can find an analytical solution for the force on a hypothetical cube surrounding one of the charges. That should be

integral(T n)dA=-q^2/(4*pi*eps0)/(2d)^2

(dl.dropbox.com/u/3102425/maxwell.png)

For my model, I have q=1e-10 C and d=1um. The resultant force should then be about 22N

However, when I calculate the integral of the M.S.T. in post processing, I get values between ~40N and >100N, depending on my mesh size. Does anyone have any insight into why this might be?

Thanks,

Alisha

Attachments:

one guess, you "box" is still influencing the filed lines, why not try a sphere ? And perhaps some INF domains

And get a dense symmetric mesh, if you compare the up and down Tex values they are rather different, but in your symmetric model I believe they should be identical

Also read carefully the entries in the knowledge base on flux calculation and Maxwell tensor integration

--

Good luck

Ivar

The sphere is a good suggestion, I don't know why I didn't think of that. With a spherical system and some finer meshing, I get values closer to what I except. However, there are still two issues: (1) the up/down integrations are slightly different, as you mention and (2) the solution does not get more accurate as I decrease mesh size. Instead, there seems to be a "sweet spot." That's not a problem in this example since I can check the integration over the "right" answer, but I obviously want to extend this to problems I don't know the answer to. Do you have any more suggestions regarding the mesh? You mentioned symmetry, although the default mesh does look symmetric, since it's a symmetric geometry.

Regarding infinite elements, I did try that as well, although I've never used those before so this is just a guess at how to set it up. The results look even more mysterious, however...

Thanks again for your help and for your extremely quick response.

Regards,

Alisha

Attachments:

Maxwell tensor integration is tricky and very mesh dependent.

One thing to try is to cut your sphere in two, symmetrically and copy-mirror the mesh to have exactly the same mesh on both sides (but I'm not sure this was implemented in v4.0 ;)

On the other side with such a symmetric model its even better to consider only 1/2 and mirror / anti mirror the physics as required

--

Good luck

Ivar

keep also in mind that Point charges give a slightly mesh-dependent electric potential as well. You might want to try with equivalent spherically-symmetric charge distributions if you desire a more accurate result.

--

Andrea Ferrario

Electromagnetics Group

COMSOL AB

I have a couple of questions related to the Maxwell stress tensor:

1. What is the difference between the upward and the downward stress tensor? Is it the stress tensor evaluated just outside and inside a boundary? I didn't find these definitions in the reference guide.

2. Is it possible to evaluate the stress tensor inside a material (i.e. not on a boundary)? Or is that only possible by calculating it manually (in MATLAB) by exporting all the fields?

Thanks!

Hi,Ivar Kjelberg:

I built a 2D micro-ring resonator model which stated in one publicated paper.

I check all the parameters I set are same as the paper mentioned. But the problem is , even in the ring, no field distributed. But when I further decrease the gap from 540nm-100nm, it occurs. however the extracted S21 is unphysical, I can not get series resonance spectrum.

I attached my mph file. Do someone can help me to find the reason.

Anyhelp is highly appreciated.

Attachments:

I have been experiencing similar difficulties in Maxwell Stress Tensor Integration. I am simulating the force on a conductive droplet (statically modelled as a solid) in air between two plates at multiple positions (digital microfluidics) and am investigating the horizontal component of the force. I notice some numerical noise in my Force Vs Position graphs (I know analytically that there should not be). It took Ivar's tip and created a uniform swept quad mesh along my surface of interest to help reduce the noise and am using a very fine element setting (~150 000 elements) but the noise only slightly decreases with an increase in element size.

I believe the issue is because the Z component of my force is very large compared to the X component (100x difference). However, I am only interested in the X component and am thinking that this could be creating numerical issues.

Any help on dealing with Z >> X issue would be greatly appreciated and I will make sure to inform the group of my findings as this seems to be a hot topic. I have also tried using both the Force Calculation Node, and manually integrating the X components of the stress tensor over the surface but both have the numerical noise (I am not completely confident that I have done the manual surface integral correctly however).

Happy New Year!!!

I am going to derive DEP forces over an spherical particle by using Maxwell stress tensor but I have no idea about integration of maxwell stress tensor. could anyone help me?

Best,

could anyone help me?

Best Regards !

Have you gotten the answer?

Do you know how to calculate the Maxwell stress tensor in MATLAB? I have no any idea ...

Do you have any referrences ?

Thanks !

Note that while COMSOL employees may participate in the discussion forum, COMSOL^{®} software users who are on-subscription should submit their questions via the Support Center for a more comprehensive response from the Technical Support team.