Discussion Forum

Posted: 1 decade ago Feb 23, 2012, 3:58 p.m. EST

Hi looks strange (I'm not by my Comsol WS so I cannot take a look) but are you in frequency domain / harmonic mode ? in which cas nJ might be an rms value and Jx a pk value ? Check the "equation view" how these are defined internally in COSMOL -- Good luck Ivar

Posted: 1 decade ago Feb 23, 2012, 7:19 p.m. EST

Hi Ivar, The simulation was in Time Dependent mode, but the result is exactly the same in Stationary mode too. ece.Jx is defined as ece.Jex+ece.Jix which are defined as 0 and ece.sigmaxx*ece.cucn1.input.Ex+ece.sigmaxy*ece.cucn1.input.Ey+ece.sigmaxz*ece.cucn1.input.Ez respectively. Where ece.sigmaxy and ece.sigmaxz are zero, so we have it reduce to ece.sigmaxx*ece.cucn1.input.Ex Decomposing further: = 1 * model.input.E1 And I can't find where model.input.E1 is defined. So we have ece.Jx = model.input.E1 for my model. Unfortunately this doesn't mean anything to me. I can't seem to find where ece.nJ is defined. In Equation View under Normal Current Density it says "Name: ece.nJ Expression: -eci.nJ". That's really the constraint that I specified, not the definition of nJ. I've been through all of the Equation Views I could find, and there is nothing that links Jx to nJ explicitly. That Jx = -nJ at the boundary, was just my common sense assumption.

Posted: 1 decade ago Feb 24, 2012, 2:16 a.m. EST

Hi first a trick to make life easier (at least i find so) you can add Definition variables, select domain 1 and write V = Ve, then add another Variable node and selct domain "2" where you write V = Vi In this way you can access the field V(x,t) by using the lettre "V" alone also for your plots, COSMOl will sort out which fomula to use and when Then to be strictly pedant, your initial conditions shoul read "9[V]*(1- x/2[m])" for both EC's where 2[m] is the total length of BOTH domains (lines, edges) I agree that it's not "obvious" to think of a "normal" current when one use a 1D physics and mean a point (COMSOL developers: what is the normal to a point ?) Then in the doc/help it's written "The normal current density is positive when the current flows inward toward the edge." But in 1D what is an edge ? It's a domain. COMSOL developers pls correct "edge" to "boundary" and PLS be strict with your naming of "entities" so that the help remains valid for 1, 2 and 3D ;) This said your model should be correct, but for one thing, it's not bidirectional you impose thevoltage on "eci" on the left boundary but you do not impose any voltage on the right boundary of "ece" just as you impose a normal current to the right boundary of "ece" but nothing on the left boundary on the current One confusion here is that COMSOL, in union mode transforms your two domains (two lines) into a section with three boundary (points) where the middle point is overlapping, but it does NOT impose any continuity as you have 2 different physics, one for each domain. Even if you use Assembly Finish mode, you will dedouble the central point into two Ids, with the same coordinates, and it will propose an default identity pair, but if you select it and apply continuity, this continuity applies ONLY within one physics. So you need to apply continuity by hand, the way you have done it, BUT in both directions. Now, if you define both current density AND voltage on both sides of the interface you will get a different result (check carefully what COMSOl does to your BCs), than if you only define a common voltage from both sides, or common current only. But I leave this to you to sort out ;) Think it over with respect to your PDE order and the minimum required BC settings required Nice model to learn :) -- Good luck Ivar

Posted: 1 decade ago Feb 26, 2012, 11:42 p.m. EST

Hi Ivar, Thank you for your detailed reply. Guess what I discovered! It seems that [b]constraining the Electric Potential also constrains its Jacobian[/b]!!! Now, what I thought is that it will be a Dirichlet BC: just constrains the value. But it seems to be a Dirichlet+Neumann BC, both at the same time: value and Jacobian. Not something I expected. I have the following clues to support this: 1. Have a look at the model attached. It is the same one as before, except I have disabled the Normal Current Density BC, and changed the Electric Potential BC from Ve (which just gave continuity) to -0.5*Ve. Now what happens: there is no current continuity across the interface, but there is Jacobian continuity alright. 2. I've tried the N ways of coupling 2 physics: a) Form Union. b) Form Assembly + Identity Pair. c) Form Assembly + Linear Extrusion. d) Form Assembly + Identity Mapping. Every time it gave the same result. 3. Switching from "Bidirectional" to "Unidirectional" constraint gives me what I want (a pure Dirichlet BC, not one of these Dirichlet+Neumann BCs). 4. Using Bidirectional, but now specifying -0.5*nojac(Ve) also gives me what I want. In other words it's telling it: don't consider the Jacobian. SOLUTION SUMMARY: I should have used "Unidirectional" instead of "Bidirectional, symmetric" in constraint settings under "Electric Potential". The reason it's called "bidirectional" is, I think, due to the algebraic equation behind it, see section "Ideal vs Non-ideal constraints" in 3.5a modeling.pdf OR "Using Weak Constraints" - > "Bidirectional and Unidirectional Constraints" in 4.2a User's Guide. The content there is too high-level for me to understand (and no references), but I think it's a good rule of thumb to say that Bidirectional couples the value of the dependent variable AND of its Jacobian, while Unidirectional couples only the value. A real useful example to understand all this is in modeling.pdf for version 3.5a, on page 354: it's called "Example — Coupling Variables and Boundary Constraints". A model is built, then ideal and non-ideal constraints are used, giving different solutions. This is explained right at the very end of that section, under the heading "Weak Non-Ideal Constraints". I was thinking that, logically, the documentation will get better and better as the version number increases, but it seems that this whole example is missing from the documentation for 4.2a, so it's worth checking the previous documentation! With regards to 1D and 2D models, I think about them as still being 3D, except that they are redundant in (Y and Z dimensions for 1D) and (Z dimension for 2D): either infinitely long in those dimensions, or, more physically, in the middle of very very long extents along those dimensions. So in 1D we have an interval, but it represents a large plate with thickness == length of interval. Also, thanks for the timesaver tips. Knowing all these things saves hours of frustration and drudgery for me and others. Evgeni

Attachments:

• nJ_Jx_discrepancy2.mph

Posted: 1 decade ago Feb 27, 2012, 5:56 a.m. EST

Hi I believe I was partly fooled in my previous reply, by the initial conditions you were using. First of all, clean up your model you have leftovers of solution 1&2 that are getting confusing, in fact I would say restart from scratch as after some time with changes, I suspect that a model can become currupt. Now, once I have reset both initial conditions to = 0[V] , I get something more according to my first estimates, by disabling the current denisty node in ece, anabling bidirectional Electric potential to Ve in ece, and keeping ONLY a GND on the eci side both currents and voltages are coherent If you engage a monodirectional Electric Potential to Vi on ece point "2" your eci remains unchanged set to 0[V] If you check the equations you will see the Constraint force change from "test(ece.V0-Ve)" to "-test(Ve)" which means Ve=0 Which you see clearly if you look at the correct Ve,Vi,ece.Jx, eci.Jx graphs This works in symmetry if you move the Electric Poitential from ece to eci (still on common point 2) And you get a coherehnt resply if you leave an electric potential, bidirectional on both ece,eci even if 1 is probybly to much -- Good luck Ivar

Posted: 1 decade ago Feb 27, 2012, 7:44 p.m. EST

Hi Ivar, With your help, I think I'm starting to understand. The problem is that I don't understand what is meant by [b]Constraint force[/b]. In Equation View for Electric Potential. Under "Constraints", there are:     Constraint: eci.V0-Vi     Constraint force: test(eci.V0-Vi)         where eci.V0 is defined as Ve This is for bidirectional, symmetric, not using weak constraints. If I change it to unidirectional:     Constraint: eci.V0-Vi     Constraint force: -test(Vi) Somehow the "Constraint force" test(eci.V0-Vi) results in a coupled Jacobian. If we check "Use weak constraints", then there is a Lagrange multiplier Vi_lm introduced, and two weak expressions.     For bidirectional:         (eci.V0-Vi)*test(-Vi_lm)         -test(eci.V0-Vi)*Vi_lm     For unidirectional:         (eci.V0-Vi)*test(-Vi_lm)         test(Vi)*Vi_lm My question is: how does the bidirectional constraint (either weak or not) result in a coupled Jacobian? The way I know it, whatever multiplies the test functions becomes zero. But I don't know what test(eci.V0-Vi) is. I assumed that Comsol was using the Galerkin method, so a test function of a dependent variable evaluated to the shape function on a given element. What does it mean to have a test function of an expression? I found a small clue in the "Constraint Forces for PDEs" topic in Help. It says for user-defined constraint types: "...enter the expression in the Constraint force expression field. In this field specify the constraint force Jacobian using the test operator." From this I learn that: • Constraint force is a Jacobian. • The test(..) operator can be used to specify a Jacobian. Looking at the definition of the test(..) operator, it seems that test(Ve-Vi) expands to:     test(Ve) * ∂/∂Ve (Ve-Vi) + test(Vi) * ∂/∂Vi (Ve-Vi) Which almost looks like something to do with a Jacobian, except it's not, because they are the derivatives with respect to Ve and Vi, not spatial variables. It seems that the answer is real close, but here is where I'm stuck. Thanks, Evgeni

Posted: 1 decade ago Feb 28, 2012, 3:53 a.m. EST

Hi I find it also delicate to understand the details of the weak form, even after a few years, And answering in detail will take me too long just now. For me the Jacobian is extracted from the "deformation mapping" of the model and is NOT directly related to the test() operators or the weak form, apart that adding a weak form influences the deformation mapping hence the Jacobian (in the opposite direction from how I read you) What I can propose is that you ask by your library for the book of Zimmermann "Multiphysics Modelling with Finite Elements Method" unfortunately there it is written with examples for 3.5 and this requires some tweaking to get it into v4 but all the stuff is therein. It's also taylored for chemistry, but the first chapters are generic and its well written. Another good book is Lennart Edsberg "Introduction to computation and Medelling for Differential Equations" You can find more on http://www.comsol.eu/support/books/ Still another book, not mentioned by COMSOL, it's brand new, and I like it for its concise and precise approach, also very close to COMSOL notations, and with a lot on the material versus spatial frame + tensor fields and mapping, how to understand the "Jacobian", is "Continuum Mechanics and Thermodynamics", CUP, 2012 by E.B.Tadmor and R.E.Miller -- Good luck Ivar

Posted: 1 decade ago Feb 28, 2012, 11:35 p.m. EST

Hi Ivar, I've looked through a few generic FEM textbooks, but couldn't find anywhere the use of Lagrange multipliers for weak form coupling in the same way as Comsol uses them. (J. N. Reddy has an example with a Lagrange multiplier constraint, but it's not quite the same.) A Comsol-specific book is the next best source to check, I guess. Trying to get to the bottom of concrete issues is a great way to learn :) I'll post the answer here if I don't give up before finding it. Thanks for your help, Evgeni

Posted: 1 decade ago Feb 29, 2012, 12:48 a.m. EST

Hi I ahave another suggestion for you: Check how COMSOL does it Take your model, save it with a new name, in the geometry Finish select Assembly and check create identity pair (this will create the selection for you). In assembly mode you have 2 domains (lines) and 4 points (not 3 as with Union) as the middle point is dedoubled, one "up" and one "down" (left and right for your case). Then get rid of the second EC the ECI, just delete it and all domains to the first EC. Then disable the potential on the middle point, and replace it by an Continuity PAIR BC condition and select the identity pair, and then take a look at the equations COMSOL is using. THe best is to do a Edit Clear All Solution and Clear All mesh then File Reset Model Save and then execute to get the results. -- Good luck Ivar